Ejercicio – Página 2

Números autodescriptivos

PorJosé A. Alonso 1-julio-202221-junio-2022

Un número n es autodescriptivo cuando para cada posición k de n (empezando a contar las posiciones a partir de 0), el dígito en la posición k es igual al número de veces que ocurre k en n. Por ejemplo, 1210 es autodescriptivo porque tiene 1 dígito igual a «0», 2 dígitos iguales a «1», 1 dígito igual a «2» y ningún dígito igual a «3».

Definir la función

   autodescriptivo :: Integer -> Bool

1	autodescriptivo :: Integer -> Bool

tal que (autodescriptivo n) se verifica si n es autodescriptivo. Por ejemplo,

   λ> autodescriptivo 1210
   True
   λ> [x | x <- [1..100000], autodescriptivo x]
   [1210,2020,21200]
   λ> autodescriptivo 9210000001000
   True

λ> autodescriptivo 1210

True

λ> [x | x <- [1..100000], autodescriptivo x]

[1210,2020,21200]

λ> autodescriptivo 9210000001000

True

Soluciones

import Data.Char (digitToInt)
import Data.List (genericLength)
import Test.QuickCheck

-- 1ª solución
-- ===========

autodescriptivo1 :: Integer -> Bool
autodescriptivo1 n = autodescriptiva (digitos n)

-- (digitos n) es la lista de los dígitos de n. Por ejemplo.
--    digitos 325 == [3,2,5]
digitos :: Integer -> [Integer]
digitos n = [read [d] | d <- show n]

-- (autodescriptiva ns) se verifica si la lista de dígitos ns es
-- autodescriptiva; es decir, si para cada posición k de ns
-- (empezando a contar las posiciones a partir de 0), el dígito en la
-- posición k es igual al número de veces que ocurre k en ns. Por
-- ejemplo, 
--    autodescriptiva [1,2,1,0] == True
--    autodescriptiva [1,2,1,1] == False
autodescriptiva :: [Integer] -> Bool
autodescriptiva ns = 
  and [x == ocurrencias k ns | (k,x) <- zip [0..] ns]

-- (ocurrencias x ys) es el número de veces que ocurre x en ys. Por
-- ejemplo, 
--    ocurrencias 1 [1,2,1,0,1] == 3
ocurrencias :: Integer -> [Integer] -> Integer
ocurrencias x ys = genericLength (filter (==x) ys)

-- 2ª solución
-- ===========

autodescriptivo2 :: Integer -> Bool
autodescriptivo2 n =
  and (zipWith (==) (map digitToInt xs)
                    [length (filter (==c) xs) |  c <- ['0'..'9']])
  where xs = show n

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_autodescriptivo :: Positive Integer -> Bool
prop_autodescriptivo (Positive n) =
  autodescriptivo1 n == autodescriptivo2 n

-- La comprobación es
--    λ> quickCheck prop_autodescriptivo
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> [x | x <- [1..3*10^5], autodescriptivo1 x]
--    [1210,2020,21200]
--    (2.59 secs, 6,560,244,696 bytes)
--    λ> [x | x <- [1..3*10^5], autodescriptivo2 x]
--    [1210,2020,21200]
--    (0.67 secs, 425,262,848 bytes)

import Data.Char (digitToInt)

import Data.List (genericLength)

import Test.QuickCheck

-- 1ª solución

-- ===========

autodescriptivo1 :: Integer -> Bool

autodescriptivo1 n = autodescriptiva (digitos n)

-- (digitos n) es la lista de los dígitos de n. Por ejemplo.

-- digitos 325 == [3,2,5]

digitos :: Integer -> [Integer]

digitos n = [read [d] | d <- show n]

-- (autodescriptiva ns) se verifica si la lista de dígitos ns es

-- autodescriptiva; es decir, si para cada posición k de ns

-- (empezando a contar las posiciones a partir de 0), el dígito en la

-- posición k es igual al número de veces que ocurre k en ns. Por

-- ejemplo,

-- autodescriptiva [1,2,1,0] == True

-- autodescriptiva [1,2,1,1] == False

autodescriptiva :: [Integer] -> Bool

autodescriptiva ns =

and [x == ocurrencias k ns | (k,x) <- zip [0..] ns]

-- (ocurrencias x ys) es el número de veces que ocurre x en ys. Por

-- ejemplo,

-- ocurrencias 1 [1,2,1,0,1] == 3

ocurrencias :: Integer -> [Integer] -> Integer

ocurrencias x ys = genericLength (filter (==x) ys)

-- 2ª solución

-- ===========

autodescriptivo2 :: Integer -> Bool

autodescriptivo2 n =

and (zipWith (==) (map digitToInt xs)

[length (filter (==c) xs) | c <- ['0'..'9']])

where xs = show n

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_autodescriptivo :: Positive Integer -> Bool

prop_autodescriptivo (Positive n) =

autodescriptivo1 n == autodescriptivo2 n

-- La comprobación es

-- λ> quickCheck prop_autodescriptivo

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> [x | x <- [1..3*10^5], autodescriptivo1 x]

-- [1210,2020,21200]

-- (2.59 secs, 6,560,244,696 bytes)

-- λ> [x | x <- [1..3*10^5], autodescriptivo2 x]

-- [1210,2020,21200]

-- (0.67 secs, 425,262,848 bytes)

El código se encuentra en GitHub.

Aproximación del número pi

PorJosé A. Alonso 30-junio-20222-julio-2022

Una forma de aproximar el número π es usando la siguiente igualdad:

    π         1     1·2     1·2·3     1·2·3·4     
   --- = 1 + --- + ----- + ------- + --------- + ....
    2         3     3·5     3·5·7     3·5·7·9

π 1 1·2 1·2·3 1·2·3·4

--- = 1 + --- + ----- + ------- + --------- + ....

2 3 3·5 3·5·7 3·5·7·9

Es decir, la serie cuyo término general n-ésimo es el cociente entre el producto de los primeros n números y los primeros n números impares:

               Π i   
   s(n) =  -----------
            Π (2*i+1)

Π i

s(n) = -----------

Π (2*i+1)

Definir la función

   aproximaPi :: Integer -> Double

1	aproximaPi :: Integer -> Double

tal que (aproximaPi n) es la aproximación del número π calculada con la serie anterior hasta el término n-ésimo. Por ejemplo,

   aproximaPi 10     == 3.1411060206
   aproximaPi 20     == 3.1415922987403397
   aproximaPi 30     == 3.1415926533011596
   aproximaPi 40     == 3.1415926535895466
   aproximaPi 50     == 3.141592653589793
   aproximaPi (10^4) == 3.141592653589793
   pi                == 3.141592653589793

aproximaPi 10 == 3.1411060206

aproximaPi 20 == 3.1415922987403397

aproximaPi 30 == 3.1415926533011596

aproximaPi 40 == 3.1415926535895466

aproximaPi 50 == 3.141592653589793

aproximaPi (10^4) == 3.141592653589793

pi == 3.141592653589793

Soluciones

import Data.Ratio ((%))
import Data.List (genericTake)
import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

-- 1ª solución
-- ===========

aproximaPi1 :: Integer -> Double
aproximaPi1 n = 
  fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

-- 2ª solución
-- ===========

aproximaPi2 :: Integer -> Double
aproximaPi2 0 = 2
aproximaPi2 n = 
  aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

-- 3ª solución
-- ===========

aproximaPi3 :: Integer -> Double
aproximaPi3 n = 
  fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,
--    λ> take 10 numeradores
--    [1,2,6,24,120,720,5040,40320,362880,3628800]
numeradores :: [Integer]
numeradores = scanl (*) 1 [2..]

-- denominadores es la sucesión de los denominadores. Por ejemplo,
--    λ> take 10 denominadores
--    [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]
denominadores :: [Integer]
denominadores = scanl (*) 3 [5, 7..]
 
-- 4ª solución
-- ===========

aproximaPi4 :: Integer -> Double
aproximaPi4 n = 
  read (x : "." ++ xs)
  where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer
aproximaPi4' n = 
  2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores))) 
  where p = 10^n

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_aproximaPi :: Property
prop_aproximaPi =
  forAll (arbitrary `suchThat` (> 3)) $ \n ->
  all (=~ aproximaPi1 n)
      [aproximaPi2 n,
       aproximaPi3 n,
       aproximaPi4 n]

(=~) :: Double -> Double -> Bool
x =~ y = abs (x - y) < 0.001

-- La comprobación es
--    λ> quickCheck prop_aproximaPi
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> aproximaPi1 3000 
--    3.141592653589793
--    (4.96 secs, 27,681,824,408 bytes)
--    λ> aproximaPi2 3000 
--    3.1415926535897922
--    (3.00 secs, 20,496,194,496 bytes)
--    λ> aproximaPi3 3000 
--    3.141592653589793
--    (3.13 secs, 13,439,528,432 bytes)
--    λ> aproximaPi4 3000 
--    3.141592653589793
--    (0.09 secs, 23,142,144 bytes)

import Data.Ratio ((%))

import Data.List (genericTake)

import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

-- 1ª solución

-- ===========

aproximaPi1 :: Integer -> Double

aproximaPi1 n =

fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

-- 2ª solución

-- ===========

aproximaPi2 :: Integer -> Double

aproximaPi2 0 = 2

aproximaPi2 n =

aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

-- 3ª solución

-- ===========

aproximaPi3 :: Integer -> Double

aproximaPi3 n =

fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,

-- λ> take 10 numeradores

-- [1,2,6,24,120,720,5040,40320,362880,3628800]

numeradores :: [Integer]

numeradores = scanl (*) 1 [2..]

-- denominadores es la sucesión de los denominadores. Por ejemplo,

-- λ> take 10 denominadores

-- [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]

denominadores :: [Integer]

denominadores = scanl (*) 3 [5, 7..]

-- 4ª solución

-- ===========

aproximaPi4 :: Integer -> Double

aproximaPi4 n =

read (x : "." ++ xs)

where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer

aproximaPi4' n =

2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores)))

where p = 10^n

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_aproximaPi :: Property

prop_aproximaPi =

forAll (arbitrary `suchThat` (> 3)) $ \n ->

all (=~ aproximaPi1 n)

[aproximaPi2 n,

aproximaPi3 n,

aproximaPi4 n]

(=~) :: Double -> Double -> Bool

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheck prop_aproximaPi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> aproximaPi1 3000

-- 3.141592653589793

-- (4.96 secs, 27,681,824,408 bytes)

-- λ> aproximaPi2 3000

-- 3.1415926535897922

-- (3.00 secs, 20,496,194,496 bytes)

-- λ> aproximaPi3 3000

-- 3.141592653589793

-- (3.13 secs, 13,439,528,432 bytes)

-- λ> aproximaPi4 3000

-- 3.141592653589793

-- (0.09 secs, 23,142,144 bytes)

El código se encuentra en GitHub.

Pandigitales primos

PorJosé A. Alonso 29-junio-202221-junio-2022

Un número con n dígitos es pandigital si contiene todos los dígitos del 1 a n exactamente una vez. Por ejemplo, 2143 es un pandigital con 4 dígitos y, además, es primo.

Definir la lista

   pandigitalesPrimos :: [Int]

1	pandigitalesPrimos :: [Int]

tal que sus elementos son los números pandigitales primos, ordenados de mayor a menor. Por ejemplo,

   take 3 pandigitalesPrimos       ==  [7652413,7642513,7641253]
   2143 `elem` pandigitalesPrimos  ==  True
   length pandigitalesPrimos       ==  538

take 3 pandigitalesPrimos == [7652413,7642513,7641253]

2143 `elem` pandigitalesPrimos == True

length pandigitalesPrimos == 538

Soluciones

import Data.List (permutations, sort)
import Data.Char (intToDigit)
import Data.Numbers.Primes (isPrime)

-- 1ª solución
-- ===========

pandigitalesPrimos1 :: [Int]
pandigitalesPrimos1 =
  concatMap nPandigitalesPrimos1 [9,8..1]

-- (nPandigitalesPrimos n) es la lista de los números pandigitales con n
-- dígitos, ordenada de mayor a menor. Por ejemplo,
--    nPandigitalesPrimos 4  ==  [4231,2341,2143,1423]
--    nPandigitalesPrimos 5  ==  []
nPandigitalesPrimos1 :: Int -> [Int]
nPandigitalesPrimos1 n = filter isPrime (pandigitales n)

-- (pandigitales n) es la lista de los números pandigitales de n dígitos
-- ordenada de mayor a menor. Por ejemplo,
--    pandigitales 3  ==  [321,312,231,213,132,123]
pandigitales :: Int -> [Int]
pandigitales n = 
    reverse $ sort $ map digitosAentero (permutations [1..n])

-- (digitosAentero ns) es el número cuyos dígitos son ns. Por ejemplo,
--    digitosAentero [3,2,5]  ==  325
digitosAentero :: [Int] -> Int
digitosAentero = read . map intToDigit

-- 2ª solución
-- ===========

pandigitalesPrimos2 :: [Int]
pandigitalesPrimos2 =
  concatMap nPandigitalesPrimos2 [9,8..1]

-- Nota. La definición de nPandigitalesPrimos1 se puede simplificar, ya
-- que la suma de los números de 1 a n es divisible por 3, entonces los
-- números  pandigitales con n dígitos también lo son y, por tanto, no
-- son primos. 
nPandigitalesPrimos2 :: Int -> [Int]
nPandigitalesPrimos2 n 
  | sum [1..n] `mod` 3 == 0 = []
  | otherwise               = filter isPrime (pandigitales n)

-- 2ª solución
-- ===========

pandigitalesPrimos3 :: [Int]
pandigitalesPrimos3 =
  concatMap nPandigitalesPrimos3 [9,8..1]

-- La definición de nPandigitales se puede simplificar, ya que
--    λ> [n | n <- [1..9], sum [1..n] `mod` 3 /= 0]
--    [1,4,7]
nPandigitalesPrimos3 :: Int -> [Int]
nPandigitalesPrimos3 n 
  | n `elem` [4,7] = filter isPrime (pandigitales n)
  | otherwise      = []

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_pandigitalesPrimos :: Bool
prop_pandigitalesPrimos =
  all (== pandigitalesPrimos1)
      [pandigitalesPrimos2,
       pandigitalesPrimos3]

-- La comprobación es
--    λ> prop_pandigitalesPrimos
--    True

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (pandigitalesPrimos1)
--    538
--    (1.44 secs, 5,249,850,032 bytes)
--    λ> length (pandigitalesPrimos2)
--    538
--    (0.14 secs, 619,249,632 bytes)
--    λ> length (pandigitalesPrimos3)
--    538
--    (0.14 secs, 619,237,464 bytes)

import Data.List (permutations, sort)

import Data.Char (intToDigit)

import Data.Numbers.Primes (isPrime)

-- 1ª solución

-- ===========

pandigitalesPrimos1 :: [Int]

pandigitalesPrimos1 =

concatMap nPandigitalesPrimos1 [9,8..1]

-- (nPandigitalesPrimos n) es la lista de los números pandigitales con n

-- dígitos, ordenada de mayor a menor. Por ejemplo,

-- nPandigitalesPrimos 4 == [4231,2341,2143,1423]

-- nPandigitalesPrimos 5 == []

nPandigitalesPrimos1 :: Int -> [Int]

nPandigitalesPrimos1 n = filter isPrime (pandigitales n)

-- (pandigitales n) es la lista de los números pandigitales de n dígitos

-- ordenada de mayor a menor. Por ejemplo,

-- pandigitales 3 == [321,312,231,213,132,123]

pandigitales :: Int -> [Int]

pandigitales n =

reverse $ sort $ map digitosAentero (permutations [1..n])

-- (digitosAentero ns) es el número cuyos dígitos son ns. Por ejemplo,

-- digitosAentero [3,2,5] == 325

digitosAentero :: [Int] -> Int

digitosAentero = read . map intToDigit

-- 2ª solución

-- ===========

pandigitalesPrimos2 :: [Int]

pandigitalesPrimos2 =

concatMap nPandigitalesPrimos2 [9,8..1]

-- Nota. La definición de nPandigitalesPrimos1 se puede simplificar, ya

-- que la suma de los números de 1 a n es divisible por 3, entonces los

-- números pandigitales con n dígitos también lo son y, por tanto, no

-- son primos.

nPandigitalesPrimos2 :: Int -> [Int]

nPandigitalesPrimos2 n

| sum [1..n] `mod` 3 == 0 = []

| otherwise = filter isPrime (pandigitales n)

-- 2ª solución

-- ===========

pandigitalesPrimos3 :: [Int]

pandigitalesPrimos3 =

concatMap nPandigitalesPrimos3 [9,8..1]

-- La definición de nPandigitales se puede simplificar, ya que

-- λ> [n | n <- [1..9], sum [1..n] `mod` 3 /= 0]

-- [1,4,7]

nPandigitalesPrimos3 :: Int -> [Int]

nPandigitalesPrimos3 n

| n `elem` [4,7] = filter isPrime (pandigitales n)

| otherwise = []

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_pandigitalesPrimos :: Bool

prop_pandigitalesPrimos =

all (== pandigitalesPrimos1)

[pandigitalesPrimos2,

pandigitalesPrimos3]

-- La comprobación es

-- λ> prop_pandigitalesPrimos

-- True

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (pandigitalesPrimos1)

-- 538

-- (1.44 secs, 5,249,850,032 bytes)

-- λ> length (pandigitalesPrimos2)

-- 538

-- (0.14 secs, 619,249,632 bytes)

-- λ> length (pandigitalesPrimos3)

-- 538

-- (0.14 secs, 619,237,464 bytes)

El código se encuentra en GitHub.

Codificación de Fibonacci

PorJosé A. Alonso 28-junio-202221-junio-2022

La codificación de Fibonacci http://bit.ly/1Lllqjv de un número n es una cadena d = d(0)d(1)…d(k-1)d(k) de ceros y unos tal que

   n = d(0)·F(2) + d(1)·F(3) +...+ d(k-1)·F(k+1) 
   d(k-1) = d(k) = 1

1 2	n = d(0)·F(2) + d(1)·F(3) +...+ d(k-1)·F(k+1) d(k-1) = d(k) = 1

donde F(i) es el i-ésimo término de la sucesión de Fibonacci

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...

Por ejemplo, la codificación de Fibonacci de 4 es «1011» ya que los dos últimos elementos son iguales a 1 y

   1·F(2) + 0·F(3) + 1·F(4) = 1·1 + 0·2 + 1·3 = 4

1	1·F(2) + 0·F(3) + 1·F(4) = 1·1 + 0·2 + 1·3 = 4

La codificación de Fibonacci de los primeros números se muestra en la siguiente tabla

    1  = 1     = F(2)           ≡       11
    2  = 2     = F(3)           ≡      011
    3  = 3     = F(4)           ≡     0011
    4  = 1+3   = F(2)+F(4)      ≡     1011
    5  = 5     = F(5)           ≡    00011
    6  = 1+5   = F(2)+F(5)      ≡    10011
    7  = 2+5   = F(3)+F(5)      ≡    01011
    8  = 8     = F(6)           ≡   000011
    9  = 1+8   = F(2)+F(6)      ≡   100011
   10  = 2+8   = F(3)+F(6)      ≡   010011
   11  = 3+8   = F(4)+F(6)      ≡   001011
   12  = 1+3+8 = F(2)+F(4)+F(6) ≡   101011
   13  = 13    = F(7)           ≡  0000011
   14  = 1+13  = F(2)+F(7)      ≡  1000011

1 = 1 = F(2) ≡ 11

2 = 2 = F(3) ≡ 011

3 = 3 = F(4) ≡ 0011

4 = 1+3 = F(2)+F(4) ≡ 1011

5 = 5 = F(5) ≡ 00011

6 = 1+5 = F(2)+F(5) ≡ 10011

7 = 2+5 = F(3)+F(5) ≡ 01011

8 = 8 = F(6) ≡ 000011

9 = 1+8 = F(2)+F(6) ≡ 100011

10 = 2+8 = F(3)+F(6) ≡ 010011

11 = 3+8 = F(4)+F(6) ≡ 001011

12 = 1+3+8 = F(2)+F(4)+F(6) ≡ 101011

13 = 13 = F(7) ≡ 0000011

14 = 1+13 = F(2)+F(7) ≡ 1000011

Definir la función

   codigoFib :: Integer -> String

1	codigoFib :: Integer -> String

tal que (codigoFib n) es la codificación de Fibonacci del número n. Por ejemplo,

   λ> codigoFib 65
   "0100100011"
   λ> [codigoFib n | n <- [1..7]]
   ["11","011","0011","1011","00011","10011","01011"]

λ> codigoFib 65

"0100100011"

λ> [codigoFib n | n <- [1..7]]

["11","011","0011","1011","00011","10011","01011"]

Comprobar con QuickCheck las siguientes propiedades:

Todo entero positivo se puede descomponer en suma de números de la sucesión de Fibonacci.
Las codificaciones de Fibonacci tienen como mínimo 2 elementos.
En las codificaciones de Fibonacci, la cadena «11» sólo aparece una vez y la única vez que aparece es al final.

Soluciones

import Data.List (isInfixOf)
import Data.Array (Array, accumArray, elems)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

codigoFib1 :: Integer -> String
codigoFib1 = concatMap show . codificaFibLista

-- (codificaFibLista n) es la lista correspondiente a la codificación de
-- Fibonacci del número n. Por ejemplo,
--    λ> codificaFibLista 65
--    [0,1,0,0,1,0,0,0,1,1]
--    λ> [codificaFibLista n | n <- [1..7]]
--    [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]
codificaFibLista :: Integer -> [Integer]
codificaFibLista n = map f [2..head xs] ++ [1]
  where xs = map fst (descomposicion n)
        f i | i `elem` xs = 1
            | otherwise = 0

-- (descomposicion n) es la lista de pares (i,f) tales que f es el
-- i-ésimo número de Fibonacci y las segundas componentes es una
-- sucesión decreciente de números de Fibonacci cuya suma es n. Por
-- ejemplo, 
--    descomposicion 65  ==  [(10,55),(6,8),(3,2)]
--    descomposicion 66  ==  [(10,55),(6,8),(4,3)]
descomposicion :: Integer -> [(Integer, Integer)]
descomposicion 0 = []
descomposicion 1 = [(2,1)]
descomposicion n = (i,x) : descomposicion (n-x)
  where (i,x) = fibAnterior n

-- (fibAnterior n) es el mayor número de Fibonacci menor o igual que
-- n. Por ejemplo,
--    fibAnterior 33  ==  (8,21)
--    fibAnterior 34  ==  (9,34)
fibAnterior :: Integer -> (Integer, Integer)
fibAnterior n = last (takeWhile p fibsConIndice)
  where p (_,x) = x <= n

-- fibsConIndice es la sucesión de los números de Fibonacci junto con
-- sus índices. Por ejemplo,
--    λ> take 10 fibsConIndice
--    [(0,0),(1,1),(2,1),(3,2),(4,3),(5,5),(6,8),(7,13),(8,21),(9,34)]
fibsConIndice :: [(Integer, Integer)]
fibsConIndice = zip [0..] fibs

-- fibs es la sucesión de Fibonacci. Por ejemplo, 
--    take 10 fibs  ==  [0,1,1,2,3,5,8,13,21,34]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

--- 2ª solución
-- ============

codigoFib2 :: Integer -> String
codigoFib2 = concatMap show . elems . codificaFibVec

-- (codificaFibVec n) es el vector correspondiente a la codificación de
-- Fibonacci del número n. Por ejemplo,
--    λ> codificaFibVec 65
--    array (0,9) [(0,0),(1,1),(2,0),(3,0),(4,1),(5,0),(6,0),(7,0),(8,1),(9,1)]
--    λ> [elems (codificaFibVec n) | n <- [1..7]]
--    [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]
codificaFibVec :: Integer -> Array Integer Integer
codificaFibVec n = accumArray (+) 0 (0,a+1) ((a+1,1):is) 
  where is = [(i-2,1) | (i,_) <- descomposicion n]
        a  = fst (head is)

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_codigoFib :: Positive Integer -> Bool 
prop_codigoFib (Positive n) =
  codigoFib1 n == codigoFib2 n
  
-- La comprobación es
--    λ> quickCheck prop_codigoFib
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> head [n | n <- [1..], length (codigoFib1 n) > 25]
--    121393
--    (4.30 secs, 3,031,108,104 bytes)
--    λ> head [n | n <- [1..], length (codigoFib2 n) > 25]
--    121393
--    (3.46 secs, 2,505,869,616 bytes)

-- Propiedades
-- ===========

-- Usaremos la 2ª definición
codigoFib :: Integer -> String
codigoFib = codigoFib2

-- Prop.: La función descomposicion es correcta:
prop_descomposicion_correcta :: Positive Integer -> Bool
prop_descomposicion_correcta (Positive n) =
  n == sum (map snd (descomposicion n))

-- La comprobación es
--    λ> quickCheck prop_descomposicion_correcta
--    +++ OK, passed 100 tests.

-- Prop.: Todo entero positivo se puede descomponer en suma de números de
-- la sucesión de Fibonacci.
prop_descomposicion :: Positive Integer -> Bool
prop_descomposicion (Positive n) =
  not (null (descomposicion n))

-- La comprobación es
--    λ> quickCheck prop_descomposicion
--    +++ OK, passed 100 tests.

-- Prop.: Las codificaciones de Fibonacci tienen como mínimo 2 elementos.
prop_length_codigoFib :: Positive Integer -> Bool
prop_length_codigoFib (Positive n) =
  length (codigoFib n) >= 2

-- La comprobación es
--    λ> quickCheck prop_length_codigoFib
--    +++ OK, passed 100 tests.

