product – Exercitium

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.

Número de divisores

PorJosé A. Alonso 13-mayo-2022

Definir la función

   numeroDivisores :: Integer -> Integer

1	numeroDivisores :: Integer -> Integer

tal que (numeroDivisores x) es el número de divisores de x. Por ejemplo,

   numeroDivisores 12  ==  6
   numeroDivisores 25  ==  3
   length (show (numeroDivisores (product [1..3*10^4])))  ==  1948

numeroDivisores 12 == 6

numeroDivisores 25 == 3

length (show (numeroDivisores (product [1..3*10^4]))) == 1948

Soluciones

import Data.List (genericLength, group, inits)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

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

numeroDivisores1 :: Integer -> Integer
numeroDivisores1 x =
  genericLength [y | y <- [1..x], x `mod` y == 0]

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

numeroDivisores2 :: Integer -> Integer
numeroDivisores2 1 = 1
numeroDivisores2 x
  | esCuadrado x = 2 * genericLength [y | y <- [1..raizEntera x], x `mod` y == 0] - 1
  | otherwise    = 2 * genericLength [y | y <- [1..raizEntera x], x `mod` y == 0]

-- (raizEntera x) es el mayor número entero cuyo cuadrado es menor o
-- igual que x. Por ejemplo,
--    raizEntera 3  ==  1
--    raizEntera 4  ==  2
--    raizEntera 5  ==  2
--    raizEntera 8  ==  2
--    raizEntera 9  ==  3
raizEntera :: Integer -> Integer
raizEntera x = floor (sqrt (fromInteger x))

-- (esCuadrado x) se verifica si x es un cuadrado perfecto. Por ejemplo,
--    esCuadrado 9  ==  True
--    esCuadrado 7  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x =
  x == (raizEntera x)^2

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

numeroDivisores3 :: Integer -> Integer
numeroDivisores3 =
  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 = map (product . concat)
          . productoCartesiano
          . map inits
          . group
          . primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos
-- xss. Por ejemplo,
--    λ> productoCartesiano [[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]

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

numeroDivisores4 :: Integer -> Integer
numeroDivisores4 = genericLength
                 . map (product . concat)
                 . sequence
                 . map inits
                 . group
                 . primeFactors

-- 5ª solución
-- ===========

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

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

-- La propiedad es
prop_numeroDivisores :: Positive Integer -> Bool
prop_numeroDivisores (Positive x) =
  all (== numeroDivisores1 x)
      [ numeroDivisores2 x
      , numeroDivisores3 x
      , numeroDivisores4 x
      , numeroDivisores5 x]

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

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

-- La comparación es
--    λ> numeroDivisores1 (product [1..10])
--    270
--    (1.67 secs, 726,327,208 bytes)
--    λ> numeroDivisores2 (product [1..10])
--    270
--    (0.01 secs, 929,000 bytes)
--
--    λ> numeroDivisores2 (product [1..16])
--    5376
--    (2.10 secs, 915,864,664 bytes)
--    λ> numeroDivisores3 (product [1..16])
--    5376
--    (0.01 secs, 548,472 bytes)
--
--    λ> numeroDivisores3 (product [1..30])
--    2332800
--    (3.80 secs, 4,149,811,688 bytes)
--    λ> numeroDivisores4 (product [1..30])
--    2332800
--    (0.59 secs, 722,253,848 bytes)
--    λ> numeroDivisores5 (product [1..30])
--    2332800
--    (0.00 secs, 587,856 bytes)

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, inits)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

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

numeroDivisores1 :: Integer -> Integer

numeroDivisores1 x =

genericLength [y | y <- [1..x], x `mod` y == 0]

-- 2ª solución

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

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 1 = 1

numeroDivisores2 x

| esCuadrado x = 2 * genericLength [y | y <- [1..raizEntera x], x `mod` y == 0] - 1

| otherwise = 2 * genericLength [y | y <- [1..raizEntera x], x `mod` y == 0]

-- (raizEntera x) es el mayor número entero cuyo cuadrado es menor o

-- igual que x. Por ejemplo,

-- raizEntera 3 == 1

-- raizEntera 4 == 2

-- raizEntera 5 == 2

-- raizEntera 8 == 2

-- raizEntera 9 == 3

raizEntera :: Integer -> Integer

raizEntera x = floor (sqrt (fromInteger x))

-- (esCuadrado x) se verifica si x es un cuadrado perfecto. Por ejemplo,

-- esCuadrado 9 == True

-- esCuadrado 7 == False

esCuadrado :: Integer -> Bool

esCuadrado x =

x == (raizEntera x)^2

-- 3ª solución

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

numeroDivisores3 :: Integer -> Integer

numeroDivisores3 =

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 = map (product . concat)

. productoCartesiano

. map inits

. group

. primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos

-- xss. Por ejemplo,

-- λ> productoCartesiano [[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]

-- 4ª solución

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

numeroDivisores4 :: Integer -> Integer

numeroDivisores4 = genericLength

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

-- 5ª solución

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

numeroDivisores5 :: Integer -> Integer

numeroDivisores5 =

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_numeroDivisores :: Positive Integer -> Bool

prop_numeroDivisores (Positive x) =

all (== numeroDivisores1 x)

[ numeroDivisores2 x

, numeroDivisores3 x

, numeroDivisores4 x

, numeroDivisores5 x]

-- La comprobación es

-- λ> quickCheck prop_numeroDivisores

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> numeroDivisores1 (product [1..10])

-- 270

-- (1.67 secs, 726,327,208 bytes)

-- λ> numeroDivisores2 (product [1..10])

-- 270

-- (0.01 secs, 929,000 bytes)

-- λ> numeroDivisores2 (product [1..16])

-- 5376

-- (2.10 secs, 915,864,664 bytes)

-- λ> numeroDivisores3 (product [1..16])

-- 5376

-- (0.01 secs, 548,472 bytes)

-- λ> numeroDivisores3 (product [1..30])

-- 2332800

-- (3.80 secs, 4,149,811,688 bytes)

-- λ> numeroDivisores4 (product [1..30])

-- 2332800

-- (0.59 secs, 722,253,848 bytes)

-- λ> numeroDivisores5 (product [1..30])

-- 2332800

-- (0.00 secs, 587,856 bytes)

El código se encuentra en GitHub.

Conjunto de divisores

PorJosé A. Alonso 12-mayo-202212-mayo-2022

Definir la función

   divisores :: Integer -> [Integer]

1	divisores :: Integer -> [Integer]

tal que (divisores x) es el conjunto de divisores de x. Por ejemplo,

  divisores 30  ==  [1,2,3,5,6,10,15,30]
  length (divisores (product [1..10]))  ==  270
  length (divisores (product [1..25]))  ==  340032

divisores 30 == [1,2,3,5,6,10,15,30]

length (divisores (product [1..10])) == 270

length (divisores (product [1..25])) == 340032

Soluciones

import Data.List (group, inits, nub, sort, subsequences)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

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

divisores1 :: Integer -> [Integer]
divisores1 n = [x | x <- [1..n], n `rem` x == 0]

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

divisores2 :: Integer -> [Integer]
divisores2 n = filter ((== 0) . mod n) [1..n]

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

divisores3 :: Integer -> [Integer]
divisores3 =
  nub . sort . map product . subsequences . primeFactors

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

divisores4 :: Integer -> [Integer]
divisores4 =
  sort
  . map (product . concat)
  . productoCartesiano
  . map inits
  . group
  . primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos
-- xss. Por ejemplo,
--    λ> productoCartesiano [[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]

-- 5ª solución
-- ===========

divisores5 :: Integer -> [Integer]
divisores5 = sort
           . map (product . concat)
           . sequence
           . map inits
           . group
           . primeFactors

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

-- La propiedad es
prop_divisores :: Positive Integer -> Bool
prop_divisores (Positive n) =
  all (== divisores1 n)
      [ divisores2 n
      , divisores3 n
      , divisores4 n
      , divisores5 n
      ]

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

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

--    λ> length (divisores (product [1..11]))
--    540
--    (12.51 secs, 7,983,499,736 bytes)
--    λ> length (divisores2 (product [1..11]))
--    540
--    (4.81 secs, 4,790,146,656 bytes)
--    λ> length (divisores3 (product [1..11]))
--    540
--    (0.10 secs, 107,339,848 bytes)
--    λ> length (divisores4 (product [1..11]))
--    540
--    (0.02 secs, 1,702,616 bytes)
--    λ> length (divisores5 (product [1..11]))
--    540
--    (0.02 secs, 1,205,824 bytes)
--
--    λ> length (divisores3 (product [1..14]))
--    2592
--    (7.89 secs, 9,378,454,912 bytes)
--    λ> length (divisores4 (product [1..14]))
--    2592
--    (0.03 secs, 9,426,528 bytes)
--    λ> length (divisores5 (product [1..14]))
--    2592
--    (0.02 secs, 6,636,608 bytes)
--    
--    λ> length (divisores4 (product [1..25]))
--    340032
--    (1.65 secs, 2,055,558,208 bytes)
--    λ> length (divisores5 (product [1..25]))
--    340032
--    (0.88 secs, 1,532,515,304 bytes)

100

101

102

103

104

105

106

107

import Data.List (group, inits, nub, sort, subsequences)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

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

divisores1 :: Integer -> [Integer]

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

-- 2ª solución

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

divisores2 :: Integer -> [Integer]

divisores2 n = filter ((== 0) . mod n) [1..n]

-- 3ª solución

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

divisores3 :: Integer -> [Integer]

divisores3 =

nub . sort . map product . subsequences . primeFactors

-- 4ª solución

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

divisores4 :: Integer -> [Integer]

divisores4 =

sort

. map (product . concat)

. productoCartesiano

. map inits

. group

. primeFactors

-- (productoCartesiano xss) es el producto cartesiano de los conjuntos

-- xss. Por ejemplo,

-- λ> productoCartesiano [[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]

-- 5ª solución

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

divisores5 :: Integer -> [Integer]

divisores5 = sort

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

-- Comprobación de equivalencia

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

-- La propiedad es

prop_divisores :: Positive Integer -> Bool

prop_divisores (Positive n) =

all (== divisores1 n)

[ divisores2 n

, divisores3 n

, divisores4 n

, divisores5 n

]

-- La comprobación es

-- λ> quickCheck prop_divisores

-- +++ OK, passed 100 tests.

-- Comparación de la eficiencia

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

-- λ> length (divisores (product [1..11]))

-- 540

-- (12.51 secs, 7,983,499,736 bytes)

-- λ> length (divisores2 (product [1..11]))

-- 540

-- (4.81 secs, 4,790,146,656 bytes)

-- λ> length (divisores3 (product [1..11]))

-- 540

-- (0.10 secs, 107,339,848 bytes)

-- λ> length (divisores4 (product [1..11]))

-- 540

-- (0.02 secs, 1,702,616 bytes)

-- λ> length (divisores5 (product [1..11]))

-- 540

-- (0.02 secs, 1,205,824 bytes)

-- λ> length (divisores3 (product [1..14]))

-- 2592

-- (7.89 secs, 9,378,454,912 bytes)

-- λ> length (divisores4 (product [1..14]))

-- 2592

-- (0.03 secs, 9,426,528 bytes)

-- λ> length (divisores5 (product [1..14]))

-- 2592

-- (0.02 secs, 6,636,608 bytes)

-- λ> length (divisores4 (product [1..25]))

-- 340032

-- (1.65 secs, 2,055,558,208 bytes)

-- λ> length (divisores5 (product [1..25]))

-- 340032

-- (0.88 secs, 1,532,515,304 bytes)

El código se encuentra en GitHub.

Mayor producto de las ramas de un árbol

PorJosé A. Alonso 28-marzo-20224-abril-2022

Los árboles se pueden representar mediante el siguiente tipo de datos

   data Arbol a = N a [Arbol a]
     deriving Show

1 2	data Arbol a = N a [Arbol a] deriving Show

Por ejemplo, los árboles

      1              3
    /  \            /|\
   2   3           / | \
       |          5  4  7
       4          |     /\
                  6    2  1

1 3

/ \ /|\

2 3 / | \

| 5 4 7

4 | /\

6 2 1

se representan por

   ej1, ej2 :: Arbol Int
   ej1 = N 1 [N 2 [],N 3 [N 4 []]]
   ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

ej1, ej2 :: Arbol Int

ej1 = N 1 [N 2 [],N 3 [N 4 []]]

ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

Definir la función

   mayorProducto :: (Ord a, Num a) => Arbol a -> a

1	mayorProducto :: (Ord a, Num a) => Arbol a -> a

tal que (mayorProducto a) es el mayor producto de las ramas del árbol a. Por ejemplo,

   λ> mayorProducto (N 1 [N 2 [], N  3 []])
   3
   λ> mayorProducto (N 1 [N 8 [], N  4 [N 3 []]])
   12
   λ> mayorProducto (N 1 [N 2 [],N 3 [N 4 []]])
   12
   λ> mayorProducto (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])
   90
   λ> mayorProducto (N (-8) [N 0 [N (-9) []],N 6 []])
   0
   λ> a = N (-4) [N (-7) [],N 14 [N 19 []],N (-1) [N (-6) [],N 21 []],N (-4) []]
   λ> mayorProducto a
   84

λ> mayorProducto (N 1 [N 2 [], N 3 []])

λ> mayorProducto (N 1 [N 8 [], N 4 [N 3 []]])

λ> mayorProducto (N 1 [N 2 [],N 3 [N 4 []]])

λ> mayorProducto (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])

λ> mayorProducto (N (-8) [N 0 [N (-9) []],N 6 []])

λ> a = N (-4) [N (-7) [],N 14 [N 19 []],N (-1) [N (-6) [],N 21 []],N (-4) []]

λ> mayorProducto a

Soluciones

import Test.QuickCheck

data Arbol a = N a [Arbol a]
  deriving Show

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

mayorProducto1 :: (Ord a, Num a) => Arbol a -> a
mayorProducto1 a = maximum [product xs | xs <- ramas a]

-- (ramas a) es la lista de las ramas del árbol a. Por ejemplo,
--    λ> ramas (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])
--    [[3,5,6],[3,4],[3,7,2],[3,7,1]]
ramas :: Arbol b -> [[b]]
ramas (N x []) = [[x]]
ramas (N x as) = [x : xs | a <- as, xs <- ramas a]

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