-- Prop.: En las codificaciones de Fibonacci, la cadena "11" sólo
-- aparece una vez y la única vez que aparece es al final.
prop3_cadena_11_en_codigoFib :: Positive Integer -> Bool
prop3_cadena_11_en_codigoFib (Positive n) = 
  take 2 xs == "11" && not ("11" `isInfixOf` drop 2 xs)
  where xs = reverse (codigoFib n)

-- La comprobación es
--    λ> quickCheck prop3_cadena_11_en_codigoFib
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

import Data.List (isInfixOf)

import Data.Array (Array, accumArray, elems)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

codigoFib1 :: Integer -> String

codigoFib1 = concatMap show . codificaFibLista

-- (codificaFibLista n) es la lista correspondiente a la codificación de

-- Fibonacci del número n. Por ejemplo,

-- λ> codificaFibLista 65

-- [0,1,0,0,1,0,0,0,1,1]

-- λ> [codificaFibLista n | n <- [1..7]]

-- [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]

codificaFibLista :: Integer -> [Integer]

codificaFibLista n = map f [2..head xs] ++ [1]

where xs = map fst (descomposicion n)

f i | i `elem` xs = 1

| otherwise = 0

-- (descomposicion n) es la lista de pares (i,f) tales que f es el

-- i-ésimo número de Fibonacci y las segundas componentes es una

-- sucesión decreciente de números de Fibonacci cuya suma es n. Por

-- ejemplo,

-- descomposicion 65 == [(10,55),(6,8),(3,2)]

-- descomposicion 66 == [(10,55),(6,8),(4,3)]

descomposicion :: Integer -> [(Integer, Integer)]

descomposicion 0 = []

descomposicion 1 = [(2,1)]

descomposicion n = (i,x) : descomposicion (n-x)

where (i,x) = fibAnterior n

-- (fibAnterior n) es el mayor número de Fibonacci menor o igual que

-- n. Por ejemplo,

-- fibAnterior 33 == (8,21)

-- fibAnterior 34 == (9,34)

fibAnterior :: Integer -> (Integer, Integer)

fibAnterior n = last (takeWhile p fibsConIndice)

where p (_,x) = x <= n

-- fibsConIndice es la sucesión de los números de Fibonacci junto con

-- sus índices. Por ejemplo,

-- λ> take 10 fibsConIndice

-- [(0,0),(1,1),(2,1),(3,2),(4,3),(5,5),(6,8),(7,13),(8,21),(9,34)]

fibsConIndice :: [(Integer, Integer)]

fibsConIndice = zip [0..] fibs

-- fibs es la sucesión de Fibonacci. Por ejemplo,

-- take 10 fibs == [0,1,1,2,3,5,8,13,21,34]

fibs :: [Integer]

fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

--- 2ª solución

-- ============

codigoFib2 :: Integer -> String

codigoFib2 = concatMap show . elems . codificaFibVec

-- (codificaFibVec n) es el vector correspondiente a la codificación de

-- Fibonacci del número n. Por ejemplo,

-- λ> codificaFibVec 65

-- array (0,9) [(0,0),(1,1),(2,0),(3,0),(4,1),(5,0),(6,0),(7,0),(8,1),(9,1)]

-- λ> [elems (codificaFibVec n) | n <- [1..7]]

-- [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]

codificaFibVec :: Integer -> Array Integer Integer

codificaFibVec n = accumArray (+) 0 (0,a+1) ((a+1,1):is)

where is = [(i-2,1) | (i,_) <- descomposicion n]

a = fst (head is)

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_codigoFib :: Positive Integer -> Bool

prop_codigoFib (Positive n) =

codigoFib1 n == codigoFib2 n

-- La comprobación es

-- λ> quickCheck prop_codigoFib

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> head [n | n <- [1..], length (codigoFib1 n) > 25]

-- 121393

-- (4.30 secs, 3,031,108,104 bytes)

-- λ> head [n | n <- [1..], length (codigoFib2 n) > 25]

-- 121393

-- (3.46 secs, 2,505,869,616 bytes)

-- Propiedades

-- ===========

-- Usaremos la 2ª definición

codigoFib :: Integer -> String

codigoFib = codigoFib2

-- Prop.: La función descomposicion es correcta:

prop_descomposicion_correcta :: Positive Integer -> Bool

prop_descomposicion_correcta (Positive n) =

n == sum (map snd (descomposicion n))

-- La comprobación es

-- λ> quickCheck prop_descomposicion_correcta

-- +++ OK, passed 100 tests.

-- Prop.: Todo entero positivo se puede descomponer en suma de números de

-- la sucesión de Fibonacci.

prop_descomposicion :: Positive Integer -> Bool

prop_descomposicion (Positive n) =

not (null (descomposicion n))

-- La comprobación es

-- λ> quickCheck prop_descomposicion

-- +++ OK, passed 100 tests.

-- Prop.: Las codificaciones de Fibonacci tienen como mínimo 2 elementos.

prop_length_codigoFib :: Positive Integer -> Bool

prop_length_codigoFib (Positive n) =

length (codigoFib n) >= 2

-- La comprobación es

-- λ> quickCheck prop_length_codigoFib

-- +++ OK, passed 100 tests.

-- Prop.: En las codificaciones de Fibonacci, la cadena "11" sólo

-- aparece una vez y la única vez que aparece es al final.

prop3_cadena_11_en_codigoFib :: Positive Integer -> Bool

prop3_cadena_11_en_codigoFib (Positive n) =

take 2 xs == "11" && not ("11" `isInfixOf` drop 2 xs)

where xs = reverse (codigoFib n)

-- La comprobación es

-- λ> quickCheck prop3_cadena_11_en_codigoFib

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

La sucesión del reloj astronómico de Praga

PorJosé A. Alonso 27-junio-202221-junio-2022

La cadena infinita «1234321234321234321…», formada por la repetición de los dígitos 123432, tiene una propiedad (en la que se basa el funcionamiento del reloj astronómico de Praga: la cadena se puede partir en una sucesión de números, de forma que la suma de los dígitos de dichos números es la sucesión de los números naturales, como se observa a continuación:

    1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, ...
    1, 2, 3, 4,  5,   6,  7,    8,   9,   10,    11, ...

1 2	1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...

Definir la lista

   reloj :: [Integer]

1	reloj :: [Integer]

cuyos elementos son los términos de la sucesión anterior. Por ejemplo,

   λ> take 11 reloj
   [1,2,3,4,32,123,43,2123,432,1234,32123]
   λ> (reloj !! 1000) `mod` (10^50)
   23432123432123432123432123432123432123432123432123

λ> take 11 reloj

[1,2,3,4,32,123,43,2123,432,1234,32123]

λ> (reloj !! 1000) `mod` (10^50)

23432123432123432123432123432123432123432123432123

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución
-- ===========

reloj1 :: [Integer]
reloj1 = aux [1..] (cycle "123432")
  where aux (n:ns) xs = read ys : aux ns zs
          where (ys,zs) = prefijoSuma n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo
-- suma es n y el resto de xs. Por ejemplo,
--    prefijoSuma 6 "12343"  ==  ("123","43")
prefijoSuma :: Int -> String -> (String,String)
prefijoSuma n xs = 
  head [(us,vs) | (us,vs) <- zip (inits xs) (tails xs)
                , sumaD us == n]

-- (sumaD xs) es la suma de los dígitos de xs. Por ejemplo,
--    sumaD "123"  ==  6
sumaD :: String -> Int
sumaD = sum . map digitToInt

-- 2ª solución
-- ===========

reloj2 :: [Integer]
reloj2 = aux [1..] (cycle "123432")
  where aux (n:ns) xs = read ys : aux ns zs
          where (ys,zs) = prefijoSuma2 n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo
-- suma es n y el resto de xs. Por ejemplo,
--    prefijoSuma2 6 "12343"  ==  ("123","43")
prefijoSuma2 :: Int -> String -> (String,String)
prefijoSuma2 n (x:xs) 
  | y == n = ([x],xs)
  | otherwise = (x:ys,zs) 
  where y       = read [x]
        (ys,zs) = prefijoSuma2 (n-y) xs

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_reloj :: NonNegative Int -> Bool
prop_reloj (NonNegative n) =
  reloj1 !! n == reloj2 !! n
  
-- La comprobación es
--    λ> quickCheck prop_reloj
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> (reloj1 !! 1000) `mod` (10^9)
--    123432123
--    (2.47 secs, 5,797,620,784 bytes)
--    λ> (reloj2 !! 1000) `mod` (10^9)
--    123432123
--    (0.44 secs, 798,841,528 bytes)

import Data.List (inits, tails)

import Data.Char (digitToInt)

import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución

-- ===========

reloj1 :: [Integer]

reloj1 = aux [1..] (cycle "123432")

where aux (n:ns) xs = read ys : aux ns zs

where (ys,zs) = prefijoSuma n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo

-- suma es n y el resto de xs. Por ejemplo,

-- prefijoSuma 6 "12343" == ("123","43")

prefijoSuma :: Int -> String -> (String,String)

prefijoSuma n xs =

head [(us,vs) | (us,vs) <- zip (inits xs) (tails xs)

, sumaD us == n]

-- (sumaD xs) es la suma de los dígitos de xs. Por ejemplo,

-- sumaD "123" == 6

sumaD :: String -> Int

sumaD = sum . map digitToInt

-- 2ª solución

-- ===========

reloj2 :: [Integer]

reloj2 = aux [1..] (cycle "123432")

where aux (n:ns) xs = read ys : aux ns zs

where (ys,zs) = prefijoSuma2 n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo

-- suma es n y el resto de xs. Por ejemplo,

-- prefijoSuma2 6 "12343" == ("123","43")

prefijoSuma2 :: Int -> String -> (String,String)

prefijoSuma2 n (x:xs)

| y == n = ([x],xs)

| otherwise = (x:ys,zs)

where y = read [x]

(ys,zs) = prefijoSuma2 (n-y) xs

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_reloj :: NonNegative Int -> Bool

prop_reloj (NonNegative n) =

reloj1 !! n == reloj2 !! n

-- La comprobación es

-- λ> quickCheck prop_reloj

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> (reloj1 !! 1000) `mod` (10^9)

-- 123432123

-- (2.47 secs, 5,797,620,784 bytes)

-- λ> (reloj2 !! 1000) `mod` (10^9)

-- 123432123

-- (0.44 secs, 798,841,528 bytes)

El código se encuentra en GitHub.

Método de bisección para aproximar raíces de funciones

PorJosé A. Alonso 24-junio-202221-junio-2022

El método de bisección para calcular un cero de una función en el intervalo [a,b] se basa en el teorema de Bolzano:

«Si f(x) es una función continua en el intervalo [a, b], y si, además, en los extremos del intervalo la función f(x) toma valores de signo opuesto (f(a) * f(b) < 0), entonces existe al menos un valor c en (a, b) para el que f(c) = 0".

El método para calcular un cero de la función f en el intervalo [a,b] con un error menor que e consiste en tomar el punto medio del intervalo c = (a+b)/2 y considerar los siguientes casos:

Si |f(c)| < e, hemos encontrado una aproximación del punto que anula f en el intervalo con un error aceptable.
Si f(c) tiene signo distinto de f(a), repetir el proceso en el intervalo [a,c].
Si no, repetir el proceso en el intervalo [c,b].

Definir la función

   biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

1	biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

tal que (biseccion f a b e) es una aproximación del punto del intervalo [a,b] en el que se anula la función f, con un error menor que e, calculada mediante el método de la bisección. Por ejemplo,

   biseccion (\x -> x^2 - 3) 0 5 0.01             ==  1.7333984375
   biseccion (\x -> x^3 - x - 2) 0 4 0.01         ==  1.521484375
   biseccion cos 0 2 0.01                         ==  1.5625
   biseccion (\x -> log (50-x) - 4) (-10) 3 0.01  ==  -5.125

biseccion (\x -> x^2 - 3) 0 5 0.01 == 1.7333984375

biseccion (\x -> x^3 - x - 2) 0 4 0.01 == 1.521484375

biseccion cos 0 2 0.01 == 1.5625

biseccion (\x -> log (50-x) - 4) (-10) 3 0.01 == -5.125

Soluciones

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

-- 1ª solución
-- ===========

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion1 f a b e  
    | abs (f c) < e = c
    | f a * f c < 0 = biseccion1 f a c e
    | otherwise     = biseccion1 f c b e
    where c = (a+b)/2

-- 2ª solución
-- ===========

biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion2 f a b e = aux a b
  where aux a' b' | abs (f c) < e = c
                  | f a' * f c < 0 = aux a' c 
                  | otherwise     = aux c b'
          where c = (a'+b')/2

-- Comprobación de equivalencia
-- ============================

newtype Polinomio = P [Int]
  deriving Show

valorPolinomio :: Polinomio -> Double -> Double
valorPolinomio (P cs) x =
  sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio
polinomioArbitrario 0 = return (P [])
polinomioArbitrario n = do
  c <- choose (-10,10)
  (P xs) <- polinomioArbitrario (n `div` 2)
  return (P (c:xs))

instance Arbitrary Polinomio where
  arbitrary = sized polinomioArbitrario

-- La propiedad es
prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property
prop_biseccion p a b e =
  f a * f b < 0 && e > 0 ==>
  biseccion1 f a b e =~ biseccion2 f a b e
  where
    f = valorPolinomio p 
    x =~ y = abs (x - y) < 0.001

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion
--    +++ OK, passed 100 tests; 2156 discarded.

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

-- 1ª solución

-- ===========

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion1 f a b e

| abs (f c) < e = c

| f a * f c < 0 = biseccion1 f a c e

| otherwise = biseccion1 f c b e

where c = (a+b)/2

-- 2ª solución

-- ===========

biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion2 f a b e = aux a b

where aux a' b' | abs (f c) < e = c

| f a' * f c < 0 = aux a' c

| otherwise = aux c b'

where c = (a'+b')/2

-- Comprobación de equivalencia

-- ============================

newtype Polinomio = P [Int]

deriving Show

valorPolinomio :: Polinomio -> Double -> Double

valorPolinomio (P cs) x =

sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio

polinomioArbitrario 0 = return (P [])

polinomioArbitrario n = do

c <- choose (-10,10)

(P xs) <- polinomioArbitrario (n `div` 2)

return (P (c:xs))

instance Arbitrary Polinomio where

arbitrary = sized polinomioArbitrario

-- La propiedad es

prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property

prop_biseccion p a b e =

f a * f b < 0 && e > 0 ==>

biseccion1 f a b e =~ biseccion2 f a b e

where

f = valorPolinomio p

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion

-- +++ OK, passed 100 tests; 2156 discarded.

El código se encuentra en GitHub.

Números para los que mcm(1,2,…n-1) = mcm(1,2,…,n)

PorJosé A. Alonso 23-junio-202221-junio-2022

Un número n es especial si mcm(1,2,…,n-1) = mcm(1,2,…,n). Por ejemplo, el 6 es especial ya que

   mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6)

1	mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6)

Definir la sucesión

   especiales :: [Integer]

1	especiales :: [Integer]

cuyos términos son los números especiales. Por ejemplo,

   take 10 especiales     ==  [1,6,10,12,14,15,18,20,21,22]
   especiales !! 50       ==  84
   especiales !! 500      ==  638
   especiales !! 5000     ==  5806
   especiales !! 50000    ==  55746

take 10 especiales == [1,6,10,12,14,15,18,20,21,22]

especiales !! 50 == 84

especiales !! 500 == 638

especiales !! 5000 == 5806

especiales !! 50000 == 55746

Soluciones

import Test.QuickCheck (NonNegative(NonNegative), quickCheck)

-- 1ª solución
-- ===========

especiales1 :: [Integer]
especiales1 = filter especial1 [1..]

especial1 :: Integer -> Bool
especial1 n = mcm1 [1..n-1] == mcm1 [1..n]

mcm1 :: [Integer] -> Integer
mcm1 []     = 1
mcm1 (x:xs) = lcm x (mcm1 xs)

-- 2ª solución
-- ===========

especiales2 :: [Integer]
especiales2 = filter especial2 [1..]

especial2 :: Integer -> Bool
especial2 n = mcm2 [1..n-1] == mcm2 [1..n]

mcm2 :: [Integer] -> Integer
mcm2 = foldr lcm 1

-- 3ª solución
-- ===========

especiales3 :: [Integer]
especiales3 = [n | ((n,x),(_,y)) <- zip mcms (tail mcms)
                 , x == y]

mcms :: [(Integer,Integer)]
mcms = zip [1..] (scanl lcm 1 [1..])

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_especiales :: NonNegative Int -> Bool
prop_especiales (NonNegative n) =
  all (== especiales1 !! n)
      [especiales2 !! n,
       especiales3 !! n]

-- La comprobación es
--    λ> quickCheck prop_especiales
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es

-- Comparación 
--    λ> especiales1 !! 2000
--    2390
--    (3.38 secs, 4,724,497,192 bytes)
--    λ> especiales2 !! 2000
--    2390
--    (1.91 secs, 4,303,415,512 bytes)
--    λ> especiales3 !! 2000
--    2390
--    (0.01 secs, 4,209,664 bytes)

import Test.QuickCheck (NonNegative(NonNegative), quickCheck)

-- 1ª solución

-- ===========

especiales1 :: [Integer]

especiales1 = filter especial1 [1..]

especial1 :: Integer -> Bool

especial1 n = mcm1 [1..n-1] == mcm1 [1..n]

mcm1 :: [Integer] -> Integer

mcm1 [] = 1

mcm1 (x:xs) = lcm x (mcm1 xs)

-- 2ª solución

-- ===========

especiales2 :: [Integer]

especiales2 = filter especial2 [1..]

especial2 :: Integer -> Bool

especial2 n = mcm2 [1..n-1] == mcm2 [1..n]

mcm2 :: [Integer] -> Integer

mcm2 = foldr lcm 1

-- 3ª solución

-- ===========

especiales3 :: [Integer]

especiales3 = [n | ((n,x),(_,y)) <- zip mcms (tail mcms)

, x == y]

mcms :: [(Integer,Integer)]

mcms = zip [1..] (scanl lcm 1 [1..])

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_especiales :: NonNegative Int -> Bool

prop_especiales (NonNegative n) =

all (== especiales1 !! n)

[especiales2 !! n,

especiales3 !! n]

-- La comprobación es

-- λ> quickCheck prop_especiales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- Comparación

-- λ> especiales1 !! 2000

-- 2390

-- (3.38 secs, 4,724,497,192 bytes)

-- λ> especiales2 !! 2000

-- 2390

-- (1.91 secs, 4,303,415,512 bytes)

-- λ> especiales3 !! 2000

-- 2390

-- (0.01 secs, 4,209,664 bytes)

El código se encuentra en GitHub.

Cálculo de la suma 11! + 22! + 33! + … + nn!

PorJosé A. Alonso 22-junio-202221-junio-2022

Definir la función

   suma :: Integer -> Integer

1	suma :: Integer -> Integer

tal que (suma n) es la suma 1·1! + 2·2! + 3·3! + ... + n·n!. Por ejemplo,

   suma 1  ==  1
   suma 2  ==  5
   suma 3  ==  23
   suma 4  ==  119
   suma 5  ==  719
   take 9 (show (suma 70000))  ==  "823780458"

suma 1 == 1

suma 2 == 5

suma 3 == 23

suma 4 == 119

suma 5 == 719

take 9 (show (suma 70000)) == "823780458"

Soluciones

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

suma1 :: Integer -> Integer
suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer
factorial n = product [1..n]

-- 2ª solución
-- ===========

suma2 :: Integer -> Integer
suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]
factoriales = scanl (*) 1 [2..]

-- 3ª solución
-- ===========

-- Basada en los siguientes cálculos
--    λ> [suma1 n | n <- [0..10]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]
--    λ> [factorial n | n <- [0..10]]
--    [1,1,2,6,24,120,720,5040,40320,362880,3628800]
--    λ> [factorial n | n <- [1..11]]
--    [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]
--    λ> [factorial n - 1 | n <- [1..11]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer
suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_suma :: Positive Integer -> Bool
prop_suma (Positive n) =
  all (== suma1 n)
      [suma2 n,
       suma3 n]
  
-- La comprobación es
--    λ> quickCheck prop_suma
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> take 5 (show (suma1 4000))
--    "73170"
--    (5.04 secs, 16,225,195,448 bytes)
--    λ> take 5 (show (suma2 4000))
--    "73170"
--    (0.08 secs, 35,862,152 bytes)
--    λ> take 5 (show (suma3 4000))
--    "73170"
--    (0.01 secs, 12,896,968 bytes)
--    
--    
--    λ> take 5 (show (suma2 40000))
--    "83669"
--    (1.82 secs, 4,549,612,264 bytes)
--    λ> take 5 (show (suma3 40000))
--    "83669"
--    (0.24 secs, 1,620,976,984 bytes)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

suma1 :: Integer -> Integer

suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª solución

-- ===========

suma2 :: Integer -> Integer

suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]

factoriales = scanl (*) 1 [2..]

-- 3ª solución

-- ===========

-- Basada en los siguientes cálculos

-- λ> [suma1 n | n <- [0..10]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

-- λ> [factorial n | n <- [0..10]]

-- [1,1,2,6,24,120,720,5040,40320,362880,3628800]

-- λ> [factorial n | n <- [1..11]]

-- [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]

-- λ> [factorial n - 1 | n <- [1..11]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer

suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_suma :: Positive Integer -> Bool

prop_suma (Positive n) =

all (== suma1 n)

[suma2 n,

suma3 n]

-- La comprobación es

-- λ> quickCheck prop_suma

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> take 5 (show (suma1 4000))

-- "73170"

-- (5.04 secs, 16,225,195,448 bytes)

-- λ> take 5 (show (suma2 4000))

-- "73170"

-- (0.08 secs, 35,862,152 bytes)

-- λ> take 5 (show (suma3 4000))

-- "73170"

-- (0.01 secs, 12,896,968 bytes)

-- λ> take 5 (show (suma2 40000))

-- "83669"

-- (1.82 secs, 4,549,612,264 bytes)

-- λ> take 5 (show (suma3 40000))

-- "83669"

-- (0.24 secs, 1,620,976,984 bytes)

El código se encuentra en GitHub.

Cálculo aproximado de integrales definidas

PorJosé A. Alonso 21-junio-202221-junio-2022

La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo usando la fórmula

   h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

1	h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

con a+n*h+h/2 <= b < a+(n+1)*h+h/2 y usando valores pequeños para h.

Definir la función

   integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

1	integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

tal que (integral a b f h) es el valor de dicha expresión. Por ejemplo, el cálculo de la integral de f(x) = x^3 entre 0 y 1, con paso 0.01, es

   integral 0 1 (^3) 0.01  ==  0.24998750000000042

1	integral 0 1 (^3) 0.01 == 0.24998750000000042

Otros ejemplos son

   integral 0 1 (^4) 0.01                   ==  0.19998333362500048
   integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01  ==  1.9999250000000026
   log 2 - integral 1 2 (\x -> 1/x) 0.01         ==  3.124931644782336e-6
   pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01  ==  -8.333333331389525e-6

integral 0 1 (^4) 0.01 == 0.19998333362500048

integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01 == 1.9999250000000026

log 2 - integral 1 2 (\x -> 1/x) 0.01 == 3.124931644782336e-6

pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01 == -8.333333331389525e-6

Soluciones

import Test.QuickCheck.HigherOrder (quickCheck')
import Test.QuickCheck (Property, (==>), quickCheck)

-- 1ª solución
-- ===========

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral1 a b f h
  | a+h/2 > b = 0
  | otherwise = h * f (a+h/2) + integral1 (a+h) b f h

-- 2ª solución
-- ===========

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral2 a b f h = aux a where
  aux x | x+h/2 > b = 0
        | otherwise = h * f (x+h/2) + aux (x+h)

-- 3ª solución
-- ===========

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral3 a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de
--    f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    suma 2 5 (1+) (^3)  ==  224
suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a
suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista
--    [a, s(a), s(s(a), ..., s(...(s(a))...)]
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    sucesion 3 20 (+2)  ==  [3,5,7,9,11,13,15,17,19]
sucesion :: Ord a => a -> a -> (a -> a) -> [a]
sucesion a b s = takeWhile (<=b) (iterate s a)

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property
prop_integral a b f h =
  a < b && h > 0 ==> 
  all (=~ integral1 a' b' f' h')
      [integral2 a' b' f' h',
       integral3 a' b' f' h']
  where
    a' = fromIntegral a
    b' = fromIntegral b
    h' = fromIntegral h
    f' = fromIntegral . f. round
    x =~ y = abs (x - y) < 0.001
        
-- La comprobación es
--    λ> quickCheck' prop_integral
--    +++ OK, passed 100 tests; 385 discarded.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> integral1 0 10 (^3) 0.00001
--    2499.999999881125
--    (2.63 secs, 1,491,006,744 bytes)
--    λ> integral2 0 10 (^3) 0.00001
--    2499.999999881125
--    (1.93 secs, 1,419,006,696 bytes)
--    λ> integral3 0 10 (^3) 0.00001
--    2499.9999998811422
--    (1.28 secs, 817,772,216 bytes)

import Test.QuickCheck.HigherOrder (quickCheck')

import Test.QuickCheck (Property, (==>), quickCheck)

-- 1ª solución

-- ===========

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral1 a b f h

| a+h/2 > b = 0

| otherwise = h * f (a+h/2) + integral1 (a+h) b f h

-- 2ª solución

-- ===========

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral2 a b f h = aux a where

aux x | x+h/2 > b = 0

| otherwise = h * f (x+h/2) + aux (x+h)

-- 3ª solución

-- ===========

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral3 a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de

-- f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- suma 2 5 (1+) (^3) == 224

suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a

suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista

-- [a, s(a), s(s(a), ..., s(...(s(a))...)]

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- sucesion 3 20 (+2) == [3,5,7,9,11,13,15,17,19]

sucesion :: Ord a => a -> a -> (a -> a) -> [a]

sucesion a b s = takeWhile (<=b) (iterate s a)

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property

prop_integral a b f h =

a < b && h > 0 ==>

all (=~ integral1 a' b' f' h')

[integral2 a' b' f' h',

integral3 a' b' f' h']

where

a' = fromIntegral a

b' = fromIntegral b

h' = fromIntegral h

f' = fromIntegral . f. round

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheck' prop_integral

-- +++ OK, passed 100 tests; 385 discarded.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> integral1 0 10 (^3) 0.00001

-- 2499.999999881125

-- (2.63 secs, 1,491,006,744 bytes)

-- λ> integral2 0 10 (^3) 0.00001

-- 2499.999999881125

-- (1.93 secs, 1,419,006,696 bytes)

-- λ> integral3 0 10 (^3) 0.00001

-- 2499.9999998811422

-- (1.28 secs, 817,772,216 bytes)

El código se encuentra en GitHub.

Menor número con una cantidad dada de divisores

PorJosé A. Alonso 20-junio-202221-junio-2022

El menor número con 2 divisores es el 2, ya que tiene 2 divisores (el 1 y el 2) y el anterior al 2 (el 1) sólo tiene 1 divisor (el 1).

El menor número con 4 divisores es el 6, ya que tiene 4 divisores (el 1, 2, 3 y 6) y sus anteriores (el 1, 2, 3, 4 y 5) tienen menos de 4 divisores (tienen 1, 1, 1, 3 y 1, respectivamente).

El menor número con 8 divisores es el 24, ya que tiene 8 divisores (el 1, 2, 3, 4, 6, 8, 12 y 24) y sus anteriores (del 1 al 23) tienen menos de 8 divisores.

El menor número con 16 divisores es el 120, ya que tiene 16 divisores (el 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 24, 30, 40, 60 y 120) y sus anteriores (del 1 al 119) tienen menos de 16 divisores.

Definir la función

   menor :: Integer -> Integer

1	menor :: Integer -> Integer

tal que (menor n) es el menor número con 2^n divisores. Por ejemplo,

   menor 1  ==  2
   menor 2  ==  6
   menor 3  ==  24
   menor 4  ==  120
   length (show (menor (4*10^4)))  ==  207945

menor 1 == 2

menor 2 == 6

menor 3 == 24

menor 4 == 120

length (show (menor (4*10^4))) == 207945

Comprobar con QuickCheck que, para todo k >=0, (menor (2^k)) es un divisor de (menor (2^(k+1))).

Nota: Este ejercicio está basado en el problema N1 de la Olimpíada Internacional de Matemáticas (IMO) del 2011.

Soluciones

import Data.List           (genericLength, genericTake, group)
import Data.Numbers.Primes (primeFactors, primes)
import Test.QuickCheck     (Positive (Positive), maxSize, stdArgs,
                            quickCheck, quickCheckWith)

-- 1ª solución
-- ===========

menor1 :: Integer -> Integer
menor1 n =
  head [x | x <- [1..], numeroDivisores x == 2^n]

-- (numeroDivisores n) es el número de divisores de n. Por ejemplo, 
--    numeroDivisores 12  ==  6
numeroDivisores :: Integer -> Integer
numeroDivisores =
  genericLength . divisores

-- (divisores x) es la lista de los divisores de x. Por ejemplo,
--    divisores 12  ==  [1,3,2,6,4,12]
--    divisores 25  ==  [1,5,25]
divisores :: Integer -> [Integer]
divisores n =
  [x | x <- [1..n], n `rem` x == 0]

-- 2ª solución
-- ===========

menor2 :: Integer -> Integer
menor2 n =
  head [x | x <- [1..], numeroDivisores2 x == 2^n]

numeroDivisores2 :: Integer -> Integer
numeroDivisores2 =
  product . map ((+1) . genericLength) . group . primeFactors

-- 3ª solución
-- ===========

menor3 :: Integer -> Integer
menor3 n = product (genericTake n potencias)

-- potencias es la sucesión de las potencias de la forma p^(2^k),
-- donde p es un número primo y k es un número natural, ordenadas de
-- menor a mayor. Por ejemplo,
--    take 14 potencias    ==  [2,3,4,5,7,9,11,13,16,17,19,23,25,29]
potencias :: [Integer]
potencias = 2 : mezcla (tail primes) (map (^2) potencias)

-- (mezcla xs ys) es la lista obtenida mezclando las dos listas xs e ys,
-- que se suponen ordenadas y disjuntas. Por ejemplo,
--    λ> take 15 (mezcla [2^n | n <- [1..]] [3^n | n <- [1..]])
--    [2,3,4,8,9,16,27,32,64,81,128,243,256,512,729]
mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)
                     | x > y = y : mezcla (x:xs) ys

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_menor :: Positive Integer -> Bool
prop_menor (Positive n) =
  all (== menor1 n)
      [menor2 n,
       menor3 n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=5}) prop_menor
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--   λ> menor 8
--   1081080
--   (47.69 secs, 94,764,856,352 bytes)
--   λ> menor2 8
--   1081080
--   (36.17 secs, 94,764,856,368 bytes)
--   λ> menor3 8
--   1081080
--   (0.00 secs, 116,960 bytes)

-- Definición de menor
-- ===================

-- En lo que sigue, usaremos menor3 como menor.
menor :: Integer -> Integer
menor = menor3

-- Propiedad
-- =========

-- La propiedad es
prop_menor_divide :: Positive Integer -> Bool
prop_menor_divide (Positive n) =
  menor (n+1) `mod` menor n == 0

-- La comprobación es
--    λ> quickCheck prop_menor_divide
--    +++ OK, passed 100 tests.

100

101

102

103

import Data.List (genericLength, genericTake, group)

import Data.Numbers.Primes (primeFactors, primes)

import Test.QuickCheck (Positive (Positive), maxSize, stdArgs,

quickCheck, quickCheckWith)

-- 1ª solución

-- ===========

menor1 :: Integer -> Integer

menor1 n =

head [x | x <- [1..], numeroDivisores x == 2^n]

-- (numeroDivisores n) es el número de divisores de n. Por ejemplo,

-- numeroDivisores 12 == 6

numeroDivisores :: Integer -> Integer

numeroDivisores =

genericLength . divisores

-- (divisores x) es la lista de los divisores de x. Por ejemplo,

-- divisores 12 == [1,3,2,6,4,12]

-- divisores 25 == [1,5,25]

divisores :: Integer -> [Integer]

divisores n =

[x | x <- [1..n], n `rem` x == 0]

-- 2ª solución

-- ===========

menor2 :: Integer -> Integer

menor2 n =

head [x | x <- [1..], numeroDivisores2 x == 2^n]

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

product . map ((+1) . genericLength) . group . primeFactors

-- 3ª solución

-- ===========

menor3 :: Integer -> Integer

menor3 n = product (genericTake n potencias)

-- potencias es la sucesión de las potencias de la forma p^(2^k),

-- donde p es un número primo y k es un número natural, ordenadas de

-- menor a mayor. Por ejemplo,

-- take 14 potencias == [2,3,4,5,7,9,11,13,16,17,19,23,25,29]

potencias :: [Integer]

potencias = 2 : mezcla (tail primes) (map (^2) potencias)

-- (mezcla xs ys) es la lista obtenida mezclando las dos listas xs e ys,

-- que se suponen ordenadas y disjuntas. Por ejemplo,

-- λ> take 15 (mezcla [2^n | n <- [1..]] [3^n | n <- [1..]])

-- [2,3,4,8,9,16,27,32,64,81,128,243,256,512,729]

mezcla :: Ord a => [a] -> [a] -> [a]

mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)

| x > y = y : mezcla (x:xs) ys

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_menor :: Positive Integer -> Bool

prop_menor (Positive n) =

all (== menor1 n)

[menor2 n,

menor3 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=5}) prop_menor

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> menor 8

-- 1081080

-- (47.69 secs, 94,764,856,352 bytes)

-- λ> menor2 8

-- 1081080

-- (36.17 secs, 94,764,856,368 bytes)

-- λ> menor3 8

-- 1081080

-- (0.00 secs, 116,960 bytes)

-- Definición de menor

-- ===================

-- En lo que sigue, usaremos menor3 como menor.

menor :: Integer -> Integer

menor = menor3

-- Propiedad

-- =========

-- La propiedad es

prop_menor_divide :: Positive Integer -> Bool

prop_menor_divide (Positive n) =

menor (n+1) `mod` menor n == 0

-- La comprobación es

-- λ> quickCheck prop_menor_divide

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Distancia esperada entre dos puntos de un cuadrado unitario

PorJosé A. Alonso 10-junio-20227-junio-2022

Definir, por simulación, la función

   distanciaEsperada :: Int -> IO Double

1	distanciaEsperada :: Int -> IO Double

tal que (distanciaEsperada n) es la distancia esperada entre n puntos del cuadrado unitario de vértices opuestos (0,0) y (1,1), elegidos aleatoriamente. Por ejemplo,

   distanciaEsperada 10     ==  0.43903617921423593
   distanciaEsperada 10     ==  0.6342350621260004
   distanciaEsperada 100    ==  0.5180418995364429
   distanciaEsperada 100    ==  0.5288261085653962
   distanciaEsperada 1000   ==  0.5143804432569616
   distanciaEsperada 10000  ==  0.5208360147922616

distanciaEsperada 10 == 0.43903617921423593

distanciaEsperada 10 == 0.6342350621260004

distanciaEsperada 100 == 0.5180418995364429

distanciaEsperada 100 == 0.5288261085653962

distanciaEsperada 1000 == 0.5143804432569616

distanciaEsperada 10000 == 0.5208360147922616

El valor exacto de la distancia esperada es

   ve = (sqrt(2) + 2 + 5*log(1+sqrt(2)))/15 = 0.5214054331647207

1	ve = (sqrt(2) + 2 + 5*log(1+sqrt(2)))/15 = 0.5214054331647207

Definir la función

   graficaDistanciaEsperada :: [Int] -> IO ()

1	graficaDistanciaEsperada :: [Int] -> IO ()

tal que (graficaDistanciaEsperada ns) dibuja las gráficas de los pares (n, distanciaEsperada n) para n en la lista creciente ns junto con la recta y = ve, donde ve es el valor exacto. Por ejemplo, (graficaDistanciaEsperada [10,30..4000]) dibuja

Soluciones

import Data.List     (genericLength)
import System.Random (newStdGen, randomRIO, randomRs)
import Control.Monad (replicateM)
import Graphics.Gnuplot.Simple

-- 1ª solución
-- ===========

-- Un punto es un par de números reales.
type Punto = (Double, Double)

-- (puntosDelCuadrado n) es una lista de n puntos del cuadrado
-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo, 
--    λ> puntosDelCuadrado 3
--    [(0.6067427807212623,0.24785843546479303),
--     (0.9579158098726746,8.047408846191773e-2),
--     (0.856758357789639,0.9814972717003113)]
--    λ> puntosDelCuadrado 3
--    [(1.9785720974027532e-2,0.6343219201012211),
--     (0.21903717179861604,0.20947986189590784),
--     (0.4739903340716357,1.2262474491489095e-2)]
puntosDelCuadrado :: Int -> IO [Punto]
puntosDelCuadrado n = do
  gen <- newStdGen
  let xs = randomRs (0,1) gen
      (as, ys) = splitAt n xs
      (bs, _)  = splitAt n ys
  return (zip as bs)
  
-- (distancia p1 p2) es la distancia entre los puntos p1 y p2. Por
-- ejemplo,
--    distancia (0,0) (3,4)  ==  5.0
distancia :: Punto -> Punto -> Double
distancia (x1,y1) (x2,y2) = sqrt ((x1-x2)^2+(y1-y2)^2)

-- (distancias ps) es la lista de las distancias entre los elementos 1º
-- y 2º, 3º y 4º, ... de ps. Por ejemplo,
--    distancias [(0,0),(3,4),(1,1),(7,9)]  ==  [5.0,10.0]
distancias :: [Punto] -> [Double]
distancias []         = []
distancias (p1:p2:ps) = distancia p1 p2 : distancias ps

-- (media xs) es la media aritmética de los elementos de xs. Por ejemplo,
--    media [1,7,1]  ==  3.0
media :: [Double] -> Double
media xs =
  sum xs / genericLength xs

-- (distanciaEsperada n) es la distancia esperada entre n puntos
-- aleatorios en el cuadrado unitario. Por ejemplo,
--    λ> distanciaEsperada 100
--    0.4712858421448363
--    λ> distanciaEsperada 100
--    0.4947745206856711
distanciaEsperada :: Int -> IO Double
distanciaEsperada n = do
  ps <- puntosDelCuadrado (2*n)
  return (media (distancias ps))

-- 2ª solución
-- ===========

distanciaEsperada2 :: Int -> IO Double
distanciaEsperada2 n = do
  ps <- puntosDelCuadrado2 (2*n)
  return (media (distancias ps))

-- (puntosDelCuadrado2 n) es una lista de n puntos del cuadrado
-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo, 
--    λ> puntosDelCuadrado2 3
--    [(0.9836699352638695,0.5143414844876929),
--     (0.8715237339877027,0.9905157772823782),
--     (0.29502946161912935,0.16889248111565192)]
--    λ> puntosDelCuadrado2 3
--    [(0.20405570457106392,0.47574116941605116),
--     (0.7128182811364226,3.201419787777959e-2),
--     (0.5576891231675457,0.9994474730919443)]
puntosDelCuadrado2 :: Int -> IO [Punto]
puntosDelCuadrado2 n =
  replicateM n puntoDelCuadrado2

-- (puntoDelCuadrado2 n) es un punto del cuadrado unitario de vértices
-- opuestos (0,0) y (1,1). Por ejemplo,  
--    λ> puntoDelCuadrado2
--    (0.7512991739803923,0.966436016138578)
--    λ> puntoDelCuadrado2
--    (0.7306826194847795,0.8984574498515252)
puntoDelCuadrado2 :: IO Punto
puntoDelCuadrado2 = do
  x <- randomRIO (0, 1.0)
  y <- randomRIO (0, 1.0)
  return (x, y)

-- 3ª solución
-- ===========

distanciaEsperada3 :: Int -> IO Double
distanciaEsperada3 n = do
  ds <- distanciasAleatorias n
  return (media ds)

-- (distanciasAleatorias n) es la lista de las distancias aleatorias
-- entre n pares de puntos del cuadrado unitario. Por ejemplo, 
--    λ> distanciasAleatorias 3
--    [0.8325589110989705,0.6803336613847881,0.1690051224111662]
--    λ> distanciasAleatorias 3
--    [0.3470124940889039,0.459002678562019,0.7665623634969365]
distanciasAleatorias :: Int -> IO [Double]
distanciasAleatorias n = 
  replicateM n distanciaAleatoria

-- distanciaAleatoria es la distancia de un par de punto del cuadrado
-- unitario elegidos aleatoriamente. Por ejemplo,
--    λ> distanciaAleatoria
--    0.8982361685460913
--    λ> distanciaAleatoria
--    0.9777207485571939
--    λ> distanciaAleatoria
--    0.6042223512347842
distanciaAleatoria :: IO Double
distanciaAleatoria = do 
  p1 <- puntoDelCuadrado2
  distancia p1 <$> puntoDelCuadrado2

-- 4ª solución
-- ===========

distanciaEsperada4 :: Int -> IO Double
distanciaEsperada4 n =
  media <$> distanciasAleatorias n

-- Gráfica
-- =======

graficaDistanciaEsperada :: [Int] -> IO ()
graficaDistanciaEsperada ns = do
  ys <- mapM distanciaEsperada ns
  let e = (sqrt 2 + 2 + 5 * log (1 + sqrt 2)) / 15
  plotLists [ Key Nothing
            -- , PNG "Distancia_esperada_entre_dos_puntos.png"
            ]
            [ zip ns ys
            , zip ns (repeat e)]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

import Data.List (genericLength)

import System.Random (newStdGen, randomRIO, randomRs)

import Control.Monad (replicateM)

import Graphics.Gnuplot.Simple

-- 1ª solución

-- ===========

-- Un punto es un par de números reales.

type Punto = (Double, Double)

-- (puntosDelCuadrado n) es una lista de n puntos del cuadrado

-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntosDelCuadrado 3

-- [(0.6067427807212623,0.24785843546479303),

-- (0.9579158098726746,8.047408846191773e-2),

-- (0.856758357789639,0.9814972717003113)]

-- λ> puntosDelCuadrado 3

-- [(1.9785720974027532e-2,0.6343219201012211),

-- (0.21903717179861604,0.20947986189590784),

-- (0.4739903340716357,1.2262474491489095e-2)]

puntosDelCuadrado :: Int -> IO [Punto]

puntosDelCuadrado n = do

gen <- newStdGen

let xs = randomRs (0,1) gen

(as, ys) = splitAt n xs

(bs, _) = splitAt n ys

return (zip as bs)

-- (distancia p1 p2) es la distancia entre los puntos p1 y p2. Por

-- ejemplo,

-- distancia (0,0) (3,4) == 5.0

distancia :: Punto -> Punto -> Double

distancia (x1,y1) (x2,y2) = sqrt ((x1-x2)^2+(y1-y2)^2)

-- (distancias ps) es la lista de las distancias entre los elementos 1º

-- y 2º, 3º y 4º, ... de ps. Por ejemplo,

-- distancias [(0,0),(3,4),(1,1),(7,9)] == [5.0,10.0]

distancias :: [Punto] -> [Double]

distancias [] = []

distancias (p1:p2:ps) = distancia p1 p2 : distancias ps

-- (media xs) es la media aritmética de los elementos de xs. Por ejemplo,

-- media [1,7,1] == 3.0

media :: [Double] -> Double

media xs =

sum xs / genericLength xs

-- (distanciaEsperada n) es la distancia esperada entre n puntos

-- aleatorios en el cuadrado unitario. Por ejemplo,

-- λ> distanciaEsperada 100

-- 0.4712858421448363

-- λ> distanciaEsperada 100

-- 0.4947745206856711

distanciaEsperada :: Int -> IO Double

distanciaEsperada n = do

ps <- puntosDelCuadrado (2*n)

return (media (distancias ps))

-- 2ª solución

-- ===========

distanciaEsperada2 :: Int -> IO Double

distanciaEsperada2 n = do

ps <- puntosDelCuadrado2 (2*n)

return (media (distancias ps))

-- (puntosDelCuadrado2 n) es una lista de n puntos del cuadrado

-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntosDelCuadrado2 3

-- [(0.9836699352638695,0.5143414844876929),

-- (0.8715237339877027,0.9905157772823782),

-- (0.29502946161912935,0.16889248111565192)]

-- λ> puntosDelCuadrado2 3

-- [(0.20405570457106392,0.47574116941605116),

-- (0.7128182811364226,3.201419787777959e-2),

-- (0.5576891231675457,0.9994474730919443)]

puntosDelCuadrado2 :: Int -> IO [Punto]

puntosDelCuadrado2 n =

replicateM n puntoDelCuadrado2

-- (puntoDelCuadrado2 n) es un punto del cuadrado unitario de vértices

-- opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntoDelCuadrado2

-- (0.7512991739803923,0.966436016138578)

-- λ> puntoDelCuadrado2

-- (0.7306826194847795,0.8984574498515252)

puntoDelCuadrado2 :: IO Punto

puntoDelCuadrado2 = do

x <- randomRIO (0, 1.0)

y <- randomRIO (0, 1.0)

return (x, y)

-- 3ª solución

-- ===========

distanciaEsperada3 :: Int -> IO Double

distanciaEsperada3 n = do

ds <- distanciasAleatorias n

return (media ds)

-- (distanciasAleatorias n) es la lista de las distancias aleatorias

-- entre n pares de puntos del cuadrado unitario. Por ejemplo,

-- λ> distanciasAleatorias 3

-- [0.8325589110989705,0.6803336613847881,0.1690051224111662]

-- λ> distanciasAleatorias 3

-- [0.3470124940889039,0.459002678562019,0.7665623634969365]

distanciasAleatorias :: Int -> IO [Double]

distanciasAleatorias n =

replicateM n distanciaAleatoria

-- distanciaAleatoria es la distancia de un par de punto del cuadrado

-- unitario elegidos aleatoriamente. Por ejemplo,

-- λ> distanciaAleatoria

-- 0.8982361685460913

-- λ> distanciaAleatoria

-- 0.9777207485571939

-- λ> distanciaAleatoria

-- 0.6042223512347842

distanciaAleatoria :: IO Double

distanciaAleatoria = do

p1 <- puntoDelCuadrado2

distancia p1 <$> puntoDelCuadrado2

-- 4ª solución

-- ===========

distanciaEsperada4 :: Int -> IO Double

distanciaEsperada4 n =

media <$> distanciasAleatorias n

-- Gráfica

-- =======

graficaDistanciaEsperada :: [Int] -> IO ()

graficaDistanciaEsperada ns = do

ys <- mapM distanciaEsperada ns

let e = (sqrt 2 + 2 + 5 * log (1 + sqrt 2)) / 15

plotLists [ Key Nothing

-- , PNG "Distancia_esperada_entre_dos_puntos.png"

]

[ zip ns ys

, zip ns (repeat e)]

El código se encuentra en GitHub.

Representación matricial de relaciones binarias

PorJosé A. Alonso 9-junio-20226-junio-2022

Dada una relación r sobre un conjunto de números enteros, la matriz asociada a r es una matriz booleana p (cuyos elementos son True o False), tal que p(i,j) = True si y sólo si i está relacionado con j mediante la relación r.

Las relaciones binarias homogéneas y las matrices booleanas se pueden representar por

   type Relacion = ([Int],[(Int,Int)])
   type Matriz = Array (Int,Int) Bool

1 2	type Relacion = ([Int],[(Int,Int)]) type Matriz = Array (Int,Int) Bool

Definir la función

   matrizRB:: Relacion -> Matriz

1	matrizRB:: Relacion -> Matriz

tal que (matrizRB r) es la matriz booleana asociada a r. Por ejemplo,

   λ> matrizRB ([1..3],[(1,1), (1,3), (3,1), (3,3)])
   array ((1,1),(3,3)) [((1,1),True) ,((1,2),False),((1,3),True),
                        ((2,1),False),((2,2),False),((2,3),False),
                        ((3,1),True) ,((3,2),False),((3,3),True)]
   λ> matrizRB ([1..3],[(1,3), (3,1)])
   array ((1,1),(3,3)) [((1,1),False),((1,2),False),((1,3),True),
                        ((2,1),False),((2,2),False),((2,3),False),
                        ((3,1),True) ,((3,2),False),((3,3),False)]
   λ> let n = 10^4 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)
   False

λ> matrizRB ([1..3],[(1,1), (1,3), (3,1), (3,3)])

array ((1,1),(3,3)) [((1,1),True) ,((1,2),False),((1,3),True),

((2,1),False),((2,2),False),((2,3),False),

((3,1),True) ,((3,2),False),((3,3),True)]

λ> matrizRB ([1..3],[(1,3), (3,1)])

array ((1,1),(3,3)) [((1,1),False),((1,2),False),((1,3),True),

((2,1),False),((2,2),False),((2,3),False),

((3,1),True) ,((3,2),False),((3,3),False)]

λ> let n = 10^4 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)

False

Soluciones

import Data.Array      (Array, accumArray, array, listArray)
import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, suchThat,
                        quickCheck)

type Relacion = ([Int],[(Int,Int)])
type Matriz   = Array (Int,Int) Bool

-- 1ª solución
-- ===========

matrizRB1 :: Relacion -> Matriz 
matrizRB1 r = 
    array ((1,1),(n,n)) 
          [((a,b), (a,b) `elem` grafo r) | a <- [1..n], b <- [1..n]]
    where n = maximum (universo r)
          universo (us,_) = us
          grafo (_,ps)    = ps

-- 2ª solución
-- ===========

matrizRB2 :: Relacion -> Matriz
matrizRB2 r = 
    listArray ((1,1),(n,n)) 
              [(a,b) `elem` snd r | a <- [1..n], b <- [1..n]]
    where n = maximum (fst r)

-- 3ª solución
-- ===========

matrizRB3 :: Relacion -> Matriz
matrizRB3 r = 
    accumArray (||) False ((1,1),(n,n)) (zip (snd r) (repeat True))
    where n = maximum (fst r)

-- Comprobación de equivalencia
-- ============================

-- Tipo de relaciones binarias
newtype RB = RB Relacion
  deriving Show

-- relacionArbitraria genera una relación arbitraria. Por ejemplo, 
--    λ> generate relacionArbitraria
--    RB ([1,2,3,4,5],[(1,4),(1,5),(2,3),(2,4),(4,2),(4,3),(4,4),(5,1),(5,2),(5,3),(5,4)])
relacionArbitraria :: Gen RB
relacionArbitraria = do
  n <- arbitrary `suchThat` (> 1)
  xs <- sublistOf [(x,y) | x <- [1..n], y <- [1..n]]
  return (RB ([1..n], xs))

-- RB es una subclase de Arbitrary
instance Arbitrary RB where
  arbitrary = relacionArbitraria

-- La propiedad es
prop_matrizRB :: RB -> Bool
prop_matrizRB (RB r) =
  all (== matrizRB1 r)
      [matrizRB2 r,
       matrizRB3 r]

-- La comprobación es
--    λ> quickCheck prop_matrixzB
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> let n = 2000 in matrizRB1 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (2.02 secs, 1,505,248,912 bytes)
--    λ> let n = 2000 in matrizRB2 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (1.92 secs, 833,232,360 bytes)
--    λ> let n = 2000 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (0.05 secs, 32,848,696 bytes)