mayorProducto2 :: (Ord a, Num a) => Arbol a -> a
mayorProducto2 a = maximum (map product (ramas a))

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

mayorProducto3 :: (Ord a, Num a) => Arbol a -> a
mayorProducto3 = maximum . map product . ramas

-- 4º solución
-- ===========

mayorProducto4 :: (Ord a, Num a) => Arbol a -> a
mayorProducto4 = maximum . productosRamas

-- (productosRamas a) es la lista de los productos de las ramas
-- del árbol a. Por ejemplo,
--    λ> productosRamas (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])
--    [90,12,42,21]
productosRamas :: (Ord a, Num a) => Arbol a -> [a]
productosRamas (N x []) = [x]
productosRamas (N x xs) = [x * y | a <- xs, y <- productosRamas a]

-- 5ª solución
-- ===========

mayorProducto5 :: (Ord a, Num a) => Arbol a -> a
mayorProducto5 (N x []) = x
mayorProducto5 (N x xs)
  | x > 0     = x * maximum (map mayorProducto5 xs)
  | x == 0    = 0
  | otherwise = x * minimum (map menorProducto xs)

-- (menorProducto a) es el menor producto de las ramas del árbol
-- a. Por ejemplo,
--    λ> menorProducto (N 1 [N 2 [], N  3 []])
--    2
--    λ> menorProducto (N 1 [N 8 [], N  4 [N 3 []]])
--    8
--    λ> menorProducto (N 1 [N 2 [],N 3 [N 4 []]])
--    2
--    λ> menorProducto (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])
--    12
menorProducto :: (Ord a, Num a) => Arbol a -> a
menorProducto (N x []) = x
menorProducto (N x xs)
  | x > 0     = x * minimum (map menorProducto xs)
  | x == 0    = 0
  | otherwise = x * maximum (map mayorProducto2 xs)

-- 6ª solución
-- ===========

mayorProducto6 :: (Ord a, Num a) => Arbol a -> a
mayorProducto6 = maximum . aux
  where aux (N a []) = [a]
        aux (N a b)  = [v,u]
          where u = maximum g
                v = minimum g
                g = map (*a) (concatMap aux b)

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

-- (arbolArbitrario n) es un árbol aleatorio de orden n. Por ejemplo,
--   > sample (arbolArbitrario 5 :: Gen (Arbol Int))
--   N 0 [N 0 []]
--   N (-2) []
--   N 4 []
--   N 2 [N 4 []]
--   N 8 []
--   N (-2) [N (-9) [],N 7 []]
--   N 11 []
--   N (-11) [N 4 [],N 14 []]
--   N 10 [N (-3) [],N 13 []]
--   N 12 [N 11 []]
--   N 20 [N (-18) [],N (-13) []]
arbolArbitrario :: Arbitrary a => Int -> Gen (Arbol a)
arbolArbitrario n = do
  x  <- arbitrary
  ms <- sublistOf [0 .. n `div` 2]
  as <- mapM arbolArbitrario ms
  return (N x as)

-- Arbol es una subclase de Arbitraria
instance Arbitrary a => Arbitrary (Arbol a) where
  arbitrary = sized arbolArbitrario

-- La propiedad es
prop_mayorProducto :: Arbol Integer -> Bool
prop_mayorProducto a =
  all (== mayorProducto1 a)
      [f a | f <- [ mayorProducto2
                  , mayorProducto3
                  , mayorProducto4
                  , mayorProducto5
                  , mayorProducto6
                  ]]

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

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

-- La comparación es
--    λ> ejArbol <- generate (arbolArbitrario 600 :: Gen (Arbol Integer))
--    λ> mayorProducto1 ejArbol
--    2419727651266241493467136000
--    (1.87 secs, 1,082,764,480 bytes)
--    λ> mayorProducto2 ejArbol
--    2419727651266241493467136000
--    (1.57 secs, 1,023,144,008 bytes)
--    λ> mayorProducto3 ejArbol
--    2419727651266241493467136000
--    (1.55 secs, 1,023,144,248 bytes)
--    λ> mayorProducto4 ejArbol
--    2419727651266241493467136000
--    (1.60 secs, 824,473,800 bytes)
--    λ> mayorProducto5 ejArbol
--    2419727651266241493467136000
--    (0.83 secs, 732,370,352 bytes)
--    λ> mayorProducto6 ejArbol
--    2419727651266241493467136000
--    (0.98 secs, 817,473,344 bytes)
--
--    λ> ejArbol2 <- generate (arbolArbitrario 700 :: Gen (Arbol Integer))
--    λ> mayorProducto5 ejArbol2
--    1044758937398026715504640000000
--    (4.94 secs, 4,170,324,376 bytes)
--    λ> mayorProducto6 ejArbol2
--    1044758937398026715504640000000
--    (5.88 secs, 4,744,782,024 bytes)

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

import Test.QuickCheck

data Arbol a = N a [Arbol a]

deriving Show

-- 1ª solución

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

mayorProducto1 :: (Ord a, Num a) => Arbol a -> a

mayorProducto1 a = maximum [product xs | xs <- ramas a]

-- (ramas a) es la lista de las ramas del árbol a. Por ejemplo,

-- λ> ramas (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])

-- [[3,5,6],[3,4],[3,7,2],[3,7,1]]

ramas :: Arbol b -> [[b]]

ramas (N x []) = [[x]]

ramas (N x as) = [x : xs | a <- as, xs <- ramas a]

-- 2ª solución

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

mayorProducto2 :: (Ord a, Num a) => Arbol a -> a

mayorProducto2 a = maximum (map product (ramas a))

-- 3ª solución

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

mayorProducto3 :: (Ord a, Num a) => Arbol a -> a

mayorProducto3 = maximum . map product . ramas

-- 4º solución

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

mayorProducto4 :: (Ord a, Num a) => Arbol a -> a

mayorProducto4 = maximum . productosRamas

-- (productosRamas a) es la lista de los productos de las ramas

-- del árbol a. Por ejemplo,

-- λ> productosRamas (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])

-- [90,12,42,21]

productosRamas :: (Ord a, Num a) => Arbol a -> [a]

productosRamas (N x []) = [x]

productosRamas (N x xs) = [x * y | a <- xs, y <- productosRamas a]

-- 5ª solución

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

mayorProducto5 :: (Ord a, Num a) => Arbol a -> a

mayorProducto5 (N x []) = x

mayorProducto5 (N x xs)

| x > 0 = x * maximum (map mayorProducto5 xs)

| x == 0 = 0

| otherwise = x * minimum (map menorProducto xs)

-- (menorProducto a) es el menor producto de las ramas del árbol

-- a. Por ejemplo,

-- λ> menorProducto (N 1 [N 2 [], N 3 []])

-- 2

-- λ> menorProducto (N 1 [N 8 [], N 4 [N 3 []]])

-- 8

-- λ> menorProducto (N 1 [N 2 [],N 3 [N 4 []]])

-- 2

-- λ> menorProducto (N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]])

-- 12

menorProducto :: (Ord a, Num a) => Arbol a -> a

menorProducto (N x []) = x

menorProducto (N x xs)

| x > 0 = x * minimum (map menorProducto xs)

| x == 0 = 0

| otherwise = x * maximum (map mayorProducto2 xs)

-- 6ª solución

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

mayorProducto6 :: (Ord a, Num a) => Arbol a -> a

mayorProducto6 = maximum . aux

where aux (N a []) = [a]

aux (N a b) = [v,u]

where u = maximum g

v = minimum g

g = map (*a) (concatMap aux b)

-- Comprobación de equivalencia

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

-- (arbolArbitrario n) es un árbol aleatorio de orden n. Por ejemplo,

-- > sample (arbolArbitrario 5 :: Gen (Arbol Int))

-- N 0 [N 0 []]

-- N (-2) []

-- N 4 []

-- N 2 [N 4 []]

-- N 8 []

-- N (-2) [N (-9) [],N 7 []]

-- N 11 []

-- N (-11) [N 4 [],N 14 []]

-- N 10 [N (-3) [],N 13 []]

-- N 12 [N 11 []]

-- N 20 [N (-18) [],N (-13) []]

arbolArbitrario :: Arbitrary a => Int -> Gen (Arbol a)

arbolArbitrario n = do