import Data.Array (Array, accumArray, array, listArray)

import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, suchThat,

quickCheck)

type Relacion = ([Int],[(Int,Int)])

type Matriz = Array (Int,Int) Bool

-- 1ª solución

-- ===========

matrizRB1 :: Relacion -> Matriz

matrizRB1 r =

array ((1,1),(n,n))

[((a,b), (a,b) `elem` grafo r) | a <- [1..n], b <- [1..n]]

where n = maximum (universo r)

universo (us,_) = us

grafo (_,ps) = ps

-- 2ª solución

-- ===========

matrizRB2 :: Relacion -> Matriz

matrizRB2 r =

listArray ((1,1),(n,n))

[(a,b) `elem` snd r | a <- [1..n], b <- [1..n]]

where n = maximum (fst r)

-- 3ª solución

-- ===========

matrizRB3 :: Relacion -> Matriz

matrizRB3 r =

accumArray (||) False ((1,1),(n,n)) (zip (snd r) (repeat True))

where n = maximum (fst r)

-- Comprobación de equivalencia

-- ============================

-- Tipo de relaciones binarias

newtype RB = RB Relacion

deriving Show

-- relacionArbitraria genera una relación arbitraria. Por ejemplo,

-- λ> generate relacionArbitraria

-- RB ([1,2,3,4,5],[(1,4),(1,5),(2,3),(2,4),(4,2),(4,3),(4,4),(5,1),(5,2),(5,3),(5,4)])

relacionArbitraria :: Gen RB

relacionArbitraria = do

n <- arbitrary `suchThat` (> 1)

xs <- sublistOf [(x,y) | x <- [1..n], y <- [1..n]]

return (RB ([1..n], xs))

-- RB es una subclase de Arbitrary

instance Arbitrary RB where

arbitrary = relacionArbitraria

-- La propiedad es

prop_matrizRB :: RB -> Bool

prop_matrizRB (RB r) =

all (== matrizRB1 r)

[matrizRB2 r,

matrizRB3 r]

-- La comprobación es

-- λ> quickCheck prop_matrixzB

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> let n = 2000 in matrizRB1 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (2.02 secs, 1,505,248,912 bytes)

-- λ> let n = 2000 in matrizRB2 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (1.92 secs, 833,232,360 bytes)

-- λ> let n = 2000 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (0.05 secs, 32,848,696 bytes)

El código se encuentra en GitHub.

Codificación de Gödel

PorJosé A. Alonso 8-junio-20225-junio-2022

Dada una lista de números naturales xs, codificación de Gödel de xs se obtiene multiplicando las potencias de los primos sucesivos, siendo los exponentes los sucesores de los elementos de xs. Por ejemplo, si xs = [6,0,4], la codificación de xs es

   2^7 * 3^1 * 5^5 = 1200000

1	2^7 * 3^1 * 5^5 = 1200000

Definir las funciones

   codificaG   :: [Integer] -> Integer
   decodificaG :: Integer -> [Integer]

1 2	codificaG :: [Integer] -> Integer decodificaG :: Integer -> [Integer]

tales que

(codificaG xs) es la codificación de Gödel de xs. Por ejemplo,

     codificaG [6,0,4]            ==  1200000
     codificaG [3,1,1]            ==  3600
     codificaG [3,1,0,0,0,0,0,1]  ==  4423058640
     codificaG [1..6]             ==  126111168580452537982500

codificaG [6,0,4] == 1200000

codificaG [3,1,1] == 3600

codificaG [3,1,0,0,0,0,0,1] == 4423058640

codificaG [1..6] == 126111168580452537982500

(decodificaG n) es la lista xs cuya codificación es n. Por ejemplo,

     decodificaG 1200000                   ==  [6,0,4]
     decodificaG 3600                      ==  [3,1,1]
     decodificaG 4423058640                ==  [3,1,0,0,0,0,0,1]
     decodificaG 126111168580452537982500  ==  [1,2,3,4,5,6]

decodificaG 1200000 == [6,0,4]

decodificaG 3600 == [3,1,1]

decodificaG 4423058640 == [3,1,0,0,0,0,0,1]

decodificaG 126111168580452537982500 == [1,2,3,4,5,6]

Comprobar con QuickCheck que ambas funciones son inversas; es decir,

   decodificaG (codificaG xs) = xs
   codificaG (decodificaG n) = n

1 2	decodificaG (codificaG xs) = xs codificaG (decodificaG n) = n

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primes, primeFactors)
import Test.QuickCheck (NonNegative (getNonNegative),
                        NonEmptyList (NonEmpty),
                        Positive (Positive, getPositive), quickCheck)

-- 1ª definición de codificaG
-- ==========================

codificaG1 :: [Integer] -> Integer
codificaG1 = codificaG1' . codificaAux

codificaAux :: [Integer] -> [Integer]
codificaAux = map succ

codificaG1' :: [Integer] -> Integer
codificaG1' = aux primes
  where aux _ []          = 1
        aux (p:ps) (x:xs) = p^x * aux ps xs

-- 2ª definición de codificaG
-- ==========================

codificaG2 :: [Integer] -> Integer
codificaG2 = codificaG2' . codificaAux

codificaG2' :: [Integer] -> Integer
codificaG2' xs = product [p^x | (p, x) <- zip primes xs]

-- 3ª definición de codificaG
-- ==========================

codificaG3 :: [Integer] -> Integer
codificaG3 = codificaG3' . codificaAux

codificaG3' :: [Integer] -> Integer
codificaG3' xs = product (zipWith (^) primes xs)

-- 4ª definición de codificaG
-- ==========================

codificaG4 :: [Integer] -> Integer
codificaG4 = codificaG4' . codificaAux

codificaG4' :: [Integer] -> Integer
codificaG4' = product . zipWith (^) primes

-- Comprobación de equivalencia de codificaG
-- =========================================

-- La propiedad es
prop_codificaG :: [NonNegative Integer] -> Bool
prop_codificaG xs =
  all (== codificaG1 ys)
      [codificaG2 ys,
       codificaG3 ys,
       codificaG4 ys]
  where ys = map getNonNegative xs

-- La comprobación es
--    λ> quickCheck prop_codificaG
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia de codificaG
-- ======================================

-- La comparación es
--    λ> length (show (codificaG1 (replicate (4*10^4) 1)))
--    208100
--    (1.73 secs, 2,312,100,136 bytes)
--    λ> length (show (codificaG2 (replicate (4*10^4) 1)))
--    208100
--    (1.63 secs, 2,327,676,832 bytes)
--    λ> length (show (codificaG3 (replicate (4*10^4) 1)))
--    208100
--    (1.62 secs, 2,323,836,832 bytes)
--    λ> length (show (codificaG4 (replicate (4*10^4) 1)))
--    208100
--    (1.54 secs, 2,147,635,680 bytes)

-- Definición de codificaG
-- =======================

-- Usaremos la 4ª
codificaG :: [Integer] -> Integer
codificaG = codificaG4

-- Definición de decodificaG
-- =========================

decodificaG :: Integer -> [Integer]
decodificaG = decodificaAux . decodificaG'

decodificaAux :: [Integer] -> [Integer]
decodificaAux = map pred

decodificaG' :: Integer -> [Integer]
decodificaG' 1 = [0]
decodificaG' n = aux primes (group (primeFactors n))
  where aux _ [] = []
        aux (x:xs) ((y:ys):yss) | x == y    = 1 + genericLength ys : aux xs yss
                                | otherwise = 0 : aux xs ((y:ys):yss)

-- Comprobación de propiedades
-- ===========================

-- La primera propiedad es
prop_decodificaG_codificaG :: NonEmptyList (NonNegative Integer) -> Bool
prop_decodificaG_codificaG (NonEmpty xs) = 
  decodificaG (codificaG ys) == ys
  where ys = map getNonNegative xs

-- La comprobación es
--    λ> quickCheck prop_decodificaG_codificaG
--    +++ OK, passed 100 tests.

-- La 2ª propiedad es
prop_codificaG_decodificaG :: Positive Integer -> Bool
prop_codificaG_decodificaG (Positive n) = 
  codificaG (decodificaG n) == n

-- la comprobación es
--    λ> quickCheck prop_codificaG_decodificaG
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

import Data.List (genericLength, group)

import Data.Numbers.Primes (primes, primeFactors)

import Test.QuickCheck (NonNegative (getNonNegative),

NonEmptyList (NonEmpty),

Positive (Positive, getPositive), quickCheck)

-- 1ª definición de codificaG

-- ==========================

codificaG1 :: [Integer] -> Integer

codificaG1 = codificaG1' . codificaAux

codificaAux :: [Integer] -> [Integer]

codificaAux = map succ

codificaG1' :: [Integer] -> Integer

codificaG1' = aux primes

where aux _ [] = 1

aux (p:ps) (x:xs) = p^x * aux ps xs

-- 2ª definición de codificaG

-- ==========================

codificaG2 :: [Integer] -> Integer

codificaG2 = codificaG2' . codificaAux

codificaG2' :: [Integer] -> Integer

codificaG2' xs = product [p^x | (p, x) <- zip primes xs]

-- 3ª definición de codificaG

-- ==========================

codificaG3 :: [Integer] -> Integer

codificaG3 = codificaG3' . codificaAux

codificaG3' :: [Integer] -> Integer

codificaG3' xs = product (zipWith (^) primes xs)

-- 4ª definición de codificaG

-- ==========================

codificaG4 :: [Integer] -> Integer

codificaG4 = codificaG4' . codificaAux

codificaG4' :: [Integer] -> Integer

codificaG4' = product . zipWith (^) primes

-- Comprobación de equivalencia de codificaG

-- =========================================

-- La propiedad es

prop_codificaG :: [NonNegative Integer] -> Bool

prop_codificaG xs =

all (== codificaG1 ys)

[codificaG2 ys,

codificaG3 ys,

codificaG4 ys]

where ys = map getNonNegative xs

-- La comprobación es

-- λ> quickCheck prop_codificaG

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de codificaG

-- ======================================

-- La comparación es

-- λ> length (show (codificaG1 (replicate (4*10^4) 1)))

-- 208100

-- (1.73 secs, 2,312,100,136 bytes)

-- λ> length (show (codificaG2 (replicate (4*10^4) 1)))

-- 208100

-- (1.63 secs, 2,327,676,832 bytes)

-- λ> length (show (codificaG3 (replicate (4*10^4) 1)))

-- 208100

-- (1.62 secs, 2,323,836,832 bytes)

-- λ> length (show (codificaG4 (replicate (4*10^4) 1)))

-- 208100

-- (1.54 secs, 2,147,635,680 bytes)

-- Definición de codificaG

-- =======================

-- Usaremos la 4ª

codificaG :: [Integer] -> Integer

codificaG = codificaG4

-- Definición de decodificaG

-- =========================

decodificaG :: Integer -> [Integer]

decodificaG = decodificaAux . decodificaG'

decodificaAux :: [Integer] -> [Integer]

decodificaAux = map pred

decodificaG' :: Integer -> [Integer]

decodificaG' 1 = [0]

decodificaG' n = aux primes (group (primeFactors n))

where aux _ [] = []

aux (x:xs) ((y:ys):yss) | x == y = 1 + genericLength ys : aux xs yss

| otherwise = 0 : aux xs ((y:ys):yss)

-- Comprobación de propiedades

-- ===========================

-- La primera propiedad es

prop_decodificaG_codificaG :: NonEmptyList (NonNegative Integer) -> Bool

prop_decodificaG_codificaG (NonEmpty xs) =

decodificaG (codificaG ys) == ys

where ys = map getNonNegative xs

-- La comprobación es

-- λ> quickCheck prop_decodificaG_codificaG

-- +++ OK, passed 100 tests.

-- La 2ª propiedad es

prop_codificaG_decodificaG :: Positive Integer -> Bool

prop_codificaG_decodificaG (Positive n) =

codificaG (decodificaG n) == n

-- la comprobación es

-- λ> quickCheck prop_codificaG_decodificaG

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos circulares

PorJosé A. Alonso 7-junio-20223-junio-2022

Un primo circular es un número tal que todas las rotaciones de dígitos producen números primos. Por ejemplo, 195 es un primo circular ya que las rotaciones de sus dígitos son 197, 971 y 719 y los tres números son primos.

Definir la lista

   circulares :: [Integer]

1	circulares :: [Integer]

cuyo valor es la lista de los números primos circulares. Por ejemplo,

   take 16 circulares == [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,197]
   circulares !! 50   == 933199

1 2	take 16 circulares == [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,197] circulares !! 50 == 933199

Soluciones

import Data.Numbers.Primes (isPrime, primes)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========
  
circulares1 :: [Integer]
circulares1 = filter esCircular1 primes

-- (esCircular1 n) se verifica si n es un número circular. Por ejemplo, 
--    esCircular1 197  ==  True
--    esCircular1 157  ==  False
esCircular1 :: Integer -> Bool
esCircular1 = all (isPrime . read) . rotaciones1 . show 

-- (rotaciones1 xs) es la lista de las rotaciones obtenidas desplazando
-- el primer elemento xs al final. Por ejemplo,
--    rotaciones1 [2,3,5]  ==  [[2,3,5],[3,5,2],[5,2,3]]
rotaciones1 :: [a] -> [[a]]
rotaciones1 xs = reverse (aux (length xs) [xs])
    where aux 1 yss      = yss
          aux n (ys:yss) = aux (n-1) (rota ys : ys :yss)

-- 2ª solución
-- ===========
  
circulares2 :: [Integer]
circulares2 = filter esCircular2 primes

esCircular2 :: Integer -> Bool
esCircular2 = all (isPrime . read) . rotaciones2 . show 

rotaciones2 :: [a] -> [[a]]
rotaciones2 xs = take (length xs) (iterate rota xs)

-- (rota xs) es la lista añadiendo el primer elemento de xs al
-- final. Por ejemplo, 
--    rota [3,2,5,7]  ==  [2,5,7,3]
rota :: [a] -> [a]
rota (x:xs) = xs ++ [x]

-- 3ª solución
-- ===========

circulares3 :: [Integer]
circulares3 = filter (all isPrime . rotaciones3) primes 

rotaciones3 :: Integer -> [Integer]
rotaciones3 n = [read (take m (drop i (cycle s))) | i <- [1..m]]
    where s = show n
          m = length s

-- 4ª definición
-- =============

-- Nota. La 4ª definición es una mejora observando que para que n sea un
-- número primo circular es necesario que todos los dígitos de n sean
-- impares, salvo para n = 2. 

circulares4 :: [Integer]
circulares4 = 2 : filter esCircular4 primes

-- (esCircular4 n) se verifica si n es un número circular. Por ejemplo, 
--    esCircular4 197  ==  True
--    esCircular4 157  ==  False
esCircular4 :: Integer -> Bool
esCircular4 n = digitosImpares n && 
                all (isPrime . read) (rotaciones2 (show n))

-- (digitosImpares n) se verifica si todos los dígitos de n son
-- impares. Por ejemplo,
--    digitosImpares 7351  ==  True
--    digitosImpares 7341  ==  False
digitosImpares :: Integer -> Bool
digitosImpares = all (`elem` "135679") . show

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_circulares :: Positive Int -> Bool
prop_circulares (Positive n) =
  all (== circulares1 !! n)
      [circulares2 !! n,
       circulares3 !! n,
       circulares4 !! n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=50}) prop_circulares
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> circulares1 !! 46
--    331999
--    (2.08 secs, 7,229,208,200 bytes)
--    λ> circulares2 !! 46
--    331999
--    (1.93 secs, 7,165,043,992 bytes)
--    λ> circulares3 !! 46
--    331999
--    (0.74 secs, 2,469,098,648 bytes)
--    λ> circulares4 !! 46
--    331999
--    (0.28 secs, 917,501,600 bytes)

100

101

102

103

104

105

106

107

import Data.Numbers.Primes (isPrime, primes)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

circulares1 :: [Integer]

circulares1 = filter esCircular1 primes

-- (esCircular1 n) se verifica si n es un número circular. Por ejemplo,

-- esCircular1 197 == True

-- esCircular1 157 == False

esCircular1 :: Integer -> Bool

esCircular1 = all (isPrime . read) . rotaciones1 . show

-- (rotaciones1 xs) es la lista de las rotaciones obtenidas desplazando

-- el primer elemento xs al final. Por ejemplo,

-- rotaciones1 [2,3,5] == [[2,3,5],[3,5,2],[5,2,3]]

rotaciones1 :: [a] -> [[a]]

rotaciones1 xs = reverse (aux (length xs) [xs])

where aux 1 yss = yss

aux n (ys:yss) = aux (n-1) (rota ys : ys :yss)

-- 2ª solución

-- ===========

circulares2 :: [Integer]

circulares2 = filter esCircular2 primes

esCircular2 :: Integer -> Bool

esCircular2 = all (isPrime . read) . rotaciones2 . show

rotaciones2 :: [a] -> [[a]]

rotaciones2 xs = take (length xs) (iterate rota xs)

-- (rota xs) es la lista añadiendo el primer elemento de xs al

-- final. Por ejemplo,

-- rota [3,2,5,7] == [2,5,7,3]

rota :: [a] -> [a]

rota (x:xs) = xs ++ [x]

-- 3ª solución

-- ===========

circulares3 :: [Integer]

circulares3 = filter (all isPrime . rotaciones3) primes

rotaciones3 :: Integer -> [Integer]

rotaciones3 n = [read (take m (drop i (cycle s))) | i <- [1..m]]

where s = show n

m = length s

-- 4ª definición

-- =============

-- Nota. La 4ª definición es una mejora observando que para que n sea un

-- número primo circular es necesario que todos los dígitos de n sean

-- impares, salvo para n = 2.

circulares4 :: [Integer]

circulares4 = 2 : filter esCircular4 primes

-- (esCircular4 n) se verifica si n es un número circular. Por ejemplo,

-- esCircular4 197 == True

-- esCircular4 157 == False

esCircular4 :: Integer -> Bool

esCircular4 n = digitosImpares n &&

all (isPrime . read) (rotaciones2 (show n))

-- (digitosImpares n) se verifica si todos los dígitos de n son

-- impares. Por ejemplo,

-- digitosImpares 7351 == True

-- digitosImpares 7341 == False

digitosImpares :: Integer -> Bool

digitosImpares = all (`elem` "135679") . show

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_circulares :: Positive Int -> Bool

prop_circulares (Positive n) =

all (== circulares1 !! n)

[circulares2 !! n,

circulares3 !! n,

circulares4 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=50}) prop_circulares

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> circulares1 !! 46

-- 331999

-- (2.08 secs, 7,229,208,200 bytes)

-- λ> circulares2 !! 46

-- 331999

-- (1.93 secs, 7,165,043,992 bytes)

-- λ> circulares3 !! 46

-- 331999

-- (0.74 secs, 2,469,098,648 bytes)

-- λ> circulares4 !! 46

-- 331999

-- (0.28 secs, 917,501,600 bytes)

El código se encuentra en GitHub.

Diccionario de frecuencias

PorJosé A. Alonso 6-junio-20222-junio-2022

Definir la función

   frecuencias :: Ord a => [a] -> Map a Int

1	frecuencias :: Ord a => [a] -> Map a Int

tal que (frecuencias xs) es el diccionario formado por los elementos de xs junto con el número de veces que aparecen en xs. Por ejemplo,

   λ> frecuencias "sosos"
   fromList [('o',2),('s',3)]
   λ> frecuencias (show (10^100))
   fromList [('0',100),('1',1)]
   λ> frecuencias (take (10^6) (cycle "abc"))
   fromList [('a',333334),('b',333333),('c',333333)]
   λ> size (frecuencias (take (10^6) (cycle [1..10^6])))
   1000000

λ> frecuencias "sosos"

fromList [('o',2),('s',3)]

λ> frecuencias (show (10^100))

fromList [('0',100),('1',1)]

λ> frecuencias (take (10^6) (cycle "abc"))

fromList [('a',333334),('b',333333),('c',333333)]

λ> size (frecuencias (take (10^6) (cycle [1..10^6])))

1000000

Soluciones

import Data.List (foldl')
import Data.Map (Map, empty, insertWith, fromList, fromListWith, size)
import Test.QuickCheck

-- 1ª solución
-- ===========

frecuencias1 :: Ord a => [a] -> Map a Int
frecuencias1 []     = empty
frecuencias1 (x:xs) = insertWith (+) x 1 (frecuencias1 xs)

-- 2ª solución
-- ===========

frecuencias2 :: Ord a => [a] -> Map a Int
frecuencias2 = foldl' (\d x-> insertWith (+) x 1 d) empty

-- 3ª solución
-- ===========

frecuencias3 :: Ord a => [a] -> Map a Int
frecuencias3 xs = fromListWith (+) (zip xs (repeat 1))

-- Equivalencia de las definiciones
-- ================================

-- La propiedad es
prop_frecuencias :: [Int] -> Bool
prop_frecuencias xs =
  all (== frecuencias1 xs)
      [ frecuencias2 xs
      , frecuencias3 xs]

-- La comprobación es
--    λ> quickCheck prop_frecuencias
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> frecuencias1 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.89 secs, 453,842,448 bytes)
--    λ> frecuencias2 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.54 secs, 274,181,128 bytes)
--    λ> frecuencias3 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.29 secs, 313,787,976 bytes)
     
--    λ> size (frecuencias1 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (3.76 secs, 2,651,926,024 bytes)
--    λ> size (frecuencias2 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (1.03 secs, 1,640,678,448 bytes)
--    λ> size (frecuencias3 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (0.88 secs, 1,672,678,536 bytes)

import Data.List (foldl')

import Data.Map (Map, empty, insertWith, fromList, fromListWith, size)

import Test.QuickCheck

-- 1ª solución

-- ===========

frecuencias1 :: Ord a => [a] -> Map a Int

frecuencias1 [] = empty

frecuencias1 (x:xs) = insertWith (+) x 1 (frecuencias1 xs)

-- 2ª solución

-- ===========

frecuencias2 :: Ord a => [a] -> Map a Int

frecuencias2 = foldl' (\d x-> insertWith (+) x 1 d) empty

-- 3ª solución

-- ===========

frecuencias3 :: Ord a => [a] -> Map a Int

frecuencias3 xs = fromListWith (+) (zip xs (repeat 1))

-- Equivalencia de las definiciones

-- ================================

-- La propiedad es

prop_frecuencias :: [Int] -> Bool

prop_frecuencias xs =

all (== frecuencias1 xs)

[ frecuencias2 xs

, frecuencias3 xs]

-- La comprobación es

-- λ> quickCheck prop_frecuencias

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> frecuencias1 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.89 secs, 453,842,448 bytes)

-- λ> frecuencias2 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.54 secs, 274,181,128 bytes)

-- λ> frecuencias3 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.29 secs, 313,787,976 bytes)

-- λ> size (frecuencias1 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (3.76 secs, 2,651,926,024 bytes)

-- λ> size (frecuencias2 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (1.03 secs, 1,640,678,448 bytes)

-- λ> size (frecuencias3 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (0.88 secs, 1,672,678,536 bytes)

El código se encuentra en GitHub.

Descomposiciones con sumandos 1 ó 2

PorJosé A. Alonso 3-junio-20225-junio-2022

Definir la funciones

   sumas  :: Int -> [[Int]]
   nSumas :: Int -> Integer

1 2	sumas :: Int -> [[Int]] nSumas :: Int -> Integer

tales que

(sumas n) es la lista de las descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por ejemplo,

      sumas 1            ==  [[1]]
      sumas 2            ==  [[1,1],[2]]
      sumas 3            ==  [[1,1,1],[1,2],[2,1]]
      sumas 4            ==  [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]
      length (sumas 26)  ==  196418
      length (sumas 33)  ==  5702887

sumas 1 == [[1]]

sumas 2 == [[1,1],[2]]

sumas 3 == [[1,1,1],[1,2],[2,1]]

sumas 4 == [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]

length (sumas 26) == 196418

length (sumas 33) == 5702887

(nSumas n) es el número de descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por ejemplo,

      nSumas 4                      ==  5
      nSumas 123                    ==  36726740705505779255899443
      length (show (nSumas 123456)) ==  25801

nSumas 4 == 5

nSumas 123 == 36726740705505779255899443

length (show (nSumas 123456)) == 25801

Soluciones

import Data.List  (genericIndex, genericLength)
import Data.Array ((!), array)
import Test.QuickCheck (Positive(Positive), quickCheckWith)

-- 1ª solución de sumas
-- ====================

sumas1 :: Int -> [[Int]]
sumas1 0 = [[]]
sumas1 1 = [[1]]
sumas1 n = [1:xs | xs <- sumas1 (n-1)] ++ [2:xs | xs <- sumas1 (n-2)]

-- 2ª solución de sumas
-- ====================

sumas2 :: Int -> [[Int]]
sumas2 0 = [[]]
sumas2 1 = [[1]]
sumas2 n = map (1:) (sumas2 (n-1)) ++ map (2:) (sumas2 (n-2))

-- 3ª solución de sumas
-- ====================

sumas3 :: Int -> [[Int]]
sumas3 n = v ! n
  where v = array (0,n) [(i, f i) | i <- [0..n]]
        f 0 = [[]]
        f 1 = [[1]]
        f k = map (1:) (v!(k-1)) ++ map (2:) (v!(k-2))
 
-- 4ª solución de sumas
-- ====================

sumas4 :: Int -> [[Int]]
sumas4 n = sucSumas !! n

-- sucSumas es la sucesión cuyo n-ésimo elemento es la lista de las
-- descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por
-- ejemplo,
--    λ> take 4 sucSumas
--    [[[]],[[1]],[[1,1],[2]],[[1,1,1],[1,2],[2,1]]]
--    λ> mapM_ print (take 5 sucSumas)
--    [[]]
--    [[1]]
--    [[1,1],[2]]
--    [[1,1,1],[1,2],[2,1]]
--    [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]
sucSumas :: [[[Int]]]
sucSumas = [[]] : [[1]] : zipWith f (tail sucSumas) sucSumas
  where f xs ys = map (1:) xs ++ map (2:) ys

-- Comprobación de equivalencia de sumas
-- =====================================

-- La propiedad es
prop_sumas :: Positive Int -> Bool
prop_sumas (Positive n) =
  all (== sumas1 n)
      [sumas2 n,
       sumas3 n,
       sumas4 n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_sumas
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumas
-- ==================================

-- La comparación es
--    λ> length (sumas1 28)
--    514229
--    (2.79 secs, 1,739,784,512 bytes)
--    λ> length (sumas2 28)
--    514229
--    (1.33 secs, 1,512,291,248 bytes)
--    λ> length (sumas3 28)
--    514229
--    (0.20 secs, 165,215,800 bytes)
--    λ> length (sumas4 28)
--    514229
--    (0.17 secs, 165,201,592 bytes)
--
--    λ> length (sumas3 33)
--    5702887
--    (2.16 secs, 1,830,761,864 bytes)
--    λ> length (sumas4 33)
--    5702887
--    (1.44 secs, 1,830,749,832 bytes)

-- Definición de sumas
-- ===================

-- La cuarta solución es más eficiente y es la que usaremos en lo
-- sucesivo:
sumas :: Int -> [[Int]]
sumas = sumas4

-- 1ª solución de nSumas
-- =====================

nSumas1 :: Int -> Integer
nSumas1 = genericLength . sumas2

-- 2ª solución de nSumas
-- =====================

nSumas2 :: Int -> Integer
nSumas2 0 = 1
nSumas2 1 = 1
nSumas2 n = nSumas2 (n-1) + nSumas2 (n-2)

-- 3ª solución de nSumas
-- =====================

nSumas3 :: Int -> Integer
nSumas3 n = v ! n
  where v = array (0,n) [(i,f i) | i <- [0..n]]
        f 0 = 1
        f 1 = 1
        f k = v ! (k-1) + v ! (k-2)

-- 4ª solución de nSumas
-- =====================

nSumas4 :: Int -> Integer
nSumas4 n = aux `genericIndex` n
  where aux = 1 : 1 : zipWith (+) aux (tail aux) 

-- Comprobación de equivalencia de nSumas
-- ======================================

-- La propiedad es
prop_nSumas :: Positive Int -> Bool
prop_nSumas (Positive n) =
  all (== nSumas1 n)
      [nSumas2 n,
       nSumas3 n,
       nSumas4 n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_nSumas
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia de nSumas
-- ===================================

-- La comparación es
--    λ> nSumas1 33
--    5702887
--    (17.32 secs, 23,140,562,600 bytes)
--    λ> nSumas2 33
--    5702887
--    (3.48 secs, 1,870,676,904 bytes)
--    λ> nSumas3 33
--    5702887
--    (0.00 secs, 152,960 bytes)
--    λ> nSumas4 33
--    5702887
--    (0.00 secs, 139,456 bytes)
--    
--    λ> length (show (nSumas3 (2*10^5)))
--    41798
--    (1.41 secs, 1,895,295,528 bytes)
--    λ> length (show (nSumas4 (2*10^5)))
--    41798
--    (2.39 secs, 1,834,998,800 bytes)

-- Nota. El valor de (nSumas n) es el n-ésimo término de la sucesión de
-- Fibonacci 1, 1, 2, 3, 5, 8, ...

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

import Data.List (genericIndex, genericLength)

import Data.Array ((!), array)

import Test.QuickCheck (Positive(Positive), quickCheckWith)

-- 1ª solución de sumas

-- ====================

sumas1 :: Int -> [[Int]]

sumas1 0 = [[]]

sumas1 1 = [[1]]

sumas1 n = [1:xs | xs <- sumas1 (n-1)] ++ [2:xs | xs <- sumas1 (n-2)]

-- 2ª solución de sumas

-- ====================

sumas2 :: Int -> [[Int]]

sumas2 0 = [[]]

sumas2 1 = [[1]]

sumas2 n = map (1:) (sumas2 (n-1)) ++ map (2:) (sumas2 (n-2))

-- 3ª solución de sumas

-- ====================

sumas3 :: Int -> [[Int]]

sumas3 n = v ! n

where v = array (0,n) [(i, f i) | i <- [0..n]]

f 0 = [[]]

f 1 = [[1]]

f k = map (1:) (v!(k-1)) ++ map (2:) (v!(k-2))

-- 4ª solución de sumas

-- ====================

sumas4 :: Int -> [[Int]]

sumas4 n = sucSumas !! n

-- sucSumas es la sucesión cuyo n-ésimo elemento es la lista de las

-- descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por

-- ejemplo,

-- λ> take 4 sucSumas

-- [[[]],[[1]],[[1,1],[2]],[[1,1,1],[1,2],[2,1]]]

-- λ> mapM_ print (take 5 sucSumas)

-- [[]]

-- [[1]]

-- [[1,1],[2]]

-- [[1,1,1],[1,2],[2,1]]

-- [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]

sucSumas :: [[[Int]]]

sucSumas = [[]] : [[1]] : zipWith f (tail sucSumas) sucSumas

where f xs ys = map (1:) xs ++ map (2:) ys

-- Comprobación de equivalencia de sumas

-- =====================================

-- La propiedad es

prop_sumas :: Positive Int -> Bool

prop_sumas (Positive n) =

all (== sumas1 n)

[sumas2 n,

sumas3 n,

sumas4 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_sumas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumas

-- ==================================

-- La comparación es

-- λ> length (sumas1 28)

-- 514229

-- (2.79 secs, 1,739,784,512 bytes)

-- λ> length (sumas2 28)

-- 514229

-- (1.33 secs, 1,512,291,248 bytes)

-- λ> length (sumas3 28)

-- 514229

-- (0.20 secs, 165,215,800 bytes)

-- λ> length (sumas4 28)

-- 514229

-- (0.17 secs, 165,201,592 bytes)

-- λ> length (sumas3 33)

-- 5702887

-- (2.16 secs, 1,830,761,864 bytes)

-- λ> length (sumas4 33)

-- 5702887

-- (1.44 secs, 1,830,749,832 bytes)

-- Definición de sumas

-- ===================

-- La cuarta solución es más eficiente y es la que usaremos en lo

-- sucesivo:

sumas :: Int -> [[Int]]

sumas = sumas4

-- 1ª solución de nSumas

-- =====================

nSumas1 :: Int -> Integer

nSumas1 = genericLength . sumas2

-- 2ª solución de nSumas

-- =====================

nSumas2 :: Int -> Integer

nSumas2 0 = 1

nSumas2 1 = 1

nSumas2 n = nSumas2 (n-1) + nSumas2 (n-2)

-- 3ª solución de nSumas

-- =====================

nSumas3 :: Int -> Integer

nSumas3 n = v ! n

where v = array (0,n) [(i,f i) | i <- [0..n]]

f 0 = 1

f 1 = 1

f k = v ! (k-1) + v ! (k-2)

-- 4ª solución de nSumas

-- =====================

nSumas4 :: Int -> Integer

nSumas4 n = aux `genericIndex` n

where aux = 1 : 1 : zipWith (+) aux (tail aux)

-- Comprobación de equivalencia de nSumas

-- ======================================

-- La propiedad es

prop_nSumas :: Positive Int -> Bool

prop_nSumas (Positive n) =

all (== nSumas1 n)

[nSumas2 n,

nSumas3 n,

nSumas4 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_nSumas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de nSumas

-- ===================================

-- La comparación es

-- λ> nSumas1 33

-- 5702887

-- (17.32 secs, 23,140,562,600 bytes)

-- λ> nSumas2 33

-- 5702887

-- (3.48 secs, 1,870,676,904 bytes)

-- λ> nSumas3 33

-- 5702887

-- (0.00 secs, 152,960 bytes)

-- λ> nSumas4 33

-- 5702887

-- (0.00 secs, 139,456 bytes)

-- λ> length (show (nSumas3 (2*10^5)))

-- 41798

-- (1.41 secs, 1,895,295,528 bytes)

-- λ> length (show (nSumas4 (2*10^5)))

-- 41798

-- (2.39 secs, 1,834,998,800 bytes)

-- Nota. El valor de (nSumas n) es el n-ésimo término de la sucesión de

-- Fibonacci 1, 1, 2, 3, 5, 8, ...

El código se encuentra en GitHub.

Suma de los elementos de las diagonales matrices espirales

PorJosé A. Alonso 2-junio-202231-mayo-2022

Empezando con el número 1 y moviéndose en el sentido de las agujas del reloj se obtienen las matrices espirales

   |1 2|   |7 8 9|   | 7  8  9 10|   |21 22 23 24 25|
   |4 3|   |6 1 2|   | 6  1  2 11|   |20  7  8  9 10|
           |5 4 3|   | 5  4  3 12|   |19  6  1  2 11|
                     |16 15 14 13|   |18  5  4  3 12|
                                     |17 16 15 14 13|

|1 2| |7 8 9| | 7 8 9 10| |21 22 23 24 25|

|4 3| |6 1 2| | 6 1 2 11| |20 7 8 9 10|

|5 4 3| | 5 4 3 12| |19 6 1 2 11|

|16 15 14 13| |18 5 4 3 12|

|17 16 15 14 13|

La suma los elementos de sus diagonales es

   + en la 2x2: 1+3+2+4               =  10
   + en la 3x3: 1+3+5+7+9             =  25
   + en la 4x4: 1+2+3+4+7+10+13+16    =  56
   + en la 5x5: 1+3+5+7+9+13+17+21+25 = 101

+ en la 2x2: 1+3+2+4 = 10

+ en la 3x3: 1+3+5+7+9 = 25

+ en la 4x4: 1+2+3+4+7+10+13+16 = 56

+ en la 5x5: 1+3+5+7+9+13+17+21+25 = 101

Definir la función

   sumaDiagonales :: Integer -> Integer

1	sumaDiagonales :: Integer -> Integer

tal que (sumaDiagonales n) es la suma de los elementos en las diagonales de la matriz espiral de orden nxn. Por ejemplo.

   sumaDiagonales 1         ==  1
   sumaDiagonales 2         ==  10
   sumaDiagonales 3         ==  25
   sumaDiagonales 4         ==  56
   sumaDiagonales 5         ==  101
   sumaDiagonales (10^6)    ==  666667166668000000
   sumaDiagonales (1+10^6)  ==  666669166671000001

   sumaDiagonales (10^2)  ==         671800
   sumaDiagonales (10^3)  ==        667168000
   sumaDiagonales (10^4)  ==       666716680000
   sumaDiagonales (10^5)  ==      666671666800000
   sumaDiagonales (10^6)  ==     666667166668000000
   sumaDiagonales (10^7)  ==    666666716666680000000
   sumaDiagonales (10^8)  ==   666666671666666800000000
   sumaDiagonales (10^9)  ==  666666667166666668000000000

sumaDiagonales 1 == 1

sumaDiagonales 2 == 10

sumaDiagonales 3 == 25

sumaDiagonales 4 == 56

sumaDiagonales 5 == 101

sumaDiagonales (10^6) == 666667166668000000

sumaDiagonales (1+10^6) == 666669166671000001

sumaDiagonales (10^2) == 671800

sumaDiagonales (10^3) == 667168000

sumaDiagonales (10^4) == 666716680000

sumaDiagonales (10^5) == 666671666800000

sumaDiagonales (10^6) == 666667166668000000

sumaDiagonales (10^7) == 666666716666680000000

sumaDiagonales (10^8) == 666666671666666800000000

sumaDiagonales (10^9) == 666666667166666668000000000

Comprobar con QuickCheck que el último dígito de (sumaDiagonales n) es 0, 4 ó 6 si n es par y es 1, 5 ó 7 en caso contrario.

Soluciones

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

sumaDiagonales1 :: Integer -> Integer
sumaDiagonales1 = sum . elementosEnDiagonales

-- (elementosEnDiagonales n) es la lista de los elementos en las
-- diagonales de la matriz espiral de orden nxn. Por ejemplo,
--    elementosEnDiagonales 1  ==  [1]
--    elementosEnDiagonales 2  ==  [1,2,3,4]
--    elementosEnDiagonales 3  ==  [1,3,5,7,9]
--    elementosEnDiagonales 4  ==  [1,2,3,4,7,10,13,16]
--    elementosEnDiagonales 5  ==  [1,3,5,7,9,13,17,21,25]
elementosEnDiagonales :: Integer -> [Integer]
elementosEnDiagonales n 
  | even n    = tail (scanl (+) 0 (concatMap (replicate 4) [1,3..n-1]))
  | otherwise = scanl (+) 1 (concatMap (replicate 4) [2,4..n-1])

-- 2ª solución
-- ===========

sumaDiagonales2 :: Integer -> Integer
sumaDiagonales2 n
  | even n    = (-1) + n `div` 2 + sum [2*k^2-k+1 | k <- [0..n]]
  | otherwise = 1 + sum [4*k^2-6*k+6 | k <- [3,5..n]]

-- 3ª solución
-- ===========

sumaDiagonales3 :: Integer -> Integer
sumaDiagonales3 n
  | even n    = n * (4*n^2 + 3*n + 8) `div` 6
  | otherwise = (4*n^3 + 3*n^2 + 8*n - 9) `div` 6

-- Equivalencia de las definiciones
-- ================================

-- La propiedad es
prop_sumaDiagonales :: Positive Integer -> Bool
prop_sumaDiagonales (Positive n) =
  all (== sumaDiagonales1 n)
      [sumaDiagonales2 n,
       sumaDiagonales3 n]

-- La comprobación es
--    λ> quickCheck prop_sumaDiagonales_equiv
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> sumaDiagonales (2*10^6)
--    5333335333336000000
--    (2.30 secs, 1,521,955,848 bytes)
--    λ> sumaDiagonales2 (2*10^6)
--    5333335333336000000
--    (2.77 secs, 1,971,411,440 bytes)
--    λ> sumaDiagonales3 (2*10^6)
--    5333335333336000000
--    (0.01 secs, 139,520 bytes)

-- Propiedad
-- =========

-- La propiedad es
prop_sumaDiagonales2 :: Positive Integer -> Bool
prop_sumaDiagonales2 (Positive n) 
  | even n    = x `elem` [0,4,6] 
  | otherwise = x `elem` [1,5,7] 
  where x = sumaDiagonales1 n `mod` 10

-- La comprobación es
--    λ> quickCheck prop_sumaDiagonales2
--    +++ OK, passed 100 tests.

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

sumaDiagonales1 :: Integer -> Integer

sumaDiagonales1 = sum . elementosEnDiagonales

-- (elementosEnDiagonales n) es la lista de los elementos en las

-- diagonales de la matriz espiral de orden nxn. Por ejemplo,

-- elementosEnDiagonales 1 == [1]

-- elementosEnDiagonales 2 == [1,2,3,4]

-- elementosEnDiagonales 3 == [1,3,5,7,9]

-- elementosEnDiagonales 4 == [1,2,3,4,7,10,13,16]

-- elementosEnDiagonales 5 == [1,3,5,7,9,13,17,21,25]

elementosEnDiagonales :: Integer -> [Integer]

elementosEnDiagonales n

| even n = tail (scanl (+) 0 (concatMap (replicate 4) [1,3..n-1]))

| otherwise = scanl (+) 1 (concatMap (replicate 4) [2,4..n-1])

-- 2ª solución

-- ===========

sumaDiagonales2 :: Integer -> Integer

sumaDiagonales2 n

| even n = (-1) + n `div` 2 + sum [2*k^2-k+1 | k <- [0..n]]

| otherwise = 1 + sum [4*k^2-6*k+6 | k <- [3,5..n]]

-- 3ª solución

-- ===========

sumaDiagonales3 :: Integer -> Integer

sumaDiagonales3 n

| even n = n * (4*n^2 + 3*n + 8) `div` 6

| otherwise = (4*n^3 + 3*n^2 + 8*n - 9) `div` 6

-- Equivalencia de las definiciones

-- ================================

-- La propiedad es

prop_sumaDiagonales :: Positive Integer -> Bool

prop_sumaDiagonales (Positive n) =

all (== sumaDiagonales1 n)

[sumaDiagonales2 n,

sumaDiagonales3 n]

-- La comprobación es

-- λ> quickCheck prop_sumaDiagonales_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> sumaDiagonales (2*10^6)

-- 5333335333336000000

-- (2.30 secs, 1,521,955,848 bytes)

-- λ> sumaDiagonales2 (2*10^6)

-- 5333335333336000000

-- (2.77 secs, 1,971,411,440 bytes)

-- λ> sumaDiagonales3 (2*10^6)

-- 5333335333336000000

-- (0.01 secs, 139,520 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_sumaDiagonales2 :: Positive Integer -> Bool

prop_sumaDiagonales2 (Positive n)

| even n = x `elem` [0,4,6]

| otherwise = x `elem` [1,5,7]

where x = sumaDiagonales1 n `mod` 10

-- La comprobación es

-- λ> quickCheck prop_sumaDiagonales2

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Término ausente en una progresión aritmética

PorJosé A. Alonso 1-junio-202230-mayo-2022

Una progresión aritmética es una sucesión de números tales que la diferencia de dos términos sucesivos cualesquiera de la sucesión es constante.

Definir la función

   ausente :: Integral a => [a] -> a

1	ausente :: Integral a => [a] -> a

tal que (ausente xs) es el único término ausente de la progresión aritmética xs. Por ejemplo,

   ausente [3,7,9,11]               ==  5
   ausente [3,5,9,11]               ==  7
   ausente [3,5,7,11]               ==  9
   ausente ([1..9]++[11..])         ==  10
   ausente ([1..10^6] ++ [2+10^6])  ==  1000001

ausente [3,7,9,11] == 5

ausente [3,5,9,11] == 7

ausente [3,5,7,11] == 9

ausente ([1..9]++[11..]) == 10

ausente ([1..10^6] ++ [2+10^6]) == 1000001

Nota. Se supone que la lista tiene al menos 3 elementos, que puede ser infinita y que sólo hay un término de la progresión aritmética que no está en la lista.

Soluciones

import Data.List (group, genericLength)
import Test.QuickCheck (Arbitrary, Gen,
                        arbitrary, frequency, suchThat, quickCheck) 

-- 1ª solución
-- ===========

ausente1 :: Integral a => [a] -> a
ausente1 (x1:xs@(x2:x3:_))
  | d1 == d2     = ausente1 xs
  | d1 == 2 * d2 = x1 + d2
  | d2 == 2 * d1 = x2 + d1
  where d1 = x2 - x1
        d2 = x3 - x2          

-- 2ª solución
-- ===========

ausente2 :: Integral a => [a] -> a
ausente2 s@(x1:x2:x3:_) 
  | x1 + x3 /= 2 * x2 = x1 + (x3 - x2)
  | otherwise         = head [a | (a,b) <- zip [x1,x2..] s
                                , a /= b]

-- 3ª solución
-- ===========

ausente3 :: Integral a => [a] -> a
ausente3  xs@(x1:x2:_) 
  | null us   = x1 + v
  | otherwise = x2 + u * genericLength (u:us) 
  where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

-- Comprobación de equivalencia
-- ============================

-- Tipo de progresiones aritméticas con un término ausente.
newtype PAconAusente = PA [Integer]
  deriving Show

-- Generación de progresiones aritméticas con un término ausente. 
progresionConAusenteArbitraria :: Gen PAconAusente
progresionConAusenteArbitraria = do
  x <- arbitrary
  d <- arbitrary `suchThat` (/= 0)
  n <- arbitrary `suchThat` (> 2)
  k <- arbitrary `suchThat` (> 0)
  let (as,_:bs) = splitAt k [x,x+d..]
  frequency [(1,return (PA (as ++ bs))),
             (1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.
instance Arbitrary PAconAusente where
  arbitrary = progresionConAusenteArbitraria

-- La propiedad es
prop_ausente :: PAconAusente -> Bool 
prop_ausente (PA xs) =
  all (== ausente1 xs)
      [ausente2 xs,
       ausente3 xs]

-- La comprobación es
--    λ> quickCheck prop_ausente
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])
--    1000001
--    (1.15 secs, 560,529,520 bytes)
--    λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])
--    1000001
--    (0.33 secs, 336,530,680 bytes)
--    λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])
--    1000001
--    (0.50 secs, 498,047,584 bytes)

import Data.List (group, genericLength)

import Test.QuickCheck (Arbitrary, Gen,

arbitrary, frequency, suchThat, quickCheck)

-- 1ª solución

-- ===========

ausente1 :: Integral a => [a] -> a

ausente1 (x1:xs@(x2:x3:_))

| d1 == d2 = ausente1 xs

| d1 == 2 * d2 = x1 + d2

| d2 == 2 * d1 = x2 + d1

where d1 = x2 - x1

d2 = x3 - x2

-- 2ª solución

-- ===========

ausente2 :: Integral a => [a] -> a

ausente2 s@(x1:x2:x3:_)

| x1 + x3 /= 2 * x2 = x1 + (x3 - x2)

| otherwise = head [a | (a,b) <- zip [x1,x2..] s

, a /= b]

-- 3ª solución

-- ===========

ausente3 :: Integral a => [a] -> a

ausente3 xs@(x1:x2:_)

| null us = x1 + v

| otherwise = x2 + u * genericLength (u:us)

where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

-- Comprobación de equivalencia

-- ============================

-- Tipo de progresiones aritméticas con un término ausente.

newtype PAconAusente = PA [Integer]

deriving Show

-- Generación de progresiones aritméticas con un término ausente.

progresionConAusenteArbitraria :: Gen PAconAusente

progresionConAusenteArbitraria = do

x <- arbitrary

d <- arbitrary `suchThat` (/= 0)

n <- arbitrary `suchThat` (> 2)

k <- arbitrary `suchThat` (> 0)

let (as,_:bs) = splitAt k [x,x+d..]

frequency [(1,return (PA (as ++ bs))),

(1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.

instance Arbitrary PAconAusente where

arbitrary = progresionConAusenteArbitraria

-- La propiedad es

prop_ausente :: PAconAusente -> Bool

prop_ausente (PA xs) =

all (== ausente1 xs)

[ausente2 xs,

ausente3 xs]

-- La comprobación es

-- λ> quickCheck prop_ausente

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])

-- 1000001

-- (1.15 secs, 560,529,520 bytes)

-- λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])

-- 1000001

-- (0.33 secs, 336,530,680 bytes)

-- λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])

-- 1000001

-- (0.50 secs, 498,047,584 bytes)

El código se encuentra en GitHub.

Polinomios de Bell

PorJosé A. Alonso 31-mayo-202229-mayo-2022

Los polinomios de Bell forman una sucesión de polinomios, definida como sigue:

B₀(x) = 1 (polinomio unidad)
Bₙ(x) = x·[Bₙ(x) + Bₙ'(x)] (donde Bₙ'(x) es la derivada de Bₙ(x))

Por ejemplo,

   B₀(x) = 1                     = 1
   B₁(x) = x·(1+0)               = x     
   B₂(x) = x·(x+1)               = x²+x         
   B₃(x) = x·(x²+x+2x+1)         = x³+3x²+x    
   B₄(x) = x·(x³+3x²+x+3x²+6x+1) = x⁴+6x³+7x²+x

B₀(x) = 1 = 1

B₁(x) = x·(1+0) = x

B₂(x) = x·(x+1) = x²+x

B₃(x) = x·(x²+x+2x+1) = x³+3x²+x

B₄(x) = x·(x³+3x²+x+3x²+6x+1) = x⁴+6x³+7x²+x

Definir la función

   polBell :: Integer -> Polinomio Integer

1	polBell :: Integer -> Polinomio Integer

tal que (polBell n) es el polinomio de Bell de grado n. Por ejemplo,

   polBell 4                    ==  x^4 + 6*x^3 + 7*x^2 + 1*x
   coeficiente 2 (polBell 4)    ==  7
   coeficiente 2 (polBell 30)   ==  536870911
   coeficiente 1 (polBell 1000) == 1
   length (show (coeficiente 9 (polBell 2000)))  ==  1903

polBell 4 == x^4 + 6*x^3 + 7*x^2 + 1*x

coeficiente 2 (polBell 4) == 7

coeficiente 2 (polBell 30) == 536870911

coeficiente 1 (polBell 1000) == 1

length (show (coeficiente 9 (polBell 2000))) == 1903

Notas: Se usa la librería I1M.PolOperaciones que se encuentra aquí y se describe aquí. Además, en el último ejemplo se usa la función coeficiente tal que (coeficiente k p) es el coeficiente del término de grado k en el polinomio p definida por

   coeficiente :: Num a => Int -> Polinomio a -> a
   coeficiente k p | k == n                 = coefLider p
                   | k > grado (restoPol p) = 0
                   | otherwise              = coeficiente k (restoPol p)
                   where n = grado p

coeficiente :: Num a => Int -> Polinomio a -> a

coeficiente k p | k == n = coefLider p

| k > grado (restoPol p) = 0

| otherwise = coeficiente k (restoPol p)

where n = grado p

Soluciones

import Data.List          (genericIndex)
import I1M.PolOperaciones (Polinomio, coefLider, consPol, derivada,
                           grado, multPol, polCero, polUnidad, restoPol,
                           sumaPol) 
import Test.QuickCheck    (Positive (Positive), quickCheck)

-- Función auxiliar
-- ================

-- (coeficiente k p) es el coeficiente del término de grado k en el
-- polinomio p.
coeficiente :: Num a => Int -> Polinomio a -> a
coeficiente k p | k == n                 = coefLider p
                | k > grado (restoPol p) = 0
                | otherwise              = coeficiente k (restoPol p)
                where n = grado p

-- 1ª solución
-- ===========

polBell1 :: Integer -> Polinomio Integer
polBell1 0 = polUnidad
polBell1 n = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))
  where p = polBell1 (n-1)