x <- arbitrary

ms <- sublistOf [0 .. n `div` 2]

as <- mapM arbolArbitrario ms

return (N x as)

-- Arbol es una subclase de Arbitraria

instance Arbitrary a => Arbitrary (Arbol a) where

arbitrary = sized arbolArbitrario

-- La propiedad es

prop_mayorProducto :: Arbol Integer -> Bool

prop_mayorProducto a =

all (== mayorProducto1 a)

[f a | f <- [ mayorProducto2

, mayorProducto3

, mayorProducto4

, mayorProducto5

, mayorProducto6

]]

-- La comprobación es

-- λ> quickCheck prop_mayorProducto

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> ejArbol <- generate (arbolArbitrario 600 :: Gen (Arbol Integer))

-- λ> mayorProducto1 ejArbol

-- 2419727651266241493467136000

-- (1.87 secs, 1,082,764,480 bytes)

-- λ> mayorProducto2 ejArbol

-- 2419727651266241493467136000

-- (1.57 secs, 1,023,144,008 bytes)

-- λ> mayorProducto3 ejArbol

-- 2419727651266241493467136000

-- (1.55 secs, 1,023,144,248 bytes)

-- λ> mayorProducto4 ejArbol

-- 2419727651266241493467136000

-- (1.60 secs, 824,473,800 bytes)

-- λ> mayorProducto5 ejArbol

-- 2419727651266241493467136000

-- (0.83 secs, 732,370,352 bytes)

-- λ> mayorProducto6 ejArbol

-- 2419727651266241493467136000

-- (0.98 secs, 817,473,344 bytes)

-- λ> ejArbol2 <- generate (arbolArbitrario 700 :: Gen (Arbol Integer))

-- λ> mayorProducto5 ejArbol2

-- 1044758937398026715504640000000

-- (4.94 secs, 4,170,324,376 bytes)

-- λ> mayorProducto6 ejArbol2

-- 1044758937398026715504640000000

-- (5.88 secs, 4,744,782,024 bytes)

El código se encuentra en GitHub.

La elaboración de las soluciones se describe en el siguiente vídeo

Nuevas soluciones

En los comentarios se pueden escribir nuevas soluciones.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

La sucesión de Sylvester

PorJosé A. Alonso 6-mayo-202013-mayo-2020

La sucesión de Sylvester es la sucesión que comienza en 2 y sus restantes términos se obtienen multiplicando los anteriores y sumándole 1.

Definir las funciones

   sylvester        :: Integer -> Integer
   graficaSylvester :: Integer -> Integer -> IO ()

1 2	sylvester :: Integer -> Integer graficaSylvester :: Integer -> Integer -> IO ()

tales que

(sylvester n) es el n-ésimo término de la sucesión de Sylvester. Por ejemplo,

     λ> [sylvester n | n <- [0..7]]
     [2,3,7,43,1807,3263443,10650056950807,113423713055421844361000443]
     λ> length (show (sylvester 25))
     6830085

λ> [sylvester n | n <- [0..7]]

[2,3,7,43,1807,3263443,10650056950807,113423713055421844361000443]

λ> length (show (sylvester 25))

6830085

(graficaSylvester d n) dibuja la gráfica de los d últimos dígitos de los n primeros términos de la sucesión de Sylvester. Por ejemplo,
- (graficaSylvester 3 30) dibuja
- (graficaSylvester 4 30) dibuja
- (graficaSylvester 5 30) dibuja

Nota: Se puede usar programación dinámica para aumentar la eficiencia.

Soluciones

import Data.List               (genericIndex)
import Data.Array              ((!), array)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG))

-- 1ª solución (por recursión)
-- ===========================

sylvester1 :: Integer -> Integer
sylvester1 0 = 2
sylvester1 n = 1 + product [sylvester1 k | k <- [0..n-1]]

-- 2ª solución (con programación dinámica)
-- =======================================

sylvester2 :: Integer -> Integer
sylvester2 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + product [v!k | k <- [0..m-1]]

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

-- Observando que
--    S(n) = 1 + S(0)*S(1)*...*S(n-2)*S(n-1)
--         = 1 + (1 + S(0)*S(1)*...*S(n-2))*S(n-1) - S(n-1)
--         = 1 + S(n-1)*S(n-1) - S(n-1)
--         = 1 + S(n-1)^2 - S(n-1)
-- se obtiene la siguiente definición.
sylvester3 :: Integer -> Integer
sylvester3 0 = 2
sylvester3 n = 1 + x^2 - x
  where x = sylvester3 (n-1)

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

sylvester4 :: Integer -> Integer
sylvester4 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + x^2 - x
    where x = v ! (m-1)

-- 5ª solución
-- ===========

sylvester5 :: Integer -> Integer
sylvester5 n = sucSylvester5 `genericIndex` n

sucSylvester5 :: [Integer]
sucSylvester5 = iterate (\x -> (x-1)*x+1) 2 

-- La comparación es
--    λ> length (show (sylvester1 23))
--    1707522
--    (6.03 secs, 4,090,415,704 bytes)
--    λ> length (show (sylvester2 23))
--    1707522
--    (0.33 secs, 109,477,296 bytes)
--    λ> length (show (sylvester3 23))
--    1707522
--    (0.35 secs, 109,395,136 bytes)
--    λ> length (show (sylvester4 23))
--    1707522
--    (0.33 secs, 109,402,440 bytes)
--    λ> length (show (sylvester5 23))
--    1707522
--    (0.30 secs, 108,676,256 bytes)

graficaSylvester :: Integer -> Integer -> IO ()
graficaSylvester d n =
  plotList [ Key Nothing
           , PNG ("La_sucesion_de_Sylvester_" ++ show (d,n) ++ ".png")
           ]
           [sylvester5 k `mod` (10^d) | k <- [0..n]]

import Data.List (genericIndex)

import Data.Array ((!), array)

import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG))

-- 1ª solución (por recursión)

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

sylvester1 :: Integer -> Integer

sylvester1 0 = 2

sylvester1 n = 1 + product [sylvester1 k | k <- [0..n-1]]

-- 2ª solución (con programación dinámica)

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

sylvester2 :: Integer -> Integer

sylvester2 n = v ! n where

v = array (0,n) [(i,f i) | i <- [0..n]]

f 0 = 2

f m = 1 + product [v!k | k <- [0..m-1]]

-- 3ª solución

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

-- Observando que

-- S(n) = 1 + S(0)*S(1)*...*S(n-2)*S(n-1)

-- = 1 + (1 + S(0)*S(1)*...*S(n-2))*S(n-1) - S(n-1)

-- = 1 + S(n-1)*S(n-1) - S(n-1)

-- = 1 + S(n-1)^2 - S(n-1)

-- se obtiene la siguiente definición.

sylvester3 :: Integer -> Integer

sylvester3 0 = 2

sylvester3 n = 1 + x^2 - x

where x = sylvester3 (n-1)

-- 4ª solución

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

sylvester4 :: Integer -> Integer

sylvester4 n = v ! n where

v = array (0,n) [(i,f i) | i <- [0..n]]

f 0 = 2

f m = 1 + x^2 - x

where x = v ! (m-1)

-- 5ª solución

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

sylvester5 :: Integer -> Integer

sylvester5 n = sucSylvester5 `genericIndex` n

sucSylvester5 :: [Integer]

sucSylvester5 = iterate (\x -> (x-1)*x+1) 2

-- La comparación es

-- λ> length (show (sylvester1 23))

-- 1707522

-- (6.03 secs, 4,090,415,704 bytes)

-- λ> length (show (sylvester2 23))

-- 1707522

-- (0.33 secs, 109,477,296 bytes)

-- λ> length (show (sylvester3 23))

-- 1707522

-- (0.35 secs, 109,395,136 bytes)

-- λ> length (show (sylvester4 23))

-- 1707522

-- (0.33 secs, 109,402,440 bytes)

-- λ> length (show (sylvester5 23))

-- 1707522

-- (0.30 secs, 108,676,256 bytes)

graficaSylvester :: Integer -> Integer -> IO ()

graficaSylvester d n =

plotList [ Key Nothing

, PNG ("La_sucesion_de_Sylvester_" ++ show (d,n) ++ ".png")

]

[sylvester5 k `mod` (10^d) | k <- [0..n]]

Cálculo de pi mediante la variante de Euler de la serie armónica

PorJosé A. Alonso 9-abril-202026-marzo-2022

En el artículo El desarrollo más bello de Pi como suma infinita, Miguel Ángel Morales comenta el desarrollo de pi publicado por Leonhard Euler en su libro «Introductio in Analysis Infinitorum» (1748).

El desarrollo es el siguiente

y se obtiene a partir de la serie armónica

modificando sólo el signo de algunos términos según el siguiente criterio:

Dejamos un + cuando el denominador de la fracción sea un 2 o un primo de la forma 4m-1.
Cambiamos a – si el denominador de la fracción es un primo de la forma 4m+1.
Si el número es compuesto ponemos el signo que quede al multiplicar los signos correspondientes a cada factor.

Por ejemplo,

la de denominador 3 = 4×1-1 lleva un +,
la de denominador 5 = 4×1+1 lleva un -,
la de denominador 13 = 4×3+1 lleva un -,
la de denominador 6 = 2×3 lleva un + (porque los dos llevan un +),
la de denominador 10 = 2×5 lleva un – (porque el 2 lleva un + y el 5 lleva un -) y
la de denominador 50 = 5x5x2 lleva un + (un – por el primer 5, otro – por el segundo 5 y un + por el 2).

Definir las funciones

  aproximacionPi :: Int -> Double
  grafica        :: Int -> IO ()

1 2	aproximacionPi :: Int -> Double grafica :: Int -> IO ()

tales que

(aproximacionPi n) es la aproximación de pi obtenida sumando los n primeros términos de la serie de Euler. Por ejemplo.

     aproximacionPi 1        ==  1.0
     aproximacionPi 10       ==  2.3289682539682537
     aproximacionPi 100      ==  2.934318000847734
     aproximacionPi 1000     ==  3.0603246224585128
     aproximacionPi 10000    ==  3.1105295744825403
     aproximacionPi 100000   ==  3.134308801935256
     aproximacionPi 1000000  ==  3.1395057903490806

aproximacionPi 1 == 1.0

aproximacionPi 10 == 2.3289682539682537

aproximacionPi 100 == 2.934318000847734

aproximacionPi 1000 == 3.0603246224585128

aproximacionPi 10000 == 3.1105295744825403

aproximacionPi 100000 == 3.134308801935256

aproximacionPi 1000000 == 3.1395057903490806

(grafica n) dibuja la gráfica de las aproximaciones de pi usando k sumando donde k toma los valores de la lista [100,110..n]. Por ejemplo, al evaluar (grafica 4000) se obtiene

Soluciones

import Data.Numbers.Primes
import Graphics.Gnuplot.Simple

-- 1ª definición
-- =============

aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [1 / fromIntegral (k * signo k) | k <- [1..n]] 

signoPrimo :: Int -> Int
signoPrimo 2 = 1
signoPrimo p | p `mod` 4 == 3 = 1
             | otherwise      = -1

signo :: Int -> Int
signo n | isPrime n = signoPrimo n
        | otherwise = product (map signoPrimo (primeFactors n))

-- 2ª definición
-- =============