-- 2ª solución
-- ===========

polBell2 :: Integer -> Polinomio Integer
polBell2 n = sucPolinomiosBell `genericIndex` n

sucPolinomiosBell :: [Polinomio Integer]
sucPolinomiosBell = iterate f polUnidad
  where f p = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_polBell :: Positive Integer -> Bool 
prop_polBell (Positive n) =
  polBell1 n == polBell2 n

-- La comprobación es
--    λ> quickCheck prop_polBell
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (show (coeficiente 9 (polBell1 2000)))
--    1903
--    (5.37 secs, 4,829,322,368 bytes)
--    λ> length (show (coeficiente 9 (polBell2 2000)))
--    1903
--    (4.03 secs, 4,825,094,064 bytes)

import Data.List (genericIndex)

import I1M.PolOperaciones (Polinomio, coefLider, consPol, derivada,

grado, multPol, polCero, polUnidad, restoPol,

sumaPol)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- Función auxiliar

-- ================

-- (coeficiente k p) es el coeficiente del término de grado k en el

-- polinomio p.

coeficiente :: Num a => Int -> Polinomio a -> a

coeficiente k p | k == n = coefLider p

| k > grado (restoPol p) = 0

| otherwise = coeficiente k (restoPol p)

where n = grado p

-- 1ª solución

-- ===========

polBell1 :: Integer -> Polinomio Integer

polBell1 0 = polUnidad

polBell1 n = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

where p = polBell1 (n-1)

-- 2ª solución

-- ===========

polBell2 :: Integer -> Polinomio Integer

polBell2 n = sucPolinomiosBell `genericIndex` n

sucPolinomiosBell :: [Polinomio Integer]

sucPolinomiosBell = iterate f polUnidad

where f p = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_polBell :: Positive Integer -> Bool

prop_polBell (Positive n) =

polBell1 n == polBell2 n

-- La comprobación es

-- λ> quickCheck prop_polBell

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (show (coeficiente 9 (polBell1 2000)))

-- 1903

-- (5.37 secs, 4,829,322,368 bytes)

-- λ> length (show (coeficiente 9 (polBell2 2000)))

-- 1903

-- (4.03 secs, 4,825,094,064 bytes)

El código se encuentra en GitHub.

Ordenación de los racionales

PorJosé A. Alonso 30-mayo-202229-mayo-2022

En este ejercicio, representamos las fracciones mediante pares de números de enteros.

Definir la función

   fraccionesOrd :: Integer -> [(Integer,Integer)]

1	fraccionesOrd :: Integer -> [(Integer,Integer)]

tal que (fraccionesOrd n) es la lista con las fracciones propias positivas ordenadas, con denominador menor o igual que n. Por ejemplo,

   λ> fraccionesOrd 4
   [(1,4),(1,3),(1,2),(2,3),(3,4)]
   λ> fraccionesOrd 5
   [(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

λ> fraccionesOrd 4

[(1,4),(1,3),(1,2),(2,3),(3,4)]

λ> fraccionesOrd 5

[(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

Soluciones

import Data.List (sort, sortBy)
import Data.Ratio ((%), numerator, denominator)
import Test.QuickCheck

-- 1ª solución
-- ===========

fraccionesOrd1 :: Integer -> [(Integer,Integer)]
fraccionesOrd1 n = 
  [(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))
                              | y <- [2..n], 
                                x <- [1..y-1], 
                                gcd x y == 1]]

-- 2ª solución
-- ===========

fraccionesOrd2 :: Integer -> [(Integer,Integer)]
fraccionesOrd2 n = 
  map snd (sort [(fromIntegral x/fromIntegral y,(x,y))
                 | y <- [2..n], 
                   x <- [1..y-1], 
                   gcd x y == 1])

-- 3ª solución
-- ===========

fraccionesOrd3 :: Integer -> [(Integer,Integer)]
fraccionesOrd3 n = 
  sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]
  where comp (a,b) (c,d) = compare (a*d) (b*c)

-- 4ª solución
-- ===========

fraccionesOrd4 :: Integer -> [(Integer,Integer)]
fraccionesOrd4 n = 
  [(numerator x, denominator x) | x <- racionalesOrd4 n]
  
-- (racionalesOrd4 n) es la lista con los racionales ordenados, con
-- denominador menor o igual que n. Por ejemplo,
--    λ> racionalesOrd4 4
--    [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]
racionalesOrd4 :: Integer -> [Rational]
racionalesOrd4 n =
  sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_fraccionesOrd :: Positive Integer -> Bool 
prop_fraccionesOrd (Positive n) =
  all (== fraccionesOrd1 n)
      [fraccionesOrd2 n,
       fraccionesOrd3 n,
       fraccionesOrd4 n]

-- La comprobación es
--    λ> quickCheck prop_fraccionesOrd
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (fraccionesOrd1 2000)
--    1216587
--    (3.65 secs, 2,879,842,368 bytes)
--    λ> length (fraccionesOrd2 2000)
--    1216587
--    (3.36 secs, 2,870,109,640 bytes)
--    λ> length (fraccionesOrd3 2000)
--    1216587
--    (8.83 secs, 5,700,519,584 bytes)
--    λ> length (fraccionesOrd4 2000)
--    1216587
--    (4.12 secs, 5,181,904,336 bytes)

import Data.List (sort, sortBy)

import Data.Ratio ((%), numerator, denominator)

import Test.QuickCheck

-- 1ª solución

-- ===========

fraccionesOrd1 :: Integer -> [(Integer,Integer)]

fraccionesOrd1 n =

[(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

x <- [1..y-1],

gcd x y == 1]]

-- 2ª solución

-- ===========

fraccionesOrd2 :: Integer -> [(Integer,Integer)]

fraccionesOrd2 n =

map snd (sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

x <- [1..y-1],

gcd x y == 1])

-- 3ª solución

-- ===========

fraccionesOrd3 :: Integer -> [(Integer,Integer)]

fraccionesOrd3 n =

sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]

where comp (a,b) (c,d) = compare (a*d) (b*c)

-- 4ª solución

-- ===========

fraccionesOrd4 :: Integer -> [(Integer,Integer)]

fraccionesOrd4 n =

[(numerator x, denominator x) | x <- racionalesOrd4 n]

-- (racionalesOrd4 n) es la lista con los racionales ordenados, con

-- denominador menor o igual que n. Por ejemplo,

-- λ> racionalesOrd4 4

-- [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]

racionalesOrd4 :: Integer -> [Rational]

racionalesOrd4 n =

sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_fraccionesOrd :: Positive Integer -> Bool

prop_fraccionesOrd (Positive n) =

all (== fraccionesOrd1 n)

[fraccionesOrd2 n,

fraccionesOrd3 n,

fraccionesOrd4 n]

-- La comprobación es

-- λ> quickCheck prop_fraccionesOrd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (fraccionesOrd1 2000)

-- 1216587

-- (3.65 secs, 2,879,842,368 bytes)

-- λ> length (fraccionesOrd2 2000)

-- 1216587

-- (3.36 secs, 2,870,109,640 bytes)

-- λ> length (fraccionesOrd3 2000)

-- 1216587

-- (8.83 secs, 5,700,519,584 bytes)

-- λ> length (fraccionesOrd4 2000)

-- 1216587

-- (4.12 secs, 5,181,904,336 bytes)

El código se encuentra en GitHub.

Polinomios cuadráticos generadores de primos

PorJosé A. Alonso 27-mayo-2022

En 1772, Euler publicó que el polinomio n² + n + 41 genera 40 números primos para todos los valores de n entre 0 y 39. Sin embargo, cuando n = 40, 40²+40+41 = 40(40+1)+41 es divisible por 41.

Usando ordenadores, se descubrió que el polinomio n² – 79n + 1601 genera 80 números primos para todos los valores de n entre 0 y 79.

Definir la función

   generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

1	generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

tal que (generadoresMaximales n) es el par (m,xs) donde

xs es la lista de pares (x,y) tales que n²+xn+y es uno de polinomios que genera un número máximo de números primos consecutivos a partir de cero entre todos los polinomios de la forma n²+an+b, con |a| ≤ n y |b| ≤ n y
m es dicho número máximo.

Por ejemplo,

   generadoresMaximales    4  ==  ( 3,[(-2,3),(-1,3),(3,3)])
   generadoresMaximales    6  ==  ( 5,[(-1,5),(5,5)])
   generadoresMaximales   41  ==  (41,[(-1,41)])
   generadoresMaximales   50  ==  (43,[(-5,47)])
   generadoresMaximales  100  ==  (48,[(-15,97)])
   generadoresMaximales  200  ==  (53,[(-25,197)])
   generadoresMaximales 1650  ==  (80,[(-79,1601)])

generadoresMaximales 4 == ( 3,[(-2,3),(-1,3),(3,3)])

generadoresMaximales 6 == ( 5,[(-1,5),(5,5)])

generadoresMaximales 41 == (41,[(-1,41)])

generadoresMaximales 50 == (43,[(-5,47)])

generadoresMaximales 100 == (48,[(-15,97)])

generadoresMaximales 200 == (53,[(-25,197)])

generadoresMaximales 1650 == (80,[(-79,1601)])

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (primes, isPrime)
import I1M.PolOperaciones (valor, consPol, polCero)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales1 n =
  (m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])
  where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos (-1) 41     ==  41
--    nPrimos (-79) 1601  ==  80
nPrimos :: Integer -> Integer -> Int
nPrimos a b =
  length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

-- 2ª solución
-- ===========

-- Notas:
-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que
--    0²+a*0+b = b es primo.
-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para
--    n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales2 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    a <- [-n..n],
                                    isPrime(1+a+b)]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- 3ª solución
-- ===========

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales3 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    p <- takeWhile (<=1+2*n) primes,
                                    let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)
-- ===========================================

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales4 n = (m,map snd zs)
  where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,
                                     p <- takeWhile (<=1+2*n) primes,
                                     let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos2 (-1) 41     ==  41
--    nPrimos2 (-79) 1601  ==  80
nPrimos2 :: Integer -> Integer -> Int
nPrimos2 a b =
  length (takeWhile isPrime [valor p n | n <- [0..]])
  where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_generadoresMaximales :: Positive Integer -> Bool
prop_generadoresMaximales (Positive n) =
  all (equivalentes (generadoresMaximales1 n'))
      [generadoresMaximales2 n',
       generadoresMaximales3 n',
       generadoresMaximales4 n']
  where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool
equivalentes (n,xs) (m,ys) =
  n == m && sort xs == sort ys

-- La comprobación es
--    λ> quickCheck prop_generadoresMaximales
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> generadoresMaximales1 300
--    (56,[(-31,281)])
--    (2.10 secs, 2,744,382,760 bytes)
--    λ> generadoresMaximales2 300
--    (56,[(-31,281)])
--    (0.17 secs, 382,103,656 bytes)
--    λ> generadoresMaximales3 300
--    (56,[(-31,281)])
--    (0.19 secs, 346,725,872 bytes)
--    λ> generadoresMaximales4 300
--    (56,[(-31,281)])
--    (0.20 secs, 388,509,808 bytes)

100

101

102

103

104

105

106

107

108

import Data.List (sort)

import Data.Numbers.Primes (primes, isPrime)

import I1M.PolOperaciones (valor, consPol, polCero)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales1 n =

(m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])

where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos (-1) 41 == 41

-- nPrimos (-79) 1601 == 80

nPrimos :: Integer -> Integer -> Int

nPrimos a b =

length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

-- 2ª solución

-- ===========

-- Notas:

-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que

-- 0²+a*0+b = b es primo.

-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para

-- n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales2 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

a <- [-n..n],

isPrime(1+a+b)]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 3ª solución

-- ===========

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales3 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)

-- ===========================================

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales4 n = (m,map snd zs)

where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos2 (-1) 41 == 41

-- nPrimos2 (-79) 1601 == 80

nPrimos2 :: Integer -> Integer -> Int

nPrimos2 a b =

length (takeWhile isPrime [valor p n | n <- [0..]])

where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_generadoresMaximales :: Positive Integer -> Bool

prop_generadoresMaximales (Positive n) =

all (equivalentes (generadoresMaximales1 n'))

[generadoresMaximales2 n',

generadoresMaximales3 n',

generadoresMaximales4 n']

where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool

equivalentes (n,xs) (m,ys) =

n == m && sort xs == sort ys

-- La comprobación es

-- λ> quickCheck prop_generadoresMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> generadoresMaximales1 300

-- (56,[(-31,281)])

-- (2.10 secs, 2,744,382,760 bytes)

-- λ> generadoresMaximales2 300

-- (56,[(-31,281)])

-- (0.17 secs, 382,103,656 bytes)

-- λ> generadoresMaximales3 300

-- (56,[(-31,281)])

-- (0.19 secs, 346,725,872 bytes)

-- λ> generadoresMaximales4 300

-- (56,[(-31,281)])

-- (0.20 secs, 388,509,808 bytes)

El código se encuentra en GitHub.

El triángulo de Floyd

PorJosé A. Alonso 26-mayo-2022

El triángulo de Floyd, llamado así en honor a Robert Floyd, es un triángulo rectángulo formado con números naturales. Para crear un triángulo de Floyd, se comienza con un 1 en la esquina superior izquierda, y se continúa escribiendo la secuencia de los números naturales de manera que cada línea contenga un número más que la anterior. Las 5 primeras líneas del triángulo de Floyd son

    1
    2   3
    4   5   6
    7   8   9  10
   11  12  13  14  15

2 3

4 5 6

7 8 9 10

11 12 13 14 15

Definir la función

   trianguloFloyd :: [[Integer]]

1	trianguloFloyd :: [[Integer]]

tal que trianguloFloyd es el triángulo de Floyd. Por ejemplo,

   λ> take 4 trianguloFloyd
   [[1],
    [2,3],
    [4,5,6],
    [7,8,9,10]]
  (trianguloFloyd !! (10^5)) !! 0  ==  5000050001
  (trianguloFloyd !! (10^6)) !! 0  ==  500000500001
  (trianguloFloyd !! (10^7)) !! 0  ==  50000005000001

λ> take 4 trianguloFloyd

[[1],

[2,3],

[4,5,6],

[7,8,9,10]]

(trianguloFloyd !! (10^5)) !! 0 == 5000050001

(trianguloFloyd !! (10^6)) !! 0 == 500000500001

(trianguloFloyd !! (10^7)) !! 0 == 50000005000001

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

trianguloFloyd1 :: [[Integer]]
trianguloFloyd1 = floyd 1 [1..]
  where floyd n xs = i : floyd (n+1) r
          where (i,r) = splitAt n xs

-- 2ª solución
-- ===========

trianguloFloyd2 :: [[Integer]]
trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el
-- triángulo de Floyd. Por ejemplo,
--    siguienteF [2,3]    ==  [4,5,6]
--    siguienteF [4,5,6]  ==  [7,8,9,10]
siguienteF :: [Integer] -> [Integer]
siguienteF xs = [a..a+n]
    where a = 1 + last xs
          n = genericLength xs

-- 3ª solución
-- ===========

trianguloFloyd3 :: [[Integer]]
trianguloFloyd3 =
  [[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución
-- ===========

trianguloFloyd4 :: [[Integer]]
trianguloFloyd4 =
  scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_trianguloFloyd :: Positive Int -> Bool
prop_trianguloFloyd (Positive n) =
  all (== (trianguloFloyd1 !! n))
      [trianguloFloyd2 !! n,
       trianguloFloyd3 !! n,
       trianguloFloyd4 !! n]

-- La comprobación es
-- λ> quickCheck prop_trianguloFloyd
-- +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> (trianguloFloyd1 !! 5000) !! 5000
--    12507501
--    (1.47 secs, 2,505,005,752 bytes)
--    λ> (trianguloFloyd2 !! 5000) !! 5000
--    12507501
--    (0.79 secs, 2,416,259,176 bytes)
--    λ> (trianguloFloyd3 !! 5000) !! 5000
--    12507501
--    (0.00 secs, 1,809,152 bytes)
--    λ> (trianguloFloyd4 !! 5000) !! 5000
--    12507501
--    (0.01 secs, 3,517,896 bytes)
--
--    λ> (trianguloFloyd3 !! (10^7)) !! 0
--    50000005000001
--    (2.45 secs, 1,656,534,080 bytes)
--    λ> (trianguloFloyd4 !! (10^7)) !! 0
--    50000005000001
--    (10.86 secs, 5,302,760,752 bytes)

import Data.List (genericLength)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

trianguloFloyd1 :: [[Integer]]

trianguloFloyd1 = floyd 1 [1..]

where floyd n xs = i : floyd (n+1) r

where (i,r) = splitAt n xs

-- 2ª solución

-- ===========

trianguloFloyd2 :: [[Integer]]

trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el

-- triángulo de Floyd. Por ejemplo,

-- siguienteF [2,3] == [4,5,6]

-- siguienteF [4,5,6] == [7,8,9,10]

siguienteF :: [Integer] -> [Integer]

siguienteF xs = [a..a+n]

where a = 1 + last xs

n = genericLength xs

-- 3ª solución

-- ===========

trianguloFloyd3 :: [[Integer]]

trianguloFloyd3 =

[[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución

-- ===========

trianguloFloyd4 :: [[Integer]]

trianguloFloyd4 =

scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_trianguloFloyd :: Positive Int -> Bool

prop_trianguloFloyd (Positive n) =

all (== (trianguloFloyd1 !! n))

[trianguloFloyd2 !! n,

trianguloFloyd3 !! n,

trianguloFloyd4 !! n]

-- La comprobación es

-- λ> quickCheck prop_trianguloFloyd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> (trianguloFloyd1 !! 5000) !! 5000

-- 12507501

-- (1.47 secs, 2,505,005,752 bytes)

-- λ> (trianguloFloyd2 !! 5000) !! 5000

-- 12507501

-- (0.79 secs, 2,416,259,176 bytes)

-- λ> (trianguloFloyd3 !! 5000) !! 5000

-- 12507501

-- (0.00 secs, 1,809,152 bytes)

-- λ> (trianguloFloyd4 !! 5000) !! 5000

-- 12507501

-- (0.01 secs, 3,517,896 bytes)

-- λ> (trianguloFloyd3 !! (10^7)) !! 0

-- 50000005000001

-- (2.45 secs, 1,656,534,080 bytes)

-- λ> (trianguloFloyd4 !! (10^7)) !! 0

-- 50000005000001

-- (10.86 secs, 5,302,760,752 bytes)

El código se encuentra en GitHub.

Numeración con base múltiple

PorJosé A. Alonso 25-mayo-202225-mayo-2022

Sea (b(i) | i ≥ 1) una sucesión infinita de números enteros mayores que 1. Entonces todo entero x mayor que cero se puede escribir de forma única como

   x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

1	x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

donde cada x(i) satisface la condición 0 ≤ x(i) < b(i+1). Se dice que [x(n),x(n-1),…,x(2),x(1),x(0)] es la representación de x en la base (b(i)). Por ejemplo, la representación de 377 en la base (2, 6, 8, …) es [7,5,0,1] ya que

   377 = 1 + 0*2 + 5*2*4 + 7*2*4*6

1	377 = 1 + 02 + 524 + 7246

y, además, 0 ≤ 1 < 2, 0 ≤ 0 < 4, 0 ≤ 5 < 6 y 0 ≤ 7 < 8.

Definir las funciones

   decimalAmultiple :: [Integer] -> Integer -> [Integer]
   multipleAdecimal :: [Integer] -> [Integer] -> Integer

1 2	decimalAmultiple :: [Integer] -> Integer -> [Integer] multipleAdecimal :: [Integer] -> [Integer] -> Integer

tales que

(decimalAmultiple bs x) es la representación del número x en la base bs. Por ejemplo,

     decimalAmultiple [2,4..] 377                      ==  [7,5,0,1]
     decimalAmultiple [2,5..] 377                      ==  [4,5,3,1]
     decimalAmultiple [2^n | n <- [1..]] 2015          ==  [1,15,3,3,1]
     decimalAmultiple (repeat 10) 2015                 ==  [2,0,1,5]

decimalAmultiple [2,4..] 377 == [7,5,0,1]

decimalAmultiple [2,5..] 377 == [4,5,3,1]

decimalAmultiple [2^n | n <- [1..]] 2015 == [1,15,3,3,1]

decimalAmultiple (repeat 10) 2015 == [2,0,1,5]

(multipleAdecimal bs cs) es el número decimal cuya representación en la base bs es cs. Por ejemplo,

     multipleAdecimal [2,4..] [7,5,0,1]                ==  377
     multipleAdecimal [2,5..] [4,5,3,1]                ==  377
     multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1]  ==  2015
     multipleAdecimal (repeat 10) [2,0,1,5]            ==  2015

multipleAdecimal [2,4..] [7,5,0,1] == 377

multipleAdecimal [2,5..] [4,5,3,1] == 377

multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1] == 2015

multipleAdecimal (repeat 10) [2,0,1,5] == 2015

Comprobar con QuickCheck que se verifican las siguientes propiedades

Para cualquier base bs y cualquier entero positivo n,

     multipleAdecimal bs (decimalAmultiple bs x) == x

1	multipleAdecimal bs (decimalAmultiple bs x) == x

Para cualquier base bs y cualquier entero positivo n, el coefiente i-ésimo de la representación múltiple de n en la base bs es un entero no negativo menos que el i-ésimo elemento de bs.

Soluciones

import Test.QuickCheck
import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple
-- ===============================

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]
decimalAmultiple1 bs n = reverse (aux bs n)
  where aux _ 0      = []
        aux (d:ds) m = r : aux ds q
          where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple
-- ===============================

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]
decimalAmultiple2 bs n = aux bs n []
  where aux _ 0  xs     = xs
        aux (d:ds) m xs = aux ds q (r:xs)
          where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple
-- ===============================

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]
decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))
  where f (_     ,0) = Nothing
        f ((y:ys),m) = Just (r,(ys,q))
                       where (q,r) = quotRem m y

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_decimalAmultiple (InfiniteList xs _) (Positive n) =
  all (== decimalAmultiple1 xs' n)
      [decimalAmultiple2 xs' n,
       decimalAmultiple3 xs' n]
  where xs' = map getPositive xs

-- Comparación de eficiencia de decimalAmultiple
-- =============================================

-- La comparación es
--    λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))
--    21731
--    (0.45 secs, 486,085,256 bytes)
--    λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))
--    21731
--    (0.32 secs, 485,563,664 bytes)
--    λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))
--    21731
--    (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal
-- ===============================

multipleAdecimal1  :: [Integer] -> [Integer] -> Integer
multipleAdecimal1 xs ns = aux xs (reverse ns)
  where aux (y:ys) (m:ms) = m + y * (aux ys ms)
        aux _ _           = 0

-- 2ª solución de multipleAdecimal
-- ===============================

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer
multipleAdecimal2 bs xs =
  sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

-- Comprobación de equivalencia de multipleAdecimal
-- ================================================

-- La propiedad es
prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool
prop_multipleAdecimal (InfiniteList xs _) ys =
  multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'
  where xs' = map getPositive xs
        ys' = map getPositive ys

-- La comprobación es
--    λ> quickCheck prop_multipleAdecimal
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia de multipleAdecimal
-- =============================================

-- La comparación es
--    λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))
--    35660
--    (0.14 secs, 179,522,152 bytes)
--    λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))
--    35660
--    (0.22 secs, 243,368,664 bytes)

-- Comprobación de las propiedades
-- ===============================

-- La primera propiedad es
prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_inversas (InfiniteList xs _) (Positive n) =
  multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n
  where xs' = map getPositive xs

-- Su comprobación es
--    λ> quickCheck prop_inversas
--    +++ OK, passed 100 tests.

-- la 2ª propiedad es
prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_coeficientes (InfiniteList xs _) (Positive n) =
  and [0 <= c && c < b | (c,b) <- zip cs xs']
  where xs' = map getPositive xs
        cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es
--    λ> quickCheck prop_coeficientes
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

import Test.QuickCheck

import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple

-- ===============================

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]

decimalAmultiple1 bs n = reverse (aux bs n)

where aux _ 0 = []

aux (d:ds) m = r : aux ds q

where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple

-- ===============================

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]

decimalAmultiple2 bs n = aux bs n []

where aux _ 0 xs = xs

aux (d:ds) m xs = aux ds q (r:xs)

where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple

-- ===============================

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]

decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))

where f (_ ,0) = Nothing

f ((y:ys),m) = Just (r,(ys,q))

where (q,r) = quotRem m y

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_decimalAmultiple (InfiniteList xs _) (Positive n) =

all (== decimalAmultiple1 xs' n)

[decimalAmultiple2 xs' n,

decimalAmultiple3 xs' n]

where xs' = map getPositive xs

-- Comparación de eficiencia de decimalAmultiple

-- =============================================

-- La comparación es

-- λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))

-- 21731

-- (0.45 secs, 486,085,256 bytes)

-- λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))

-- 21731

-- (0.32 secs, 485,563,664 bytes)

-- λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))

-- 21731

-- (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal

-- ===============================

multipleAdecimal1 :: [Integer] -> [Integer] -> Integer

multipleAdecimal1 xs ns = aux xs (reverse ns)

where aux (y:ys) (m:ms) = m + y * (aux ys ms)

aux _ _ = 0

-- 2ª solución de multipleAdecimal

-- ===============================

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer

multipleAdecimal2 bs xs =

sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

-- Comprobación de equivalencia de multipleAdecimal

-- ================================================

-- La propiedad es

prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool

prop_multipleAdecimal (InfiniteList xs _) ys =

multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'

where xs' = map getPositive xs

ys' = map getPositive ys

-- La comprobación es

-- λ> quickCheck prop_multipleAdecimal

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de multipleAdecimal

-- =============================================

-- La comparación es

-- λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))

-- 35660

-- (0.14 secs, 179,522,152 bytes)

-- λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))

-- 35660

-- (0.22 secs, 243,368,664 bytes)

-- Comprobación de las propiedades

-- ===============================

-- La primera propiedad es

prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_inversas (InfiniteList xs _) (Positive n) =

multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n

where xs' = map getPositive xs

-- Su comprobación es

-- λ> quickCheck prop_inversas

-- +++ OK, passed 100 tests.

-- la 2ª propiedad es

prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_coeficientes (InfiniteList xs _) (Positive n) =

and [0 <= c && c < b | (c,b) <- zip cs xs']

where xs' = map getPositive xs

cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es

-- λ> quickCheck prop_coeficientes

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Matriz zigzagueante

PorJosé A. Alonso 24-mayo-202224-mayo-2022

La matriz zizagueante de orden n es la matriz cuadrada con n filas y n columnas y cuyos elementos son los n² primeros números naturales colocados de manera creciente a lo largo de las diagonales secundarias. Por ejemplo, La matriz zigzagueante de orden 5 es

    0  1  5  6 14
    2  4  7 13 15
    3  8 12 16 21
    9 11 17 20 22
   10 18 19 23 24

0 1 5 6 14

2 4 7 13 15

3 8 12 16 21

9 11 17 20 22

10 18 19 23 24

La colocación de los elementos se puede ver gráficamente en esta figura

Definir la función

   zigZag :: Int -> Matrix Int

1	zigZag :: Int -> Matrix Int

tal que (zigZag n) es la matriz zigzagueante de orden n. Por ejemplo,

   λ> zigZag 5
   ┌                ┐
   │  0  1  5  6 14 │
   │  2  4  7 13 15 │
   │  3  8 12 16 21 │
   │  9 11 17 20 22 │
   │ 10 18 19 23 24 │
   └                ┘
   λ> zigZag 4
   ┌             ┐
   │  0  1  5  6 │
   │  2  4  7 12 │
   │  3  8 11 13 │
   │  9 10 14 15 │
   └             ┘
   λ> maximum (zigZag 1500)
   2249999

λ> zigZag 5

┌ ┐

│ 0 1 5 6 14 │

│ 2 4 7 13 15 │

│ 3 8 12 16 21 │

│ 9 11 17 20 22 │

│ 10 18 19 23 24 │

└ ┘

λ> zigZag 4

┌ ┐

│ 0 1 5 6 │

│ 2 4 7 12 │

│ 3 8 11 13 │

│ 9 10 14 15 │

└ ┘

λ> maximum (zigZag 1500)

2249999

Soluciones

import Data.List (sort, sortBy)
import Data.Matrix (Matrix, fromList)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

zigZag1 :: Int -> Matrix Int
zigZag1 n = fromList n n (elementosZigZag n)

-- (elementosZigZag n) es la lista de los elementos de la matriz
-- zizagueante de orden n. Por ejemplo.
--    λ> elementosZigZag 5
--    [0,1,5,6,14,2,4,7,13,15,3,8,12,16,21,9,11,17,20,22,10,18,19,23,24]
elementosZigZag :: Int -> [Int]
elementosZigZag n =
  map snd (sort (zip (ordenZigZag n) [0..]))

-- (ordenZigZag n) es la lista de puntos del cuadrado nxn recorridos en
-- zig-zag por las diagonales secundarias. Por ejemplo,
--    λ> ordenZigZag 4
--    [(1,1), (1,2),(2,1), (3,1),(2,2),(1,3), (1,4),(2,3),(3,2),(4,1),
--     (4,2),(3,3),(2,4), (3,4),(4,3), (4,4)]
ordenZigZag :: Int -> [(Int,Int)]
ordenZigZag n = concat [aux n m | m <- [2..2*n]]
    where aux k m | odd m     = [(x,m-x) | x <- [max 1 (m-k)..min k (m-1)]]
                  | otherwise = [(m-x,x) | x <- [max 1 (m-k)..min k (m-1)]]

-- 2ª solución
-- ===========

zigZag2 :: Int -> Matrix Int
zigZag2 n = fromList n n (elementosZigZag2 n)

elementosZigZag2 :: Int -> [Int]
elementosZigZag2 n =
  map snd (sort (zip (ordenZigZag2 n) [0..]))

ordenZigZag2 :: Int -> [(Int,Int)]
ordenZigZag2 n = sortBy comp [(x,y) | x <- [1..n], y <- [1..n]]
    where comp (x1,y1) (x2,y2) | x1+y1 < x2+y2 = LT
                               | x1+y1 > x2+y2 = GT
                               | even (x1+y1)  = compare y1 y2
                               | otherwise     = compare x1 x2

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_zigZag :: Positive Int -> Bool
prop_zigZag (Positive n) =
  zigZag1 n == zigZag2 n

-- La comprobación es
--    λ> quickCheck prop_zigZag
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (zigZag1 5000)
--    25000000
--    (2.57 secs, 1,800,683,952 bytes)
--    λ> length (zigZag2 5000)
--    25000000
--    (2.20 secs, 1,800,683,952 bytes)
--
--    λ> maximum (zigZag1 1100)
--    1209999
--    (2.12 secs, 1,840,095,864 bytes)
--    λ> maximum (zigZag2 1100)
--    1209999
--    (21.27 secs, 11,661,088,256 bytes)

import Data.List (sort, sortBy)

import Data.Matrix (Matrix, fromList)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

zigZag1 :: Int -> Matrix Int

zigZag1 n = fromList n n (elementosZigZag n)

-- (elementosZigZag n) es la lista de los elementos de la matriz

-- zizagueante de orden n. Por ejemplo.

-- λ> elementosZigZag 5

-- [0,1,5,6,14,2,4,7,13,15,3,8,12,16,21,9,11,17,20,22,10,18,19,23,24]

elementosZigZag :: Int -> [Int]

elementosZigZag n =

map snd (sort (zip (ordenZigZag n) [0..]))

-- (ordenZigZag n) es la lista de puntos del cuadrado nxn recorridos en

-- zig-zag por las diagonales secundarias. Por ejemplo,

-- λ> ordenZigZag 4

-- [(1,1), (1,2),(2,1), (3,1),(2,2),(1,3), (1,4),(2,3),(3,2),(4,1),

-- (4,2),(3,3),(2,4), (3,4),(4,3), (4,4)]

ordenZigZag :: Int -> [(Int,Int)]

ordenZigZag n = concat [aux n m | m <- [2..2*n]]

where aux k m | odd m = [(x,m-x) | x <- [max 1 (m-k)..min k (m-1)]]

| otherwise = [(m-x,x) | x <- [max 1 (m-k)..min k (m-1)]]

-- 2ª solución

-- ===========

zigZag2 :: Int -> Matrix Int

zigZag2 n = fromList n n (elementosZigZag2 n)

elementosZigZag2 :: Int -> [Int]

elementosZigZag2 n =

map snd (sort (zip (ordenZigZag2 n) [0..]))

ordenZigZag2 :: Int -> [(Int,Int)]

ordenZigZag2 n = sortBy comp [(x,y) | x <- [1..n], y <- [1..n]]

where comp (x1,y1) (x2,y2) | x1+y1 < x2+y2 = LT

| x1+y1 > x2+y2 = GT

| even (x1+y1) = compare y1 y2

| otherwise = compare x1 x2

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_zigZag :: Positive Int -> Bool

prop_zigZag (Positive n) =

zigZag1 n == zigZag2 n

-- La comprobación es

-- λ> quickCheck prop_zigZag

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (zigZag1 5000)

-- 25000000

-- (2.57 secs, 1,800,683,952 bytes)

-- λ> length (zigZag2 5000)

-- 25000000

-- (2.20 secs, 1,800,683,952 bytes)

-- λ> maximum (zigZag1 1100)

-- 1209999

-- (2.12 secs, 1,840,095,864 bytes)

-- λ> maximum (zigZag2 1100)

-- 1209999

-- (21.27 secs, 11,661,088,256 bytes)

El código se encuentra en GitHub.

Densidades de números abundantes, perfectos y deficientes

PorJosé A. Alonso 23-mayo-202223-mayo-2022

La n-ésima densidad de un tipo de número es el cociente entre la cantidad de los números entre 1 y n que son del tipo considerado y n. Por ejemplo, la 7-ésima densidad de los múltiplos de 3 es 2/7 ya que entre los 7 primeros números sólo 2 son múltiplos de 3.

Definir las funciones

   densidades :: Int -> (Double,Double,Double)
   graficas   :: Int -> IO ()

1 2	densidades :: Int -> (Double,Double,Double) graficas :: Int -> IO ()

tales que

(densidades n) es la terna formada por la n-ésima densidad
- de los números abundantes (es decir, para los que la suma de sus divisores propios es mayor que el número),
- de los números perfectos (es decir, para los que la suma de sus divisores propios es mayor que el número) y
- de los números deficientes (es decir, para los que la suma de sus divisores propios es menor que el número).
Por ejemplo,

    densidades 100     ==  (0.22,    2.0e-2, 0.76)
    densidades 1000    ==  (0.246,   3.0e-3, 0.751)
    densidades 10000   ==  (0.2488,  4.0e-4, 0.7508)
    densidades 100000  ==  (0.24795, 4.0e-5, 0.75201)
    densidades 1000000 ==  (0.247545,4.0e-6, 0.752451)

densidades 100 == (0.22, 2.0e-2, 0.76)

densidades 1000 == (0.246, 3.0e-3, 0.751)

densidades 10000 == (0.2488, 4.0e-4, 0.7508)

densidades 100000 == (0.24795, 4.0e-5, 0.75201)

densidades 1000000 == (0.247545,4.0e-6, 0.752451)

(graficas n) dibuja las gráficas de las k-ésimas densidades (para k entre 1 y n) de los números abundantes, de los números perfectos y de los números deficientes. Por ejemplo, (graficas 100) dibuja

y (graficas 400) dibuja

Soluciones

import Data.List (genericLength, group, partition)
import Data.Array (accumArray, assocs)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple (plotLists, Attribute (Key))
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

densidades1 :: Int -> (Double,Double,Double)
densidades1 n = (f a, f p, f d)
  where a   = nAbundantes n
        p   = nPerfectos n
        d   = n - a - p
        f x = fromIntegral x / fromIntegral n

-- (nAbundantes n) es la cantidad de números abundantes desde 1 hasta
-- n. Por ejemplo,
--    nAbundantes 100 == 22
nAbundantes :: Int -> Int
nAbundantes n = length (filter esAbundante [1..n])

-- (esAbundante n) se verifica si n es un número abundante. Por ejemplo,
--    esAbundante 12 == True
--    esAbundante 22 == False
esAbundante :: Int -> Bool
esAbundante n = sumaDivisores n > n

-- (sumaDivisores n) es la suma de los divisores propios de n. Por
-- ejemplo,
--    sumaDivisores 12 == 16
--    sumaDivisores 22 == 14
sumaDivisores :: Int -> Int
sumaDivisores = sum . divisores

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,
--    divisores 12== [1,2,3,4,6]
--    divisores 22== [1,2,11]
divisores :: Int -> [Int]
divisores n = [x | x <- [1..n-1],
                   n `mod` x == 0]

-- (nPerfectos n) es la cantidad de números perfectos desde 1 hasta
-- n. Por ejemplo,
--    nPerfectos 100 == 2
nPerfectos :: Int -> Int
nPerfectos n = length (filter esPerfecto [1..n])

-- (esPerfecto n) se verifica si n es un número perfecto. Por ejemplo,
--    esPerfecto 28 == True
--    esPerfecto 38 == False
esPerfecto :: Int -> Bool
esPerfecto n = sumaDivisores n == n

-- 2ª solución
-- ===========

densidades2 :: Int -> (Double,Double,Double)
densidades2 n = (f as, f ps, f ds)
  where (as,pds) = partition esAbundante [1..n]
        (ps,ds)  = partition esPerfecto pds
        f xs     = genericLength xs / fromIntegral n

-- 3ª solución
-- ===========

densidades3 :: Int -> (Double,Double,Double)
densidades3 n = (f as, f ps, f ds)
  where cs       = map clasificacion [1..n]
        (as,pds) = partition (== Abundante) cs
        (ps,ds)  = partition (== Perfecto) pds
        f xs     = genericLength xs / fromIntegral n

data Clase = Abundante | Perfecto | Deficiente
  deriving (Eq, Show)

-- (clasificacion n) es la clase de número de n. Por ejemplo,
--    clasificacion 12 == Abundante
--    clasificacion 22 == Deficiente
--    clasificacion 28 == Perfecto
clasificacion :: Int -> Clase
clasificacion n
  | sd > n    = Abundante
  | sd < n    = Deficiente
  | otherwise = Perfecto
  where sd = sumaDivisores n

-- 4ª solución
-- ===========

densidades4 :: Int -> (Double,Double,Double)
densidades4 n = (f as, f ps, f ds)
  where cs       = map clasificacion2 [1..n]
        (as,pds) = partition (== Abundante) cs
        (ps,ds)  = partition (== Perfecto) pds
        f xs     = genericLength xs / fromIntegral n

-- 2ª definición de clasificacion
clasificacion2 :: Int -> Clase
clasificacion2 n
  | sd > n    = Abundante
  | sd < n    = Deficiente
  | otherwise = Perfecto
  where sd = sumaDivisores2 n

-- 2ª definición de sumaDivisores
sumaDivisores2 :: Int -> Int
sumaDivisores2 x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Int -> [(Int,Int)]
factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Int)
primeroYlongitud (x:xs) =
  (x, 1 + length xs)

-- 5ª solución
-- ===========

densidades5 :: Int -> (Double,Double,Double)
densidades5 n = (f a, f p, f d)
  where (a,p,d) = distribucion n
        f x = fromIntegral x / fromIntegral n

-- (distribucion n) es la terna (a,p,d) donde a es la cantidad de
-- números abundantes de 1 a n, p la de los perfectos y d la de los
-- deficientes. Por ejemplo,
--    distribucion 100  ==  (22,2,76)
distribucion :: Int -> (Int,Int,Int)
distribucion n = aux (0,0,0) (sumaDivisoresHasta n)
  where aux (a,p,d) [] = (a,p,d)
        aux (a,p,d) ((x,y):xys)
          | x < y     = aux (1+a,p,d) xys
          | x > y     = aux (a,p,1+d) xys
          | otherwise = aux (a,1+p,d) xys

-- (sumaDivisoresHasta n) es la lista de los pares (a,b) tales que a
-- varía entre 1 y n y b es la suma de los divisores propios de a. Por
-- ejemplo,
--    λ> sumaDivisoresHasta 12
--    [(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]
sumaDivisoresHasta :: Int -> [(Int,Int)]
sumaDivisoresHasta n =
  assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a está
-- entre 2 y n y b es un divisor propio e x. Por ejemplo,
--    λ> divisoresHasta 6
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]
--    λ> divisoresHasta 8
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]
divisoresHasta :: Int -> [(Int,Int)]
divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_densidades :: Positive Int -> Bool
prop_densidades (Positive n) =
  all (== densidades1 n)
      [ densidades2 n
      , densidades3 n
      , densidades4 n
      , densidades5 n
      ]

-- La comprobación es
--    λ> quickCheck prop_densidades
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> densidades1 2000
--    (0.2465,1.5e-3,0.752)
--    (1.59 secs, 804,883,768 bytes)
--    λ> densidades2 2000
--    (0.2465,1.5e-3,0.752)
--    (1.39 secs, 704,749,552 bytes)
--    λ> densidades3 2000
--    (0.2465,1.5e-3,0.752)
--    (0.81 secs, 403,783,312 bytes)
--    λ> densidades4 2000
--    (0.2465,1.5e-3,0.752)
--    (0.04 secs, 34,364,960 bytes)
--    λ> densidades5 2000
--    (0.2465,1.5e-3,0.752)
--    (0.02 secs, 4,716,720 bytes)
--
--    λ> densidades4 100000
--    (0.24795,4.0e-5,0.75201)
--    (1.61 secs, 4,826,971,624 bytes)
--    λ> densidades5 100000
--    (0.24795,4.0e-5,0.75201)
--    (0.43 secs, 291,989,272 bytes)

-- Gráfica
-- =======

graficas :: Int -> IO ()
graficas n =
  plotLists [Key Nothing]
            [ [x | (x,_,_) <- ts]
            , [y | (_,y,_) <- ts]
            , [z | (_,_,z) <- ts]]
  where ts = [densidades5 k | k <- [1..n]]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

import Data.List (genericLength, group, partition)

import Data.Array (accumArray, assocs)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple (plotLists, Attribute (Key))

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

densidades1 :: Int -> (Double,Double,Double)

densidades1 n = (f a, f p, f d)

where a = nAbundantes n

p = nPerfectos n

d = n - a - p

f x = fromIntegral x / fromIntegral n

-- (nAbundantes n) es la cantidad de números abundantes desde 1 hasta

-- n. Por ejemplo,

-- nAbundantes 100 == 22

nAbundantes :: Int -> Int

nAbundantes n = length (filter esAbundante [1..n])

-- (esAbundante n) se verifica si n es un número abundante. Por ejemplo,

-- esAbundante 12 == True

-- esAbundante 22 == False

esAbundante :: Int -> Bool

esAbundante n = sumaDivisores n > n

-- (sumaDivisores n) es la suma de los divisores propios de n. Por

-- ejemplo,

-- sumaDivisores 12 == 16

-- sumaDivisores 22 == 14

sumaDivisores :: Int -> Int

sumaDivisores = sum . divisores

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,

-- divisores 12== [1,2,3,4,6]

-- divisores 22== [1,2,11]

divisores :: Int -> [Int]

divisores n = [x | x <- [1..n-1],

n `mod` x == 0]

-- (nPerfectos n) es la cantidad de números perfectos desde 1 hasta

-- n. Por ejemplo,

-- nPerfectos 100 == 2

nPerfectos :: Int -> Int

nPerfectos n = length (filter esPerfecto [1..n])

-- (esPerfecto n) se verifica si n es un número perfecto. Por ejemplo,

-- esPerfecto 28 == True

-- esPerfecto 38 == False

esPerfecto :: Int -> Bool

esPerfecto n = sumaDivisores n == n

-- 2ª solución

-- ===========

densidades2 :: Int -> (Double,Double,Double)

densidades2 n = (f as, f ps, f ds)

where (as,pds) = partition esAbundante [1..n]

(ps,ds) = partition esPerfecto pds

f xs = genericLength xs / fromIntegral n

-- 3ª solución

-- ===========

densidades3 :: Int -> (Double,Double,Double)

densidades3 n = (f as, f ps, f ds)

where cs = map clasificacion [1..n]

(as,pds) = partition (== Abundante) cs

(ps,ds) = partition (== Perfecto) pds

f xs = genericLength xs / fromIntegral n

data Clase = Abundante | Perfecto | Deficiente

deriving (Eq, Show)

-- (clasificacion n) es la clase de número de n. Por ejemplo,

-- clasificacion 12 == Abundante

-- clasificacion 22 == Deficiente

-- clasificacion 28 == Perfecto

clasificacion :: Int -> Clase

clasificacion n

| sd > n = Abundante

| sd < n = Deficiente

| otherwise = Perfecto

where sd = sumaDivisores n

-- 4ª solución

-- ===========

densidades4 :: Int -> (Double,Double,Double)

densidades4 n = (f as, f ps, f ds)

where cs = map clasificacion2 [1..n]

(as,pds) = partition (== Abundante) cs

(ps,ds) = partition (== Perfecto) pds

f xs = genericLength xs / fromIntegral n

-- 2ª definición de clasificacion

clasificacion2 :: Int -> Clase

clasificacion2 n

| sd > n = Abundante

| sd < n = Deficiente

| otherwise = Perfecto

where sd = sumaDivisores2 n

-- 2ª definición de sumaDivisores

sumaDivisores2 :: Int -> Int

sumaDivisores2 x =

product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la

-- descomposición prima de x. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

factorizacion :: Int -> [(Int,Int)]

factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs

-- y la longitud de xs. Por ejemplo,

-- primeroYlongitud [3,2,5,7] == (3,4)

primeroYlongitud :: [a] -> (a,Int)

primeroYlongitud (x:xs) =

(x, 1 + length xs)

-- 5ª solución

-- ===========

densidades5 :: Int -> (Double,Double,Double)

densidades5 n = (f a, f p, f d)

where (a,p,d) = distribucion n

f x = fromIntegral x / fromIntegral n

-- (distribucion n) es la terna (a,p,d) donde a es la cantidad de

-- números abundantes de 1 a n, p la de los perfectos y d la de los

-- deficientes. Por ejemplo,

-- distribucion 100 == (22,2,76)

distribucion :: Int -> (Int,Int,Int)

distribucion n = aux (0,0,0) (sumaDivisoresHasta n)

where aux (a,p,d) [] = (a,p,d)

aux (a,p,d) ((x,y):xys)

| x < y = aux (1+a,p,d) xys

| x > y = aux (a,p,1+d) xys

| otherwise = aux (a,1+p,d) xys

-- (sumaDivisoresHasta n) es la lista de los pares (a,b) tales que a

-- varía entre 1 y n y b es la suma de los divisores propios de a. Por

-- ejemplo,

-- λ> sumaDivisoresHasta 12

-- [(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]

sumaDivisoresHasta :: Int -> [(Int,Int)]

sumaDivisoresHasta n =

assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a está

-- entre 2 y n y b es un divisor propio e x. Por ejemplo,

-- λ> divisoresHasta 6

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]

-- λ> divisoresHasta 8

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]

divisoresHasta :: Int -> [(Int,Int)]

divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_densidades :: Positive Int -> Bool

prop_densidades (Positive n) =

all (== densidades1 n)

[ densidades2 n

, densidades3 n

, densidades4 n

, densidades5 n

]

-- La comprobación es

-- λ> quickCheck prop_densidades

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> densidades1 2000

-- (0.2465,1.5e-3,0.752)

-- (1.59 secs, 804,883,768 bytes)

-- λ> densidades2 2000

-- (0.2465,1.5e-3,0.752)

-- (1.39 secs, 704,749,552 bytes)

-- λ> densidades3 2000

-- (0.2465,1.5e-3,0.752)

-- (0.81 secs, 403,783,312 bytes)

-- λ> densidades4 2000

-- (0.2465,1.5e-3,0.752)

-- (0.04 secs, 34,364,960 bytes)

-- λ> densidades5 2000

-- (0.2465,1.5e-3,0.752)

-- (0.02 secs, 4,716,720 bytes)

-- λ> densidades4 100000

-- (0.24795,4.0e-5,0.75201)

-- (1.61 secs, 4,826,971,624 bytes)

-- λ> densidades5 100000

-- (0.24795,4.0e-5,0.75201)

-- (0.43 secs, 291,989,272 bytes)

-- Gráfica

-- =======

graficas :: Int -> IO ()

graficas n =

plotLists [Key Nothing]

[ [x | (x,_,_) <- ts]

, [y | (_,y,_) <- ts]

, [z | (_,_,z) <- ts]]

where ts = [densidades5 k | k <- [1..n]]

El código se encuentra en GitHub.

Sumas de divisores propios

PorJosé A. Alonso 20-mayo-2022

Definir la función

   sumaDivisoresHasta :: Integer -> [(Integer,Integer)]

1	sumaDivisoresHasta :: Integer -> [(Integer,Integer)]

tal que (sumaDivisoresHasta n) es la lista de los pares (a,b) tales que a es un número entre 1 y n y b es la suma de los divisores propios de a. Por ejemplo,

   λ> sumaDivisoresHasta 12
   [(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]
   λ> last (sumaDivisoresHasta2 (10^7))
   (10000000,14902280)

λ> sumaDivisoresHasta 12

[(1,0),(2,1),(3,1),(4,3),(5,1),(6,6),(7,1),(8,7),(9,4),(10,8),(11,1),(12,16)]

λ> last (sumaDivisoresHasta2 (10^7))

(10000000,14902280)

Soluciones

import Data.Array (accumArray, assocs)
import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

-- 1ª solución
-- ===========

sumaDivisoresHasta1 :: Integer -> [(Integer,Integer)]
sumaDivisoresHasta1 n = [(x, sum (divisores x)) | x <- [1..n]]

divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución
-- ===========

sumaDivisoresHasta2 :: Integer -> [(Integer,Integer)]
sumaDivisoresHasta2 n = [(x, sumaDivisores x) | x <- [1..n]]

sumaDivisores :: Integer -> Integer
sumaDivisores x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Integer)
primeroYlongitud (x:xs) =
  (x, 1 + genericLength xs)

-- 3ª solución
-- ===========

sumaDivisoresHasta3 :: Integer -> [(Integer,Integer)]
sumaDivisoresHasta3 n = assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a es
-- un número entre 2 y n y b es un divisor propio e x. Por ejemplo,
--    λ> divisoresHasta 6
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]
--    λ> divisoresHasta 8
--    [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]
divisoresHasta :: Integer -> [(Integer,Integer)]
divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_sumaDivisoresHasta :: Positive Integer -> Bool
prop_sumaDivisoresHasta (Positive n) =
  all (== sumaDivisoresHasta1 n)
      [ sumaDivisoresHasta2 n
      , sumaDivisoresHasta3 n
      ]