aproximacionPi2 :: Int -> Double
aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]
serieEuler =
  scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]

-- Definición de grafica
-- =====================

grafica :: Int -> IO ()
grafica n = 
    plotList [Key Nothing]
             [(k,aproximacionPi2 k) | k <- [100,110..n]]

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- 1ª definición

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

aproximacionPi :: Int -> Double

aproximacionPi n =

sum [1 / fromIntegral (k * signo k) | k <- [1..n]]

signoPrimo :: Int -> Int

signoPrimo 2 = 1

signoPrimo p | p `mod` 4 == 3 = 1

| otherwise = -1

signo :: Int -> Int

signo n | isPrime n = signoPrimo n

| otherwise = product (map signoPrimo (primeFactors n))

-- 2ª definición

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

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]

serieEuler =

scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]

-- Definición de grafica

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

grafica :: Int -> IO ()

grafica n =

plotList [Key Nothing]

[(k,aproximacionPi2 k) | k <- [100,110..n]]

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

Productos de dos y tres números consecutivos

PorJosé A. Alonso 3-abril-202010-abril-2020

Definir la función

   productos :: Integer -> Integer -> [[Integer]]

1	productos :: Integer -> Integer -> [[Integer]]

tal que (productos n x) es las listas de n elementos consecutivos cuyo producto es x. Por ejemplo,

   productos 2 6     ==  [[2,3]]
   productos 3 6     ==  [[1,2,3]]
   productos 4 1680  ==  [[5,6,7,8]]
   productos 2 5     ==  []

productos 2 6 == [[2,3]]

productos 3 6 == [[1,2,3]]

productos 4 1680 == [[5,6,7,8]]

productos 2 5 == []

Comprobar con QuickCheck que si n > 0 y x > 0, entonces

   productos n (product [x..x+n-1]) == [[x..x+n-1]]

1	productos n (product [x..x+n-1]) == [[x..x+n-1]]

Usando productos, definir la función

   productosDe2y3consecutivos :: [Integer]

1	productosDe2y3consecutivos :: [Integer]

cuyos elementos son los números naturales (no nulos) que pueden expresarse simultáneamente como producto de dos y tres números consecutivos. Por ejemplo,

   head productosDe2y3consecutivos  ==  6

1	head productosDe2y3consecutivos == 6

Nota. Según demostró Mordell en 1962, productosDe2y3consecutivos sólo tiene dos elementos.

Soluciones

import Test.QuickCheck

-- 1ª definición
productos1 :: Integer -> Integer -> [[Integer]]
productos1 n x =
  [[y..y+n-1] | y <- [1..x]
              , product [y..y+n-1] == x]

-- 2ª definición
productos2 :: Integer -> Integer -> [[Integer]]
productos2 n x =
  [[z..z+n-1] | z <- [1..y]
              , product [z..z+n-1] == x]
  where y = head (filter (\y -> y^n >= x) [2..])

productos :: Integer -> Integer -> [[Integer]]
productos = productos2

prop_productos n x =
  n > 0 && x > 0 ==> productos n (product [x..x+n-1]) == [[x..x+n-1]]

-- La comprobación es
--    λ> quickCheck prop_productos
--    +++ OK, passed 100 tests.
--    (0.10 secs, 26409644 bytes)

productosDe2y3consecutivos :: [Integer]
productosDe2y3consecutivos =
  [x | x <- [1..] 
     , (not . null) (productos 2 x)
     , (not . null) (productos 3 x)]

-- El cálculo es
--    λ> take 2 productosDe2y3consecutivos
--    [6,210]
--    λ> productos 2 210
--    [[14,15]]
--    λ> productos 3 210
--    [[5,6,7]]

import Test.QuickCheck

-- 1ª definición

productos1 :: Integer -> Integer -> [[Integer]]

productos1 n x =

[[y..y+n-1] | y <- [1..x]

, product [y..y+n-1] == x]

-- 2ª definición

productos2 :: Integer -> Integer -> [[Integer]]

productos2 n x =

[[z..z+n-1] | z <- [1..y]

, product [z..z+n-1] == x]

where y = head (filter (\y -> y^n >= x) [2..])

productos :: Integer -> Integer -> [[Integer]]

productos = productos2

prop_productos n x =

n > 0 && x > 0 ==> productos n (product [x..x+n-1]) == [[x..x+n-1]]

-- La comprobación es

-- λ> quickCheck prop_productos

-- +++ OK, passed 100 tests.

-- (0.10 secs, 26409644 bytes)

productosDe2y3consecutivos :: [Integer]

productosDe2y3consecutivos =

[x | x <- [1..]

, (not . null) (productos 2 x)

, (not . null) (productos 3 x)]

-- El cálculo es

-- λ> take 2 productosDe2y3consecutivos

-- [6,210]

-- λ> productos 2 210

-- [[14,15]]

-- λ> productos 3 210

-- [[5,6,7]]

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

Pensamiento

«El verdadero viaje de descubrimiento no consiste en buscar nuevos paisajes sino en tener nuevos ojos.»

Marcel Proust.

Números de Bell

PorJosé A. Alonso 5-febrero-202012-febrero-2020

Una partición de un conjunto A es un conjunto de subconjuntos no vacíos de A, disjuntos dos a dos y cuya unión es A. Por ejemplo, el conjunto {1, 2, 3} tiene exactamente 5 particiones:

   {{1}, {2}, {3}}
   {{1,2}, {3}}
   {{1,3}, {2}}
   {{1}, {2,3}}
   {{1,2,3}}

{{1}, {2}, {3}}

{{1,2}, {3}}

{{1,3}, {2}}

{{1}, {2,3}}

El n-ésimo número de Bell, B(n), es el número de particiones de un conjunto de n elementos. Por lo visto anteriormentem B(3) = 5.

Definir las funciones

   particiones :: [a] -> [[[a]]]
   bell :: Integer -> Integer

1 2	particiones :: [a] -> [[[a]]] bell :: Integer -> Integer

tales que

(particiones xs) es el conjunto de las particiones de xs. Por ejemplo,

     λ> particiones [1,2]
     [[[1,2]],[[1],[2]]]
     λ> particiones [1,2,3]
     [[[1,2,3]],[[1],[2,3]],[[1,2],[3]],[[2],[1,3]],[[1],[2],[3]]]
     λ> particiones "abcd"
     [["abcd"],["a","bcd"],["ab","cd"],["b","acd"],["abc","d"],["bc","ad"],
      ["ac","bd"],["c","abd"],["a","b","cd"],["a","bc","d"],["a","c","bd"],
      ["ab","c","d"],["b","ac","d"],["b","c","ad"],["a","b","c","d"]]

λ> particiones [1,2]

[[[1,2]],[[1],[2]]]

λ> particiones [1,2,3]

[[[1,2,3]],[[1],[2,3]],[[1,2],[3]],[[2],[1,3]],[[1],[2],[3]]]

λ> particiones "abcd"

[["abcd"],["a","bcd"],["ab","cd"],["b","acd"],["abc","d"],["bc","ad"],

["ac","bd"],["c","abd"],["a","b","cd"],["a","bc","d"],["a","c","bd"],

["ab","c","d"],["b","ac","d"],["b","c","ad"],["a","b","c","d"]]

(bell n) es el n-ésimo número de Bell. Por ejemplo,

     λ> bell 3
     5
     λ> map bell [0..10]
     [1,1,2,5,15,52,203,877,4140,21147,115975]

λ> bell 3

λ> map bell [0..10]

[1,1,2,5,15,52,203,877,4140,21147,115975]

Comprobar con QuickCheck que (bell n) es equivalente a la función B(n) definida por

$B(0) = 1$
$B(n) = \displaystyle \sum_{k=0}^{n-1} \binom{n-1}{k} B(k)$

Soluciones

import Data.List (genericLength)
import Test.QuickCheck

-- Definición de particiones
-- =========================

particiones :: [a] -> [[[a]]]
particiones [] = [[]]
particiones (x:xs) =
  concat [([x] : yss) : inserta x yss | yss <- ysss]
  where ysss = particiones xs

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los
-- elementos de yss. Por ejemplo, 
--    λ> inserta 1 [[2,3],[4],[5,6,7]]
--    [[[1,2,3],[4],[5,6,7]],[[2,3],[1,4],[5,6,7]],[[2,3],[4],[1,5,6,7]]]
inserta :: a -> [[a]] -> [[[a]]]
inserta _ []       = []
inserta x (ys:yss) = ((x:ys):yss) : [ys : zs | zs <- inserta x yss] 

-- Definición de Bell
-- ==================

bell :: Integer -> Integer
bell n = genericLength (particiones [1..n])

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

prop_Bell :: Integer -> Property
prop_Bell n =
  n >= 0 ==> bell n == b n

b :: Integer -> Integer
b 0 = 1
b n = sum [comb (n-1) k * b k | k <- [0..n-1]]

comb :: Integer -> Integer -> Integer
comb n k = product [n-k+1..n] `div` product [1..k]

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

import Data.List (genericLength)

import Test.QuickCheck

-- Definición de particiones

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

particiones :: [a] -> [[[a]]]

particiones [] = [[]]

particiones (x:xs) =

concat [([x] : yss) : inserta x yss | yss <- ysss]

where ysss = particiones xs

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los

-- elementos de yss. Por ejemplo,

-- λ> inserta 1 [[2,3],[4],[5,6,7]]

-- [[[1,2,3],[4],[5,6,7]],[[2,3],[1,4],[5,6,7]],[[2,3],[4],[1,5,6,7]]]

inserta :: a -> [[a]] -> [[[a]]]

inserta _ [] = []

inserta x (ys:yss) = ((x:ys):yss) : [ys : zs | zs <- inserta x yss]

-- Definición de Bell

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

bell :: Integer -> Integer

bell n = genericLength (particiones [1..n])

-- Propiedad

-- =========

prop_Bell :: Integer -> Property

prop_Bell n =

n >= 0 ==> bell n == b n

b :: Integer -> Integer

b 0 = 1

b n = sum [comb (n-1) k * b k | k <- [0..n-1]]

comb :: Integer -> Integer -> Integer

comb n k = product [n-k+1..n] `div` product [1..k]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_Bell

-- +++ OK, passed 100 tests.

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

Pensamiento

«Cambiemos nuestra actitud tradicional en la construcción de programas. En lugar de imaginar que nuestra tarea principal es indicarle a una computadora lo que debe hacer, concentrémonos más bien en explicarle a los seres humanos lo que queremos que haga una computadora.»

Donald Knuth.

Acotación del primorial

PorJosé A. Alonso 31-enero-20207-febrero-2020

El primorial de un número natural n es el producto de todos los números primos menores o iguales a n. Por ejemplo, el primorial de 5 es 30 porque el producto de los primos menores o iguales que 5 es

   2 * 3 * 5 = 30

1	2 * 3 * 5 = 30

La propiedad de Erdös de acotación de los primoriales afirma que

Para todo número natural n, su primorial es menor o igual que 4ⁿ.

Definir las funciones

   primorial :: Integer -> Integer
   primoriales :: [Integer]

1 2	primorial :: Integer -> Integer primoriales :: [Integer]

tales que

(primorial n) es el primorial de n. Por ejemplo,

     primorial 3  ==  6
     primorial 5  ==  30
     primorial 8  ==  210

primorial 3 == 6

primorial 5 == 30

primorial 8 == 210

primoriales es la sucesión de los primoriales. Por ejemplo,

   λ> take 15 primoriales
   [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

1 2	λ> take 15 primoriales [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

Comprobar con QuickCheck la propiedad de Erdös de acotación de los primoriales.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

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

primorial :: Integer -> Integer
primorial n = product (takeWhile (<= n) primes)

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

primorial2 :: Integer -> Integer
primorial2 0 = 1
primorial2 n | gcd n x == 1 = n*x
             | otherwise    = x
  where x = primorial2 (n-1)

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

--    λ> length (show (primorial (5*10^5)))
--    216852
--    (1.65 secs, 2,472,977,584 bytes)
--    λ> length (show (primorial2 (5*10^5)))
--    216852
--    (3.56 secs, 2,719,162,272 bytes)

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

--    λ> take 15 primoriales
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales :: [Integer]
primoriales = map primorial [0..]

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

--    λ> take 15 primoriales2
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales2 :: [Integer]
primoriales2 = map primorial2 [0..]

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

--    λ> take 15 primoriales3
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales3 :: [Integer]
primoriales3 = scanl1 f [1..]
  where f x n | gcd n x == 1 = n*x
              | otherwise    = x

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

--    λ> minimum (take 5000 primoriales)
--    1
--    (1.56 secs, 4,857,760,464 bytes)
--    λ> minimum (take 5000 primoriales2)
--    1
--    (9.39 secs, 10,942,848,240 bytes)
--    λ> minimum (take 5000 primoriales3)
--    1
--    (0.01 secs, 5,575,024 bytes)
--    
--    λ> minimum (take 6000 primoriales)
--    1
--    (2.22 secs, 7,013,937,248 bytes)
--    λ> minimum (take 6000 primoriales3)
--    1
--    (0.01 secs, 6,737,328 bytes)

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

prop_primorial :: Integer -> Property
prop_primorial n =
  n >= 0 ==> primorial n <= 4^n

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

import Data.Numbers.Primes

import Test.QuickCheck

-- 1ª definición de primorial

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

primorial :: Integer -> Integer

primorial n = product (takeWhile (<= n) primes)

-- 2ª definición de primorial

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

primorial2 :: Integer -> Integer

primorial2 0 = 1

primorial2 n | gcd n x == 1 = n*x

| otherwise = x

where x = primorial2 (n-1)

-- Comparación de eficiencia

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

-- λ> length (show (primorial (5*10^5)))

-- 216852

-- (1.65 secs, 2,472,977,584 bytes)

-- λ> length (show (primorial2 (5*10^5)))

-- 216852

-- (3.56 secs, 2,719,162,272 bytes)

-- 1ª definición de primoriales

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

-- λ> take 15 primoriales

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales :: [Integer]

primoriales = map primorial [0..]

-- 2ª definición de primoriales

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

-- λ> take 15 primoriales2

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales2 :: [Integer]

primoriales2 = map primorial2 [0..]

-- 3ª definición de primoriales

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

-- λ> take 15 primoriales3

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales3 :: [Integer]

primoriales3 = scanl1 f [1..]

where f x n | gcd n x == 1 = n*x

| otherwise = x

-- Comparación de eficiencia

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

-- λ> minimum (take 5000 primoriales)

-- 1

-- (1.56 secs, 4,857,760,464 bytes)

-- λ> minimum (take 5000 primoriales2)

-- 1

-- (9.39 secs, 10,942,848,240 bytes)

-- λ> minimum (take 5000 primoriales3)

-- 1

-- (0.01 secs, 5,575,024 bytes)

-- λ> minimum (take 6000 primoriales)

-- 1

-- (2.22 secs, 7,013,937,248 bytes)

-- λ> minimum (take 6000 primoriales3)

-- 1

-- (0.01 secs, 6,737,328 bytes)

-- Propiedad

-- =========

prop_primorial :: Integer -> Property

prop_primorial n =

n >= 0 ==> primorial n <= 4^n

-- La comprobación es

-- λ> quickCheck prop_primorial

-- +++ OK, passed 100 tests.

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

Pensamiento

«Las matemáticas son la reina de las ciencias y la teoría de los números es la reina de las matemáticas.»

Carl Friedrich Gauss.

Productos de sumas de cuatro cuadrados

PorJosé A. Alonso 24-enero-202031-enero-2020

Definir la función

   productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

1	productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

tal que (productoSuma4Cuadrados as bs cs ds) es el producto de las sumas de los cuadrados de cada una de las listas que ocupan la misma posición (hasta que alguna se acaba). Por ejemplo,

   productoSuma4Cuadrados [2,3] [1,5] [4,6] [0,3,9]
   = (2² + 1² + 4² + 0²) * (3² + 5² + 6² + 3²)
   = (4 +  1 + 16  + 0)  * (9 + 25 + 36  + 9)
   = 1659

productoSuma4Cuadrados [2,3] [1,5] [4,6] [0,3,9]

= (2² + 1² + 4² + 0²) * (3² + 5² + 6² + 3²)

= (4 + 1 + 16 + 0) * (9 + 25 + 36 + 9)

= 1659

Comprobar con QuickCheckWith que si as, bs cs y ds son listas no vacías de enteros positivos, entonces (productoSuma4Cuadrados as bs cs ds) se puede escribir como la suma de los cuadrados de cuatro enteros positivos.

Soluciones

import Data.List (zip4)
import Test.QuickCheck

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

productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a
productoSuma4Cuadrados (a:as) (b:bs) (c:cs) (d:ds) =
  (a^2+b^2+c^2+d^2) * productoSuma4Cuadrados as bs cs ds
productoSuma4Cuadrados _ _ _ _ = 1

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

productoSuma4Cuadrados2 :: Integral a => [a] -> [a] -> [a] -> [a] -> a
productoSuma4Cuadrados2 as bs cs ds =
  product [a^2 + b^2 + c^2 + d^2 | (a,b,c,d) <- zip4 as bs cs ds]

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

-- La propiedad es
prop_productoSuma4Cuadrados ::
  [Integer] -> [Integer] -> [Integer] -> [Integer] -> Property
prop_productoSuma4Cuadrados as bs cs ds =
  all (not . null) [as, bs, cs, ds]
  ==> 
  esSuma4Cuadrados (productoSuma4Cuadrados as' bs' cs' ds')
  where as' = [1 + abs a | a <- as]
        bs' = [1 + abs b | b <- bs]
        cs' = [1 + abs c | c <- cs]
        ds' = [1 + abs d | d <- ds]

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por
-- ejemplo, 
--    esSuma4Cuadrados 42  ==  True
--    esSuma4Cuadrados 11  ==  False
--    esSuma4Cuadrados 41  ==  False
esSuma4Cuadrados :: Integer -> Bool
esSuma4Cuadrados = not . null . sumas4Cuadrados

-- (sumas4Cuadrados n) es la lista de las descomposiciones de n como
-- sumas de 4 cuadrados. Por ejemplo,
--    sumas4Cuadrados 42  ==  [(16,16,9,1),(25,9,4,4),(36,4,1,1)]
sumas4Cuadrados :: Integer -> [(Integer,Integer,Integer,Integer)]
sumas4Cuadrados n =
  [(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]
                   , b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]
                   , c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]
                   , let d = n - a^2 - b^2 - c^2
                   , c^2 <= d 
                   , esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))

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

import Data.List (zip4)

import Test.QuickCheck

-- 1ª solución

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

productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

productoSuma4Cuadrados (a:as) (b:bs) (c:cs) (d:ds) =

(a^2+b^2+c^2+d^2) * productoSuma4Cuadrados as bs cs ds

productoSuma4Cuadrados _ _ _ _ = 1

-- 2ª solución

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

productoSuma4Cuadrados2 :: Integral a => [a] -> [a] -> [a] -> [a] -> a

productoSuma4Cuadrados2 as bs cs ds =

product [a^2 + b^2 + c^2 + d^2 | (a,b,c,d) <- zip4 as bs cs ds]

-- Propiedad

-- =========

-- La propiedad es

prop_productoSuma4Cuadrados ::

[Integer] -> [Integer] -> [Integer] -> [Integer] -> Property

prop_productoSuma4Cuadrados as bs cs ds =

all (not . null) [as, bs, cs, ds]

==>

esSuma4Cuadrados (productoSuma4Cuadrados as' bs' cs' ds')

where as' = [1 + abs a | a <- as]

bs' = [1 + abs b | b <- bs]

cs' = [1 + abs c | c <- cs]

ds' = [1 + abs d | d <- ds]

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por

-- ejemplo,

-- esSuma4Cuadrados 42 == True

-- esSuma4Cuadrados 11 == False

-- esSuma4Cuadrados 41 == False

esSuma4Cuadrados :: Integer -> Bool

esSuma4Cuadrados = not . null . sumas4Cuadrados

-- (sumas4Cuadrados n) es la lista de las descomposiciones de n como

-- sumas de 4 cuadrados. Por ejemplo,

-- sumas4Cuadrados 42 == [(16,16,9,1),(25,9,4,4),(36,4,1,1)]

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

sumas4Cuadrados n =

[(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]

, b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]

, c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]

, let d = n - a^2 - b^2 - c^2

, c^2 <= d

, esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- La comprobación es

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

-- +++ OK, passed 100 tests.

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

Pensamiento

¿Vivir? Sencillamente:
la sed y el agua cerca …
o el agua lejos, más, la sed y el agua,
un poco de cansancio ¡y a beberla!.

Antonio Machado

Teorema de Liouville sobre listas CuCu

PorJosé A. Alonso 3-enero-202010-enero-2020

Una lista CuCu es una lista de números enteros positivos tales que la suma de sus Cubos es igual al Cuadrado de su suma. Por ejemplo, [1, 2, 3, 2, 4, 6] es una lista CuCu ya que

   1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

1	1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

La lista de Liouville correspondiente al número entero positivo n es la lista formada por el número de divisores de cada divisor de n. Por ejemplo, para el número 20 se tiene que sus divisores son

   1, 2, 4, 5, 10, 20

1	1, 2, 4, 5, 10, 20

puesto que el número de sus divisores es

El 1 tiene 1 divisor (el 1 solamente).
El 2 tiene 2 divisores (el 1 y el 2).
El 4 tiene 3 divisores (el 1, el 2 y el 4).
El 5 tiene 2 divisores (el 1 y el 5).
El 10 tiene 4 divisores (el 1, el 2, el 5 y el 10).
El 20 tiene 6 divisores (el 1, el 2, el 4, el 5, el 10 y el 20).

la lista de Liouville de 20 es [1, 2, 3, 2, 4, 6] que, como se comentó anteriormente, es una lista CuCu.

El teorema de Lioville afirma que todas las lista de Lioville son CuCu.

Definir las funciones

   esCuCu :: [Integer] -> Bool
   liouville :: Integer -> [Integer]

1 2	esCuCu :: [Integer] -> Bool liouville :: Integer -> [Integer]

tales que

(esCuCu xs) se verifica si la lista xs es CuCu; es decir, la suma de los cubos de sus elementos es igual al cuadrado de su suma. Por ejemplo,

     esCuCu [1,2,3]        ==  True
     esCuCu [1,2,3,2]      ==  False
     esCuCu [1,2,3,2,4,6]  ==  True

esCuCu [1,2,3] == True

esCuCu [1,2,3,2] == False

esCuCu [1,2,3,2,4,6] == True

(liouville n) es la lista de Lioville correspondiente al número n. Por ejemplo,

     liouville 20  ==  [1,2,3,2,4,6]
     liouville 60  ==  [1,2,2,3,2,4,4,6,4,6,8,12]
     length (liouville (product [1..25]))  ==  340032

liouville 20 == [1,2,3,2,4,6]

liouville 60 == [1,2,2,3,2,4,4,6,4,6,8,12]

length (liouville (product [1..25])) == 340032

Comprobar con QuickCheck

que para todo entero positivo n, (liouville (2^n)) es la lista [1,2,3,…,n+1] y
el teorema de Lioville; es decir, para todo entero positivo n, (liouville n) es una lista CuCu.

Nota: Este ejercicio está basado en Cómo generar conjuntos CuCu de Gaussianos.

Soluciones

import Data.List (genericLength, group, inits, sort)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

esCuCu :: [Integer] -> Bool
esCuCu xs = sum (map (^3) xs) == (sum xs)^2

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

liouville :: Integer -> [Integer]
liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. Por ejemplo, 
--   divisores 30  ==  [1,2,3,5,6,10,15,30]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n], n `mod` x == 0]

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

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

liouville2 :: Integer -> [Integer]
liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones
-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj
-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo, 
--   divisores2 30  ==  [1,2,3,5,6,10,15,30]
divisores2 :: Integer -> [Integer]
divisores2 = sort
           . map (product . concat)
           . sequence
           . map inits
           . group
           . primeFactors

-- (numeroDivisores2 x) es el número de divisores de x. Por ejemplo, 
--    numeroDivisores2 12  ==  6
--    numeroDivisores2 25  ==  3
numeroDivisores2 :: Integer -> Integer
numeroDivisores2 =
  product . map ((+1) . genericLength) . group . primeFactors

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

-- La comparación es
--    λ> length (liouville (product [1..11]))
--    540
--    (13.66 secs, 7,983,550,640 bytes)
--    λ> length (liouville2 (product [1..11]))
--    540
--    (0.01 secs, 1,255,328 bytes)

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

-- La propiedad es
prop_Liouville :: Integer -> Property
prop_Liouville n =
  n > 0 ==> liouville2 (2^n) == [1..n+1]

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

-- Teorema de Liouville
-- ====================

-- La propiedad es
teorema_Liouville :: Integer -> Property
teorema_Liouville n =
  n > 0 ==> esCuCu (liouville n)

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

import Data.List (genericLength, group, inits, sort)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

esCuCu :: [Integer] -> Bool

esCuCu xs = sum (map (^3) xs) == (sum xs)^2

-- 1ª definición de liouville

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

liouville :: Integer -> [Integer]

liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. Por ejemplo,

-- divisores 30 == [1,2,3,5,6,10,15,30]

divisores :: Integer -> [Integer]

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

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

-- numeroDivisores 12 == 6

-- numeroDivisores 25 == 3

numeroDivisores :: Integer -> Integer

numeroDivisores n = genericLength (divisores n)

-- 2ª definición de liouville

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

liouville2 :: Integer -> [Integer]

liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones

-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj

-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo,

-- divisores2 30 == [1,2,3,5,6,10,15,30]

divisores2 :: Integer -> [Integer]

divisores2 = sort

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

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

-- numeroDivisores2 12 == 6

-- numeroDivisores2 25 == 3

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

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

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (liouville (product [1..11]))

-- 540

-- (13.66 secs, 7,983,550,640 bytes)

-- λ> length (liouville2 (product [1..11]))

-- 540

-- (0.01 secs, 1,255,328 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_Liouville :: Integer -> Property

prop_Liouville n =

n > 0 ==> liouville2 (2^n) == [1..n+1]

-- La comprobación es

-- λ> quickCheck prop_Liouville

-- +++ OK, passed 100 tests.

-- Teorema de Liouville

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

-- La propiedad es

teorema_Liouville :: Integer -> Property

teorema_Liouville n =

n > 0 ==> esCuCu (liouville n)

-- La comprobación es

-- λ> quickCheck teorema_Liouville

-- +++ OK, passed 100 tests.

Pensamiento

¡Oh, tarde viva y quieta
que opuso al panta rhei su nada corre.

Antonio Machado

Cálculo de pi usando la fórmula de Vieta

PorJosé A. Alonso 16-diciembre-201923-diciembre-2019

La fórmula de Vieta para el cálculo de pi es la siguiente

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

1 2	aproximacionPi :: Int -> Double errorPi :: Double -> Int

tales que

(aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,

     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795

aproximacionPi 5 == 3.140331156954753

aproximacionPi 10 == 3.1415914215112

aproximacionPi 15 == 3.141592652386592

aproximacionPi 20 == 3.1415926535886207

aproximacionPi 25 == 3.141592653589795

(errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

errorPi 0.1 == 2

errorPi 0.01 == 4

errorPi 0.001 == 6

errorPi 0.0001 == 7

errorPi 1e-4 == 7

errorPi 1e-14 == 24

pi == 3.141592653589793

aproximacionPi 24 == 3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

-- 1ª definición de aproximacionPi

aproximacionPi :: Int -> Double

aproximacionPi n = product [2 / aux x | x <- [0..n]]

where

aux 0 = 1

aux 1 = sqrt 2

aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = product [2/x | x <- 1 : xs]

where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi

aproximacionPi3 :: Int -> Double

aproximacionPi3 n = product (2 : take n (map (2/) xs))

where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi

errorPi :: Double -> Int

errorPi x = head [n | n <- [1..]

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x = until aceptable (+1) 1

where aceptable n = abs (pi - aproximacionPi n) < x

Pensamiento

El tiempo que la barba me platea,
cavó mis ojos y agrandó mi frente,
va siendo en mi recuerdo transparente,
y mientras más al fondo, más clarea.

Antonio Machado

Pares definidos por su MCD y su MCM

PorJosé A. Alonso 13-diciembre-201920-diciembre-2019

Definir las siguientes funciones

   pares  :: Integer -> Integer -> [(Integer,Integer)]
   nPares :: Integer -> Integer -> Integer

1 2	pares :: Integer -> Integer -> [(Integer,Integer)] nPares :: Integer -> Integer -> Integer

tales que

(pares a b) es la lista de los pares de números enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     pares 3 3  == [(3,3)]
     pares 4 12 == [(4,12),(12,4)]
     pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]
     pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
     pares 2 7  == []
     pares 12 3  ==  []
     length (pares 3 (product [3,5..91]))  ==  8388608

pares 3 3 == [(3,3)]

pares 4 12 == [(4,12),(12,4)]

pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]

pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

pares 2 7 == []

pares 12 3 == []

length (pares 3 (product [3,5..91])) == 8388608

(nPares a b) es el número de pares de enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     nPares 3 3   ==  1
     nPares 4 12  ==  2
     nPares 2 12  ==  4
     nPares 2 60  ==  8
     nPares 2 7   ==  0
     nPares 12 3  ==  0
     nPares 3 (product [3..3*10^4]) `mod` (10^12)  ==  477999992832
     length (show (nPares 3 (product [3..3*10^4])))  ==  977

nPares 3 3 == 1

nPares 4 12 == 2

nPares 2 12 == 4

nPares 2 60 == 8

nPares 2 7 == 0

nPares 12 3 == 0

nPares 3 (product [3..3*10^4]) `mod` (10^12) == 477999992832

length (show (nPares 3 (product [3..3*10^4]))) == 977

Soluciones

import Data.Numbers.Primes (primeFactors)
import Data.List (genericLength, group, nub, sort, subsequences)
import Test.QuickCheck

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

pares1 :: Integer -> Integer -> [(Integer,Integer)]
pares1 a b = [(x,y) | x <- [1..b]
                    , y <- [1..b]
                    , gcd x y == a
                    , lcm x y == b]

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

pares2 :: Integer -> Integer -> [(Integer,Integer)]
pares2 a b = [(x,y) | x <- [a,a+a..b]
                    , y <- [a,a+a..b]
                    , gcd x y == a
                    , lcm x y == b]

-- Comparación de eficiencia
--    λ> length (pares1 3 (product [3,5..11]))
--    16
--    (95.12 secs, 86,534,165,528 bytes)
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)

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

pares3 :: Integer -> Integer -> [(Integer,Integer)]
pares3 a b = [(x,y) | x <- [a,a+a..b]
                    , c `rem` x == 0
                    , let y = c `div` x
                    , gcd x y == a
                    ]
  where c = a * b

-- Comparacioń de eficiencia
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)
--    λ> length (pares3 3 (product [3,5..11]))
--    16
--    (0.01 secs, 878,104 bytes)

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

-- Para la cuarta definición de pares se observa la relación con los
-- factores primos
--    λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]
--    [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]
--    λ> [primeFactors x | (x,y) <- pares1 2 12]
--    [[2],[2,2],[2,3],[2,2,3]]
--    λ> [primeFactors x | (x,y) <- pares1 2 60]
--    [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 6 60]
--    [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 2 24]
--    [[2],[2,3],[2,2,2],[2,2,2,3]]
-- Se observa que cada pares se obtiene de uno de los subconjuntos de los
-- divisores primos de b/a. Por ejemplo,
--    λ> (a,b) = (2,24)
--    λ> b `div` a
--    12
--    λ> primeFactors it
--    [2,2,3]
--    λ> group it
--    [[2,2],[3]]
--    λ> subsequences it
--    [[],[[2,2]],[[3]],[[2,2],[3]]]
--    λ> map concat it
--    [[],[2,2],[3],[2,2,3]]
--    λ> map product it
--    [1,4,3,12]
--    λ> [(a * x, b `div` x) | x <- it]
--    [(2,24),(8,6),(6,8),(24,2)]
-- A partir de la observación se construye la siguiente definición

pares4 :: Integer -> Integer -> [(Integer,Integer)]
pares4 a b
  | b `mod` a /= 0 = []
  | otherwise =
    [(a * x, b `div` x)
    | x <- map (product . concat)
               ((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,
-- pero no necesariamente en el mismo orden. Por ejemplo,
--    λ> pares3 2 60 
--    [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
--    λ> pares4 2 60 
--    [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]
--    λ> pares3 2 60 == sort (pares4 2 60)
--    True

-- Comparacioń de eficiencia
--    λ> length (pares3 3 (product [3,5..17]))
--    64
--    (4.44 secs, 2,389,486,440 bytes)
--    λ> length (pares4 3 (product [3,5..17]))
--    64
--    (0.00 secs, 177,704 bytes)
  
-- Propiedades de equivalencia de pares
-- ====================================

prop_pares :: Integer -> Integer -> Property
prop_pares a b =
  a > 0 && b > 0 ==>
  all (== pares1 a b)
      [sort (f a b) | f <- [ pares2
                           , pares3
                           , pares4
                           ]]
  
prop_pares2 :: Integer -> Integer -> Property
prop_pares2 a b =
  a > 0 && b > 0 ==>
  all (== pares1 a (a * b))
      [sort (f a (a * b)) | f <- [ pares2
                                 , pares3
                                 , pares4
                                 ]]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2
--    +++ OK, passed 100 tests.
             
-- 1ª definición de nPares
-- =======================

nPares1 :: Integer -> Integer -> Integer
nPares1 a b = genericLength (pares4 a b)

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

nPares2 :: Integer -> Integer -> Integer
nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia
--    λ> nPares1 3 (product [3,5..91])
--    8388608
--    (4.68 secs, 4,178,295,920 bytes)
--    λ> nPares2 3 (product [3,5..91])
--    8388608
--    (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares
-- ===================================

prop_nPares :: Integer -> Integer -> Property
prop_nPares a b =
  a > 0 && b > 0 ==>
  nPares1 a (a * b) == nPares2 a (a * b)
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares
--    +++ 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

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

import Data.Numbers.Primes (primeFactors)

import Data.List (genericLength, group, nub, sort, subsequences)

import Test.QuickCheck

-- 1ª definición de pares

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

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

pares1 a b = [(x,y) | x <- [1..b]

, y <- [1..b]

, gcd x y == a

, lcm x y == b]

-- 2ª definición de pares

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

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

pares2 a b = [(x,y) | x <- [a,a+a..b]

, y <- [a,a+a..b]

, gcd x y == a

, lcm x y == b]

-- Comparación de eficiencia

-- λ> length (pares1 3 (product [3,5..11]))

-- 16

-- (95.12 secs, 86,534,165,528 bytes)

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- 3ª definición de pares

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

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

pares3 a b = [(x,y) | x <- [a,a+a..b]

, c `rem` x == 0

, let y = c `div` x

, gcd x y == a

]

where c = a * b

-- Comparacioń de eficiencia

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- λ> length (pares3 3 (product [3,5..11]))

-- 16

-- (0.01 secs, 878,104 bytes)

-- 4ª definición de pares

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

-- Para la cuarta definición de pares se observa la relación con los

-- factores primos

-- λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]

-- [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]

-- λ> [primeFactors x | (x,y) <- pares1 2 12]

-- [[2],[2,2],[2,3],[2,2,3]]

-- λ> [primeFactors x | (x,y) <- pares1 2 60]

-- [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 6 60]

-- [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 2 24]

-- [[2],[2,3],[2,2,2],[2,2,2,3]]

-- Se observa que cada pares se obtiene de uno de los subconjuntos de los

-- divisores primos de b/a. Por ejemplo,

-- λ> (a,b) = (2,24)

-- λ> b `div` a

-- 12

-- λ> primeFactors it

-- [2,2,3]

-- λ> group it

-- [[2,2],[3]]

-- λ> subsequences it

-- [[],[[2,2]],[[3]],[[2,2],[3]]]

-- λ> map concat it

-- [[],[2,2],[3],[2,2,3]]

-- λ> map product it

-- [1,4,3,12]

-- λ> [(a * x, b `div` x) | x <- it]

-- [(2,24),(8,6),(6,8),(24,2)]

-- A partir de la observación se construye la siguiente definición

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

pares4 a b

| b `mod` a /= 0 = []

| otherwise =

[(a * x, b `div` x)

| x <- map (product . concat)

((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,

-- pero no necesariamente en el mismo orden. Por ejemplo,

-- λ> pares3 2 60

-- [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

-- λ> pares4 2 60

-- [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]

-- λ> pares3 2 60 == sort (pares4 2 60)

-- True

-- Comparacioń de eficiencia

-- λ> length (pares3 3 (product [3,5..17]))

-- 64

-- (4.44 secs, 2,389,486,440 bytes)

-- λ> length (pares4 3 (product [3,5..17]))

-- 64

-- (0.00 secs, 177,704 bytes)

-- Propiedades de equivalencia de pares

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

prop_pares :: Integer -> Integer -> Property

prop_pares a b =

a > 0 && b > 0 ==>

all (== pares1 a b)

[sort (f a b) | f <- [ pares2

, pares3

, pares4

]]

prop_pares2 :: Integer -> Integer -> Property

prop_pares2 a b =

a > 0 && b > 0 ==>

all (== pares1 a (a * b))

[sort (f a (a * b)) | f <- [ pares2

, pares3

, pares4

]]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2

-- +++ OK, passed 100 tests.

-- 1ª definición de nPares

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

nPares1 :: Integer -> Integer -> Integer

nPares1 a b = genericLength (pares4 a b)

-- 2ª definición de nPares

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

nPares2 :: Integer -> Integer -> Integer

nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia

-- λ> nPares1 3 (product [3,5..91])

-- 8388608

-- (4.68 secs, 4,178,295,920 bytes)

-- λ> nPares2 3 (product [3,5..91])

-- 8388608

-- (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares

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

prop_nPares :: Integer -> Integer -> Property

prop_nPares a b =

a > 0 && b > 0 ==>

nPares1 a (a * b) == nPares2 a (a * b)

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares

-- +++ OK, passed 100 tests.

Pensamiento

Largo es el camino de la enseñanza por medio de teorías; breve y eficaz por medio de ejemplos. ~ Séneca

Sublistas con producto dado

PorJosé A. Alonso 5-diciembre-201917-diciembre-2019

Definir las funciones

   sublistasConProducto :: Integer -> [Integer] -> [[Integer]]
   unifactorizables :: [Integer]

1 2	sublistasConProducto :: Integer -> [Integer] -> [[Integer]] unifactorizables :: [Integer]

tales que

(sublistasConProducto n xs) es la lista de las sublistas de la lista ordenada estrictamente creciente xs (cuyos elementos son enteros mayores que 1) cuyo producto es el número entero n (con n mayor que 1). Por ejemplo,

     λ> sublistasConProducto 72 [2,3,4,5,6,7,9,10,16]
     [[2,4,9],[3,4,6]]
     λ> sublistasConProducto 720 [2,3,4,5,6,7,9,10,16]
     [[2,3,4,5,6],[2,4,9,10],[3,4,6,10],[5,9,16]]
     λ> sublistasConProducto 2 [4,7]
     []
     λ> length (sublistasConProducto 1234567 [1..1234567])
     4

λ> sublistasConProducto 72 [2,3,4,5,6,7,9,10,16]

[[2,4,9],[3,4,6]]

λ> sublistasConProducto 720 [2,3,4,5,6,7,9,10,16]

[[2,3,4,5,6],[2,4,9,10],[3,4,6,10],[5,9,16]]

λ> sublistasConProducto 2 [4,7]

[]

λ> length (sublistasConProducto 1234567 [1..1234567])

unifactorizables es la lísta de los números enteros mayores que 1 que se pueden escribir sólo de una forma única como producto de enteros distintos mayores que uno. Por ejemplo,

     λ> take 20 unifactorizables
     [2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]
     λ> unifactorizables !! 300
     1873

λ> take 20 unifactorizables

[2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]

λ> unifactorizables !! 300

1873

Soluciones

import Test.QuickCheck
import Data.List (nub, sort, subsequences)

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

sublistasConProducto :: Integer -> [Integer] -> [[Integer]]
sublistasConProducto n xs =
  [ys | ys <- subsequences xs
      , product ys == n]

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

sublistasConProducto2 :: Integer -> [Integer] -> [[Integer]]
sublistasConProducto2 _ [] = []
sublistasConProducto2 n (x:xs)
  | x > n     = []
  | x == n    = [[x]]
  | r == 0    = map (x:) (sublistasConProducto2 q xs)
                ++ sublistasConProducto2 n xs
  | otherwise = sublistasConProducto2 n xs
  where (q,r) = quotRem n x

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

-- La propiedad es
prop_sublistasConProducto :: Integer -> [Integer] -> Bool
prop_sublistasConProducto n xs =
  sort (sublistasConProducto n' xs') == sublistasConProducto2 n' xs'
  where n'  = 2 + abs n
        xs' = (nub . sort . map ((+2) . abs)) xs

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

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

--    λ> sublistasConProducto 15 [1..23]
--    [[3,5],[1,3,5],[15],[1,15]]
--    (3.44 secs, 7,885,411,472 bytes)
--    λ> sublistasConProducto2 15 [1..23]
--    [[1,3,5],[1,15],[3,5],[15]]
--    (0.01 secs, 135,056 bytes)
--
--    λ> length (sublistasConProducto2 1234567 [1..1234567])
--    4
--    (1.49 secs, 1,054,380,480 bytes)

-- Definición de unifactorizables
-- ==============================

unifactorizables :: [Integer]
unifactorizables =
  [n | n <- [2..]
     , length (sublistasConProducto2 n [2..n]) == 1]

import Test.QuickCheck

import Data.List (nub, sort, subsequences)

-- 1ª solución

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

sublistasConProducto :: Integer -> [Integer] -> [[Integer]]

sublistasConProducto n xs =

[ys | ys <- subsequences xs

, product ys == n]

-- 2ª solución

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

sublistasConProducto2 :: Integer -> [Integer] -> [[Integer]]

sublistasConProducto2 _ [] = []

sublistasConProducto2 n (x:xs)

| x > n = []

| x == n = [[x]]

| r == 0 = map (x:) (sublistasConProducto2 q xs)

++ sublistasConProducto2 n xs

| otherwise = sublistasConProducto2 n xs

where (q,r) = quotRem n x

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sublistasConProducto :: Integer -> [Integer] -> Bool

prop_sublistasConProducto n xs =

sort (sublistasConProducto n' xs') == sublistasConProducto2 n' xs'

where n' = 2 + abs n

xs' = (nub . sort . map ((+2) . abs)) xs

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=30}) prop_sublistasConProducto

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> sublistasConProducto 15 [1..23]

-- [[3,5],[1,3,5],[15],[1,15]]

-- (3.44 secs, 7,885,411,472 bytes)

-- λ> sublistasConProducto2 15 [1..23]

-- [[1,3,5],[1,15],[3,5],[15]]

-- (0.01 secs, 135,056 bytes)

-- λ> length (sublistasConProducto2 1234567 [1..1234567])

-- 4

-- (1.49 secs, 1,054,380,480 bytes)

-- Definición de unifactorizables

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

unifactorizables :: [Integer]

unifactorizables =

[n | n <- [2..]

, length (sublistasConProducto2 n [2..n]) == 1]

Pensamiento

Y en el encinar,
¡luna redonda y beata,
siempre conmigo a la par!
Cerca de Úbeda la grande,
cuyos cerros nadie verá,
me iba siguiendo la luna
sobre el olivar.
Una luna jadeante,
siempre conmigo a la par.

Antonio Machado

Mayor producto de n dígitos consecutivos de un número

PorJosé A. Alonso 29-noviembre-20196-diciembre-2019

Definir la función

   mayorProducto :: Int -> Integer -> Integer

1	mayorProducto :: Int -> Integer -> Integer

tal que (mayorProducto n x) es el mayor producto de n dígitos consecutivos del número x (suponiendo que x tiene al menos n dígitos). Por ejemplo,

   mayorProducto 2 325                  ==  10
   mayorProducto 5 11111                ==  1
   mayorProducto 5 113111               ==  3
   mayorProducto 5 110111               ==  0
   mayorProducto 5 10151112             ==  10
   mayorProducto 5 101511124            ==  10
   mayorProducto 5 (product [1..1000])  ==  41472

mayorProducto 2 325 == 10

mayorProducto 5 11111 == 1

mayorProducto 5 113111 == 3

mayorProducto 5 110111 == 0

mayorProducto 5 10151112 == 10

mayorProducto 5 101511124 == 10

mayorProducto 5 (product [1..1000]) == 41472

Nota: Este ejercicio está basado en el problema 8 del Proyecto Euler

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)

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

mayorProducto :: Int -> Integer -> Integer
mayorProducto n x =
  maximum [product xs | xs <- segmentos n (digitos x)]

-- (digitos x) es la lista de las digitos del número x. Por ejemplo, 
--    digitos 325  ==  [3,2,5]
digitos :: Integer -> [Integer]
digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la
-- lista xs. Por ejemplo,
--    segmentos 2 [3,5,4,6]  ==  [[3,5],[5,4],[4,6]]
segmentos :: Int -> [Integer] -> [[Integer]]
segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

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

mayorProducto2 :: Int -> Integer -> Integer
mayorProducto2 n x = maximum (aux ns)
    where ns     = [read [d] | d <- show x]
          aux xs | length xs < n = []
                 | otherwise     = product (take n xs) : aux (tail xs)

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

mayorProducto3 :: Int -> Integer -> Integer
mayorProducto3 n = maximum
                 . map (product . take n)
                 . filter ((>=n) . length) 
                 . tails
                 . digitos

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

mayorProducto4 :: Int -> Integer -> Integer
mayorProducto4 n = maximum  
                 . map (product . map (fromIntegral . digitToInt)) 
                 . filter ((==n) . length) 
                 . concatMap inits
                 . tails 
                 . show

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

--    λ> mayorProducto 5 (product [1..500])
--    28224
--    (0.01 secs, 1,645,256 bytes)
--    λ> mayorProducto2 5 (product [1..500])
--    28224
--    (0.03 secs, 5,848,416 bytes)
--    λ> mayorProducto3 5 (product [1..500])
--    28224
--    (0.03 secs, 1,510,640 bytes)
--    λ> mayorProducto4 5 (product [1..500])
--    28224
--    (1.85 secs, 10,932,551,216 bytes)
--    
--    λ> mayorProducto 5 (product [1..7000])
--    46656
--    (0.10 secs, 68,590,808 bytes)
--    λ> mayorProducto2 5 (product [1..7000])
--    46656
--    (1.63 secs, 157,031,432 bytes)
--    λ> mayorProducto3 5 (product [1..7000])
--    46656
--    (1.55 secs, 65,727,176 bytes)

import Data.List (inits, tails)

import Data.Char (digitToInt)

-- 1ª solución

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

mayorProducto :: Int -> Integer -> Integer

mayorProducto n x =

maximum [product xs | xs <- segmentos n (digitos x)]

-- (digitos x) es la lista de las digitos del número x. Por ejemplo,

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

digitos :: Integer -> [Integer]

digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la

-- lista xs. Por ejemplo,

-- segmentos 2 [3,5,4,6] == [[3,5],[5,4],[4,6]]

segmentos :: Int -> [Integer] -> [[Integer]]

segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

-- 2ª solución

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

mayorProducto2 :: Int -> Integer -> Integer

mayorProducto2 n x = maximum (aux ns)

where ns = [read [d] | d <- show x]

aux xs | length xs < n = []

| otherwise = product (take n xs) : aux (tail xs)

-- 3ª solución

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

mayorProducto3 :: Int -> Integer -> Integer

mayorProducto3 n = maximum

. map (product . take n)

. filter ((>=n) . length)

. tails

. digitos

-- 4ª solución

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

mayorProducto4 :: Int -> Integer -> Integer

mayorProducto4 n = maximum

. map (product . map (fromIntegral . digitToInt))

. filter ((==n) . length)

. concatMap inits

. tails

. show

-- Comparación de eficiencia

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

-- λ> mayorProducto 5 (product [1..500])

-- 28224

-- (0.01 secs, 1,645,256 bytes)

-- λ> mayorProducto2 5 (product [1..500])

-- 28224

-- (0.03 secs, 5,848,416 bytes)

-- λ> mayorProducto3 5 (product [1..500])

-- 28224

-- (0.03 secs, 1,510,640 bytes)

-- λ> mayorProducto4 5 (product [1..500])

-- 28224

-- (1.85 secs, 10,932,551,216 bytes)

-- λ> mayorProducto 5 (product [1..7000])

-- 46656

-- (0.10 secs, 68,590,808 bytes)

-- λ> mayorProducto2 5 (product [1..7000])

-- 46656

-- (1.63 secs, 157,031,432 bytes)

-- λ> mayorProducto3 5 (product [1..7000])

-- 46656

-- (1.55 secs, 65,727,176 bytes)

Pensamiento

«El control de la complejidad es la esencia de la programación.» ~ B.W. Kernigan

Menor número triangular con más de n divisores

PorJosé A. Alonso 27-noviembre-20194-diciembre-2019

La sucesión de los números triangulares se obtiene sumando los números naturales.

   *     *      *        *         *   
        * *    * *      * *       * *  
              * * *    * * *     * * * 
                      * * * *   * * * *
                               * * * * * 
   1     3      6        10        15

* * * * *

* * * * * * * *

* * * * * * * * *

* * * * * * * *

* * * * *

1 3 6 10 15

Así, el 7º número triangular es

   1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

1	1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

Los primeros 10 números triangulares son

   1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

1	1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Los divisores de los primeros 7 números triangulares son:

    1: 1
    3: 1,3
    6: 1,2,3,6
   10: 1,2,5,10
   15: 1,3,5,15
   21: 1,3,7,21
   28: 1,2,4,7,14,28

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

Como se puede observar, 28 es el menor número triangular con más de 5 divisores.

Definir la función

   menorTriangularConAlMenosNDivisores :: Int -> Integer

1	menorTriangularConAlMenosNDivisores :: Int -> Integer

tal que (menorTriangularConAlMenosNDivisores n) es el menor número triangular que tiene al menos n divisores. Por ejemplo,

   menorTriangularConAlMenosNDivisores 5    ==  28
   menorTriangularConAlMenosNDivisores 50   ==  25200
   menorTriangularConAlMenosNDivisores 500  ==  76576500

menorTriangularConAlMenosNDivisores 5 == 28

menorTriangularConAlMenosNDivisores 50 == 25200

menorTriangularConAlMenosNDivisores 500 == 76576500

Nota: Este ejercicio está basado en el problema 12 del Proyecto Euler

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)

menorTriangularConAlMenosNDivisores :: Int -> Integer
menorTriangularConAlMenosNDivisores n = 
  head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones
-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y
-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

-- triangulares es la sucesión de los números triangulares. Por ejemplo,
--    take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
triangulares :: [Integer]
triangulares = scanl1 (+) [1..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,
--    nDivisores 28  ==  6
nDivisores :: Integer -> Int
nDivisores = product . map ((+1) . length) . group . primeFactors

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

menorTriangularConAlMenosNDivisores :: Int -> Integer

menorTriangularConAlMenosNDivisores n =

head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones

-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y

-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

-- triangulares es la sucesión de los números triangulares. Por ejemplo,

-- take 10 triangulares == [1,3,6,10,15,21,28,36,45,55]

triangulares :: [Integer]

triangulares = scanl1 (+) [1..]

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

-- nDivisores 28 == 6

nDivisores :: Integer -> Int

nDivisores = product . map ((+1) . length) . group . primeFactors

Pensamiento

«La Matemática es una ciencia experimental y la computación es el experimento.» ~ Rivin

Último dígito no nulo del factorial

PorJosé A. Alonso 12-noviembre-201920-noviembre-2019

El factorial de 7 es

   7! = 1 * 2 * 3 * 4 * 5 * 6 * 7 = 5040

1	7! = 1 * 2 * 3 * 4 * 5 * 6 * 7 = 5040

por tanto, el último dígito no nulo del factorial de 7 es 4.

Definir la función

   ultimoNoNuloFactorial :: Integer -> Integer

1	ultimoNoNuloFactorial :: Integer -> Integer

tal que (ultimoNoNuloFactorial n) es el último dígito no nulo del factorial de n. Por ejemplo,

   ultimoNoNuloFactorial  7  == 4
   ultimoNoNuloFactorial 10  == 8
   ultimoNoNuloFactorial 12  == 6
   ultimoNoNuloFactorial 97  == 2
   ultimoNoNuloFactorial  0  == 1

ultimoNoNuloFactorial 7 == 4

ultimoNoNuloFactorial 10 == 8

ultimoNoNuloFactorial 12 == 6

ultimoNoNuloFactorial 97 == 2

ultimoNoNuloFactorial 0 == 1

Comprobar con QuickCheck que si n es mayor que 4, entonces el último dígito no nulo del factorial de n es par.

Soluciones

import Data.Char (digitToInt)
import Test.QuickCheck

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

ultimoNoNuloFactorial :: Integer -> Integer
ultimoNoNuloFactorial = ultimoNoNulo . factorial

-- (factorial n) es el factorial de n. Por ejemplo,
--    factorial 7  ==  5040
factorial :: Integer -> Integer
factorial n = product [1..n]

-- (ultimoNoNulo n) es el último dígito no nulo de n. Por ejemplo,
--    ultimoNoNulo 5040  ==  4
ultimoNoNulo :: Integer -> Integer
ultimoNoNulo n | r /= 0    = r
               | otherwise = ultimoNoNulo q
  where (q,r) = n `quotRem` 10

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

ultimoNoNuloFactorial2 :: Integer -> Integer
ultimoNoNuloFactorial2 = last . filter (/= 0) . digitos . factorial

digitos :: Integer -> [Integer]
digitos n = [read [x] | x <- show n]

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

ultimoNoNuloFactorial3 :: Integer -> Integer
ultimoNoNuloFactorial3 = last . filter (/= 0) . digitos3 . factorial3

digitos3 :: Integer -> [Integer]
digitos3 = map (fromIntegral . digitToInt) . show

factorial3 :: Integer -> Integer
factorial3 = product . enumFromTo 1

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

ultimoNoNulo4 :: Integer -> Integer
ultimoNoNulo4 n = read [head (dropWhile (=='0') (reverse (show n)))]

-- 5ª solución
-- ===========

ultimoNoNulo5 :: Integer -> Integer
ultimoNoNulo5 =
  read . return . head . dropWhile ('0' ==) . reverse . show

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

-- La propiedad es
prop_ultimoNoNuloFactorial :: Integer -> Property
prop_ultimoNoNuloFactorial n = 
  n > 4 ==> even (ultimoNoNuloFactorial n)
                  
-- La comprobación es
--    ghci> quickCheck prop_ultimoNoNuloFactorial
--    +++ OK, passed 100 tests.

import Data.Char (digitToInt)

import Test.QuickCheck

-- 1ª solución

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

ultimoNoNuloFactorial :: Integer -> Integer

ultimoNoNuloFactorial = ultimoNoNulo . factorial

-- (factorial n) es el factorial de n. Por ejemplo,

-- factorial 7 == 5040

factorial :: Integer -> Integer

factorial n = product [1..n]

-- (ultimoNoNulo n) es el último dígito no nulo de n. Por ejemplo,

-- ultimoNoNulo 5040 == 4

ultimoNoNulo :: Integer -> Integer

ultimoNoNulo n | r /= 0 = r

| otherwise = ultimoNoNulo q

where (q,r) = n `quotRem` 10

-- 2ª solución

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