-- La comprobación es
--    λ> quickCheck prop_sumaDivisoresHasta
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> last (sumaDivisoresHasta1 (10^6))
--    (1000000,1480437)
--    (0.47 secs, 308,542,392 bytes)
--    λ> last (sumaDivisoresHasta2 (10^6))
--    (1000000,2480437)
--    (0.26 secs, 208,548,944 bytes)
--    λ> last (sumaDivisoresHasta3 (10^6))
--    (1000000,1480437)
--    (6.65 secs, 3,249,831,856 bytes)
--
--    λ> last (sumaDivisoresHasta1 (5*10^6))
--    (5000000,7402312)
--    (2.27 secs, 1,540,543,352 bytes)
--    λ> last (sumaDivisoresHasta2 (5*10^6))
--    (5000000,12402312)
--    (1.19 secs, 1,040,549,800 bytes)

import Data.Array (accumArray, assocs)

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

-- ===========

sumaDivisoresHasta1 :: Integer -> [(Integer,Integer)]

sumaDivisoresHasta1 n = [(x, sum (divisores x)) | x <- [1..n]]

divisores :: Integer -> [Integer]

divisores n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución

-- ===========

sumaDivisoresHasta2 :: Integer -> [(Integer,Integer)]

sumaDivisoresHasta2 n = [(x, sumaDivisores x) | x <- [1..n]]

sumaDivisores :: Integer -> Integer

sumaDivisores x =

product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la

-- descomposición prima de x. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs

-- y la longitud de xs. Por ejemplo,

-- primeroYlongitud [3,2,5,7] == (3,4)

primeroYlongitud :: [a] -> (a,Integer)

primeroYlongitud (x:xs) =

(x, 1 + genericLength xs)

-- 3ª solución

-- ===========

sumaDivisoresHasta3 :: Integer -> [(Integer,Integer)]

sumaDivisoresHasta3 n = assocs (accumArray (+) 0 (1,n) (divisoresHasta n))

-- (divisoresHasta n) es la lista de los pares (a,b) tales que a es

-- un número entre 2 y n y b es un divisor propio e x. Por ejemplo,

-- λ> divisoresHasta 6

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]

-- λ> divisoresHasta 8

-- [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]

divisoresHasta :: Integer -> [(Integer,Integer)]

divisoresHasta n = [(a,b) | b <- [1..n `div` 2], a <- [b*2, b*3..n]]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_sumaDivisoresHasta :: Positive Integer -> Bool

prop_sumaDivisoresHasta (Positive n) =

all (== sumaDivisoresHasta1 n)

[ sumaDivisoresHasta2 n

, sumaDivisoresHasta3 n

]

-- La comprobación es

-- λ> quickCheck prop_sumaDivisoresHasta

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> last (sumaDivisoresHasta1 (10^6))

-- (1000000,1480437)

-- (0.47 secs, 308,542,392 bytes)

-- λ> last (sumaDivisoresHasta2 (10^6))

-- (1000000,2480437)

-- (0.26 secs, 208,548,944 bytes)

-- λ> last (sumaDivisoresHasta3 (10^6))

-- (1000000,1480437)

-- (6.65 secs, 3,249,831,856 bytes)

-- λ> last (sumaDivisoresHasta1 (5*10^6))

-- (5000000,7402312)

-- (2.27 secs, 1,540,543,352 bytes)

-- λ> last (sumaDivisoresHasta2 (5*10^6))

-- (5000000,12402312)

-- (1.19 secs, 1,040,549,800 bytes)

El código se encuentra en GitHub.

Parejas de números y divisores

PorJosé A. Alonso 19-mayo-2022

Definir la función

   divisoresHasta :: Int -> [(Int,Int)]

1	divisoresHasta :: Int -> [(Int,Int)]

tal que (divisoresHasta n) es la lista de los pares (a,b) tales que a es un número entre 2 y n y b es un divisor propio de a. Por ejemplo,

   λ> divisoresHasta 6
   [(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]
   λ> divisoresHasta 8
   [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]
   λ> length (divisoresHasta 1234567)
   16272448

λ> divisoresHasta 6

[(2,1),(3,1),(4,1),(5,1),(6,1),(4,2),(6,2),(6,3)]

λ> divisoresHasta 8

[(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),(8,1),(4,2),(6,2),(8,2),(6,3),(8,4)]

λ> length (divisoresHasta 1234567)

16272448

Soluciones

import Data.List (sort, sortBy)
import Data.Ord (comparing)
import Data.Tuple (swap)
import Test.QuickCheck

-- 1ª solución
-- ===========

divisoresHasta1 :: Int -> [(Int,Int)]
divisoresHasta1 n =
  ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])
  where ordena ps = [intercambia p | p <- sort (map intercambia ps)]
        intercambia (x,y) = (y,x)

divisoresPropios :: Int -> [Int]
divisoresPropios n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución
-- ===========

divisoresHasta2 :: Int -> [(Int,Int)]
divisoresHasta2 n =
  ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])
  where ordena           = sortBy comp
        comp (x,y) (u,v) = compare (y,x) (v,u)

-- 3ª solución
-- ===========

divisoresHasta3 :: Int -> [(Int,Int)]
divisoresHasta3 n =
  ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])
  where ordena = sortBy (comparing swap)

-- 4ª solución
-- ===========

divisoresHasta4 :: Int -> [(Int,Int)]
divisoresHasta4 n =
  [(a,b) | b <- [1..n `div` 2], a <- [2*b, 3*b..n]]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_divisoresHasta :: Positive Int -> Bool
prop_divisoresHasta (Positive n) =
  all (== divisoresHasta1 n)
      [ divisoresHasta2 n
      , divisoresHasta3 n
      , divisoresHasta4 n
      ]

-- La comprobación es
--    λ> quickCheck prop_divisoresHasta
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (divisoresHasta1 4000)
--    29805
--    (1.78 secs, 843,584,472 bytes)
--    λ> length (divisoresHasta2 4000)
--    29805
--    (1.78 secs, 879,421,056 bytes)
--    λ> length (divisoresHasta3 4000)
--    29805
--    (1.70 secs, 876,475,584 bytes)
--    λ> length (divisoresHasta4 4000)
--    29805
--    (0.02 secs, 6,332,016 bytes)

import Data.List (sort, sortBy)

import Data.Ord (comparing)

import Data.Tuple (swap)

import Test.QuickCheck

-- 1ª solución

-- ===========

divisoresHasta1 :: Int -> [(Int,Int)]

divisoresHasta1 n =

ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])

where ordena ps = [intercambia p | p <- sort (map intercambia ps)]

intercambia (x,y) = (y,x)

divisoresPropios :: Int -> [Int]

divisoresPropios n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución

-- ===========

divisoresHasta2 :: Int -> [(Int,Int)]

divisoresHasta2 n =

ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])

where ordena = sortBy comp

comp (x,y) (u,v) = compare (y,x) (v,u)

-- 3ª solución

-- ===========

divisoresHasta3 :: Int -> [(Int,Int)]

divisoresHasta3 n =

ordena (concat [[(a,b) | b <- divisoresPropios a] | a <- [2..n]])

where ordena = sortBy (comparing swap)

-- 4ª solución

-- ===========

divisoresHasta4 :: Int -> [(Int,Int)]

divisoresHasta4 n =

[(a,b) | b <- [1..n `div` 2], a <- [2*b, 3*b..n]]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_divisoresHasta :: Positive Int -> Bool

prop_divisoresHasta (Positive n) =

all (== divisoresHasta1 n)

[ divisoresHasta2 n

, divisoresHasta3 n

, divisoresHasta4 n

]

-- La comprobación es

-- λ> quickCheck prop_divisoresHasta

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (divisoresHasta1 4000)

-- 29805

-- (1.78 secs, 843,584,472 bytes)

-- λ> length (divisoresHasta2 4000)

-- 29805

-- (1.78 secs, 879,421,056 bytes)

-- λ> length (divisoresHasta3 4000)

-- 29805

-- (1.70 secs, 876,475,584 bytes)

-- λ> length (divisoresHasta4 4000)

-- 29805

-- (0.02 secs, 6,332,016 bytes)

El código se encuentra en GitHub.

Sumas de 4 primos

PorJosé A. Alonso 18-mayo-202218-mayo-2022

La conjetura de Waring sobre los números primos establece que todo número impar es primo o la suma de tres primos. La conjetura de Goldbach afirma que todo par mayor que 2 es la suma de dos números primos. Ambos ha estado abiertos durante más de 200 años. En este problema no se propone su solución, sino una tarea más simple: buscar una manera de expresar los enteros mayores que 7 como suma de exactamente cuatro números primos; es decir, definir la función

   suma4primos :: Integer -> [(Integer,Integer,Integer,Integer)]

1	suma4primos :: Integer -> [(Integer,Integer,Integer,Integer)]

tal que (suma4primos n) es la lista de las cuádruplas crecientes (a,b,c,d) de números primos cuya suma es n (que se supone mayor que 7). Por ejemplo,

   suma4primos 18             == [(2,2,3,11),(2,2,7,7),(3,3,5,7),(3,5,5,5)]
   head (suma4primos (10^14)) == (2,2,23,99999999999973)

1 2	suma4primos 18 == [(2,2,3,11),(2,2,7,7),(3,3,5,7),(3,5,5,5)] head (suma4primos (10^14)) == (2,2,23,99999999999973)

Comprobar con QuickCheck que todo entero mayor que 7 se puede escribir como suma de exactamente cuatro números primos.

Soluciones

import Data.Numbers.Primes (isPrime, primes)
import Test.QuickCheck

-- 1ª solución
-- ===========

suma4primos1 :: Integer -> [(Integer, Integer, Integer, Integer)]
suma4primos1 n =
  [(a,b,c,d) | let as = takeWhile (< n) primes,
               a <- as,
               let bs = takeWhile (< n-a) as,
               b <- bs, a <= b,
               let cs = takeWhile (< n-a-b) bs,
               c <- cs, b <= c,
               let d = n-a-b-c, c <= d,
               isPrime d]

-- 2ª solución
-- ===========

suma4primos2 :: Integer -> [(Integer, Integer, Integer, Integer)]
suma4primos2 n =
  [(a,b,c,d) | let as = takeWhile (< n) primes,
               a <- as,
               let bs = takeWhile (< n-a) (dropWhile (< a) as),
               b <- bs,
               let cs = takeWhile (<n-a-b) (dropWhile (< b) bs),
               c <- cs,
               let d = n-a-b-c,
               c <= d,
               isPrime d]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_suma4primos :: Positive Integer -> Bool
prop_suma4primos (Positive n) =
  suma4primos1 n == suma4primos2 n

-- La comprobación es
--    λ> quickCheck prop_suma4primos
--    +++ OK, passed 100 tests; 526 discarded.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> length (suma4primos1 2000)
--    90219
--    (2.98 secs, 4,517,620,744 bytes)
--    λ> length (suma4primos2 2000)
--    90219
--    (2.22 secs, 4,223,251,928 bytes)
--
--    λ> head (suma4primos1 (10^14))
--    (2,2,23,99999999999973)
--    (1.67 secs, 5,963,327,168 bytes)
--    λ> head (suma4primos2 (10^14))
--    (2,2,23,99999999999973)
--    (1.70 secs, 5,963,326,848 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_suma4primos2 :: Integer -> Property
prop_suma4primos2 n =
  n > 7 ==> not (null (suma4primos1 n))

-- La comprobación es
--    λ> quickCheck prop_suma4primos2
--    +++ OK, passed 100 tests; 582 discarded.

import Data.Numbers.Primes (isPrime, primes)

import Test.QuickCheck

-- 1ª solución

-- ===========

suma4primos1 :: Integer -> [(Integer, Integer, Integer, Integer)]

suma4primos1 n =

[(a,b,c,d) | let as = takeWhile (< n) primes,

a <- as,

let bs = takeWhile (< n-a) as,

b <- bs, a <= b,

let cs = takeWhile (< n-a-b) bs,

c <- cs, b <= c,

let d = n-a-b-c, c <= d,

isPrime d]

-- 2ª solución

-- ===========

suma4primos2 :: Integer -> [(Integer, Integer, Integer, Integer)]

suma4primos2 n =

[(a,b,c,d) | let as = takeWhile (< n) primes,

a <- as,

let bs = takeWhile (< n-a) (dropWhile (< a) as),

b <- bs,

let cs = takeWhile (<n-a-b) (dropWhile (< b) bs),

c <- cs,

let d = n-a-b-c,

c <= d,

isPrime d]

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_suma4primos :: Positive Integer -> Bool

prop_suma4primos (Positive n) =

suma4primos1 n == suma4primos2 n

-- La comprobación es

-- λ> quickCheck prop_suma4primos

-- +++ OK, passed 100 tests; 526 discarded.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> length (suma4primos1 2000)

-- 90219

-- (2.98 secs, 4,517,620,744 bytes)

-- λ> length (suma4primos2 2000)

-- 90219

-- (2.22 secs, 4,223,251,928 bytes)

-- λ> head (suma4primos1 (10^14))

-- (2,2,23,99999999999973)

-- (1.67 secs, 5,963,327,168 bytes)

-- λ> head (suma4primos2 (10^14))

-- (2,2,23,99999999999973)

-- (1.70 secs, 5,963,326,848 bytes)

-- Comprobación de la propiedad

-- ============================

-- La propiedad es

prop_suma4primos2 :: Integer -> Property

prop_suma4primos2 n =

n > 7 ==> not (null (suma4primos1 n))

-- La comprobación es

-- λ> quickCheck prop_suma4primos2

-- +++ OK, passed 100 tests; 582 discarded.

El código se encuentra en GitHub.

Sucesión de sumas de dos números abundantes

PorJosé A. Alonso 17-mayo-202217-mayo-2022

Un número n es abundante si la suma de los divisores propios de n es mayor que n. El primer número abundante es el 12 (cuyos divisores propios son 1, 2, 3, 4 y 6 cuya suma es 16). Por tanto, el menor número que es la suma de dos números abundantes es el 24.

Definir la sucesión

   sumasDeDosAbundantes :: [Integer]

1	sumasDeDosAbundantes :: [Integer]

cuyos elementos son los números que se pueden escribir como suma de dos números abundantes. Por ejemplo,

   take 10 sumasDeDosAbundantes  ==  [24,30,32,36,38,40,42,44,48,50]

1	take 10 sumasDeDosAbundantes == [24,30,32,36,38,40,42,44,48,50]

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

-- 1ª solución
-- ===========

sumasDeDosAbundantes1 :: [Integer]
sumasDeDosAbundantes1 = [n | n <- [1..], esSumaDeDosAbundantes n]

-- (esSumaDeDosAbundantes n) se verifica si n es suma de dos números
-- abundantes. Por ejemplo,
--    esSumaDeDosAbundantes 24           ==  True
--    any esSumaDeDosAbundantes [1..22]  ==  False
esSumaDeDosAbundantes :: Integer -> Bool
esSumaDeDosAbundantes n = (not . null) [x | x <- xs, n-x `elem` xs]
  where xs = takeWhile (<n) abundantes

-- abundantes es la lista de los números abundantes. Por ejemplo,
--    take 10 abundantes  ==  [12,18,20,24,30,36,40,42,48,54]
abundantes :: [Integer]
abundantes = [n | n <- [2..], abundante n]

-- (abundante n) se verifica si n es abundante. Por ejemplo,
--    abundante 12  ==  True
--    abundante 11  ==  False
abundante :: Integer -> Bool
abundante n = sum (divisores n) > n

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,
--    divisores 12  ==  [1,2,3,4,6]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución
-- ===========

sumasDeDosAbundantes2 :: [Integer]
sumasDeDosAbundantes2 = filter esSumaDeDosAbundantes2 [1..]

esSumaDeDosAbundantes2 :: Integer -> Bool
esSumaDeDosAbundantes2 n = (not . null) [x | x <- xs, n-x `elem` xs]
  where xs = takeWhile (<n) abundantes2

abundantes2 :: [Integer]
abundantes2 = filter abundante2 [2..]

abundante2 :: Integer -> Bool
abundante2 n = sumaDivisores n > n

sumaDivisores :: Integer -> Integer
sumaDivisores x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Integer)
primeroYlongitud (x:xs) =
  (x, 1 + genericLength xs)

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_sumasDeDosAbundantes :: Positive Int -> Bool
prop_sumasDeDosAbundantes (Positive n) =
  sumasDeDosAbundantes1 !! n == sumasDeDosAbundantes2 !! n

-- La comprobación es
--    λ> quickCheck prop_sumasDeDosAbundantes
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> sumasDeDosAbundantes1 !! (2*10^3)
--    2887
--    (2.54 secs, 516,685,168 bytes)
--    λ> sumasDeDosAbundantes2 !! (2*10^3)
--    2887
--    (1.43 secs, 141,606,136 bytes)

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

-- ===========

sumasDeDosAbundantes1 :: [Integer]

sumasDeDosAbundantes1 = [n | n <- [1..], esSumaDeDosAbundantes n]

-- (esSumaDeDosAbundantes n) se verifica si n es suma de dos números

-- abundantes. Por ejemplo,

-- esSumaDeDosAbundantes 24 == True

-- any esSumaDeDosAbundantes [1..22] == False

esSumaDeDosAbundantes :: Integer -> Bool

esSumaDeDosAbundantes n = (not . null) [x | x <- xs, n-x `elem` xs]

where xs = takeWhile (<n) abundantes

-- abundantes es la lista de los números abundantes. Por ejemplo,

-- take 10 abundantes == [12,18,20,24,30,36,40,42,48,54]

abundantes :: [Integer]

abundantes = [n | n <- [2..], abundante n]

-- (abundante n) se verifica si n es abundante. Por ejemplo,

-- abundante 12 == True

-- abundante 11 == False

abundante :: Integer -> Bool

abundante n = sum (divisores n) > n

-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,

-- divisores 12 == [1,2,3,4,6]

divisores :: Integer -> [Integer]

divisores n = [x | x <- [1..n `div` 2], n `mod` x == 0]

-- 2ª solución

-- ===========

sumasDeDosAbundantes2 :: [Integer]

sumasDeDosAbundantes2 = filter esSumaDeDosAbundantes2 [1..]

esSumaDeDosAbundantes2 :: Integer -> Bool

esSumaDeDosAbundantes2 n = (not . null) [x | x <- xs, n-x `elem` xs]

where xs = takeWhile (<n) abundantes2

abundantes2 :: [Integer]

abundantes2 = filter abundante2 [2..]

abundante2 :: Integer -> Bool

abundante2 n = sumaDivisores n > n

sumaDivisores :: Integer -> Integer

sumaDivisores x =

product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x] - x

-- (factorizacion x) es la lista de las bases y exponentes de la

-- descomposición prima de x. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs

-- y la longitud de xs. Por ejemplo,

-- primeroYlongitud [3,2,5,7] == (3,4)

primeroYlongitud :: [a] -> (a,Integer)

primeroYlongitud (x:xs) =

(x, 1 + genericLength xs)

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_sumasDeDosAbundantes :: Positive Int -> Bool

prop_sumasDeDosAbundantes (Positive n) =

sumasDeDosAbundantes1 !! n == sumasDeDosAbundantes2 !! n

-- La comprobación es

-- λ> quickCheck prop_sumasDeDosAbundantes

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> sumasDeDosAbundantes1 !! (2*10^3)

-- 2887

-- (2.54 secs, 516,685,168 bytes)

-- λ> sumasDeDosAbundantes2 !! (2*10^3)

-- 2887

-- (1.43 secs, 141,606,136 bytes)

El código se encuentra en GitHub.

Suma de divisores

PorJosé A. Alonso 16-mayo-2022

Definir la función

   sumaDivisores :: Integer -> Integer

1	sumaDivisores :: Integer -> Integer

tal que (sumaDivisores x) es la suma de los divisores de x. Por ejemplo,

   sumaDivisores 12  ==  28
   sumaDivisores 25  ==  31
   sumaDivisores (product [1..25])  ==  93383273455325195473152000
   length (show (sumaDivisores (product [1..30000])))  ==  121289
   maximum (map sumaDivisores [1..10^5])  ==  403200

sumaDivisores 12 == 28

sumaDivisores 25 == 31

sumaDivisores (product [1..25]) == 93383273455325195473152000

length (show (sumaDivisores (product [1..30000]))) == 121289

maximum (map sumaDivisores [1..10^5]) == 403200

Soluciones

import Data.List (genericLength, group, inits)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

-- 1ª solución
-- ===========

sumaDivisores1 :: Integer -> Integer
sumaDivisores1 = sum . divisores

-- (divisores x) es la lista de los divisores de x. Por ejemplo,
--    divisores 60  ==  [1,5,3,15,2,10,6,30,4,20,12,60]
divisores :: Integer -> [Integer]
divisores = map (product . concat)
          . productoCartesiano
          . map inits
          . group
          . primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos
-- xss. Por ejemplo,
--    λ> producto [[1,3],[2,5],[6,4]]
--    [[1,2,6],[1,2,4],[1,5,6],[1,5,4],[3,2,6],[3,2,4],[3,5,6],[3,5,4]]
productoCartesiano :: [[a]] -> [[a]]
productoCartesiano []       = [[]]
productoCartesiano (xs:xss) =
  [x:ys | x <- xs, ys <- productoCartesiano xss]

-- 2ª solución
-- ===========

sumaDivisores2 :: Integer -> Integer
sumaDivisores2 = sum
               . map (product . concat)
               . sequence
               . map inits
               . group
               . primeFactors

-- 3ª solución
-- ===========

-- Si la descomposición de x en factores primos es
--    x = p(1)^e(1) . p(2)^e(2) . .... . p(n)^e(n)
-- entonces la suma de los divisores de x es
--    p(1)^(e(1)+1) - 1     p(2)^(e(2)+1) - 1       p(n)^(e(2)+1) - 1
--   ------------------- . ------------------- ... -------------------
--        p(1)-1                p(2)-1                  p(n)-1
-- Ver la demostración en http://bit.ly/2zUXZPc

sumaDivisores3 :: Integer -> Integer
sumaDivisores3 x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Integer)
primeroYlongitud (x:xs) =
  (x, 1 + genericLength xs)

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_sumaDivisores :: Positive Integer -> Bool
prop_sumaDivisores (Positive x) =
  all (== sumaDivisores1 x)
      [ sumaDivisores2 x
      , sumaDivisores3 x
      ]

-- La comprobación es
--    λ> quickCheck prop_sumaDivisores
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--   λ> sumaDivisores1 251888923423315469521109880000000
--   1471072204661054993275791673480320
--   (10.63 secs, 10,614,618,080 bytes)
--   λ> sumaDivisores2 251888923423315469521109880000000
--   1471072204661054993275791673480320
--   (2.51 secs, 5,719,399,056 bytes)
--   λ> sumaDivisores3 251888923423315469521109880000000
--   1471072204661054993275791673480320
--   (0.01 secs, 177,480 bytes)

import Data.List (genericLength, group, inits)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

-- ===========

sumaDivisores1 :: Integer -> Integer

sumaDivisores1 = sum . divisores

-- (divisores x) es la lista de los divisores de x. Por ejemplo,

-- divisores 60 == [1,5,3,15,2,10,6,30,4,20,12,60]

divisores :: Integer -> [Integer]

divisores = map (product . concat)

. productoCartesiano

. map inits

. group

. primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos

-- xss. Por ejemplo,

-- λ> producto [[1,3],[2,5],[6,4]]

-- [[1,2,6],[1,2,4],[1,5,6],[1,5,4],[3,2,6],[3,2,4],[3,5,6],[3,5,4]]

productoCartesiano :: [[a]] -> [[a]]

productoCartesiano [] = [[]]

productoCartesiano (xs:xss) =

[x:ys | x <- xs, ys <- productoCartesiano xss]

-- 2ª solución

-- ===========

sumaDivisores2 :: Integer -> Integer

sumaDivisores2 = sum

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

-- 3ª solución

-- ===========

-- Si la descomposición de x en factores primos es

-- x = p(1)^e(1) . p(2)^e(2) . .... . p(n)^e(n)

-- entonces la suma de los divisores de x es

-- p(1)^(e(1)+1) - 1 p(2)^(e(2)+1) - 1 p(n)^(e(2)+1) - 1

-- ------------------- . ------------------- ... -------------------

-- p(1)-1 p(2)-1 p(n)-1

-- Ver la demostración en http://bit.ly/2zUXZPc

sumaDivisores3 :: Integer -> Integer

sumaDivisores3 x =

product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de la

-- descomposición prima de x. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion = map primeroYlongitud . group . primeFactors

-- (primeroYlongitud xs) es el par formado por el primer elemento de xs

-- y la longitud de xs. Por ejemplo,

-- primeroYlongitud [3,2,5,7] == (3,4)

primeroYlongitud :: [a] -> (a,Integer)

primeroYlongitud (x:xs) =

(x, 1 + genericLength xs)

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_sumaDivisores :: Positive Integer -> Bool

prop_sumaDivisores (Positive x) =

all (== sumaDivisores1 x)

[ sumaDivisores2 x

, sumaDivisores3 x

]

-- La comprobación es

-- λ> quickCheck prop_sumaDivisores

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> sumaDivisores1 251888923423315469521109880000000

-- 1471072204661054993275791673480320

-- (10.63 secs, 10,614,618,080 bytes)

-- λ> sumaDivisores2 251888923423315469521109880000000

-- 1471072204661054993275791673480320

-- (2.51 secs, 5,719,399,056 bytes)

-- λ> sumaDivisores3 251888923423315469521109880000000

-- 1471072204661054993275791673480320

-- (0.01 secs, 177,480 bytes)

El código se encuentra en GitHub.

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico