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

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

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.

Pensamiento

Busca el tu esencial,
que no está en ninguna parte
y en todas partes está.

Antonio Machado

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

PorJosé A. Alonso 14-mayo-201925-abril-2021

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

Soluciones

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

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

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

-- (cifras x) es la lista de las cifras del número x. Por ejemplo, 
--    cifras 325  ==  [3,2,5]
cifras :: Integer -> [Integer]
cifras 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
                 . cifras

-- 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 soluciones                                          --
-- ---------------------------------------------------------------------

-- Tiempo (en segundos) del cálculo de (mayorProducto 5 (product [1..n]))
-- 
--    | Def | 500  | 1000 | 5000 | 10000 
--    | 1   | 0.01 | 0.02 | 0.06 | 0.11
--    | 2   | 0.01 | 0.03 | 0.80 | 3.98
--    | 3   | 0.01 | 0.03 | 0.82 | 4.13
--    | 4   | 2.77 |      |      |

import Test.QuickCheck

import Data.List (inits, tails)

import Data.Char (digitToInt)

-- 1ª solución

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

mayorProducto1 :: Int -> Integer -> Integer

mayorProducto1 n x =

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

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

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

cifras :: Integer -> [Integer]

cifras 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

. cifras

-- 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 soluciones --

-- ---------------------------------------------------------------------

-- Tiempo (en segundos) del cálculo de (mayorProducto 5 (product [1..n]))

-- | Def | 500 | 1000 | 5000 | 10000

-- | 1 | 0.01 | 0.02 | 0.06 | 0.11

-- | 2 | 0.01 | 0.03 | 0.80 | 3.98

-- | 3 | 0.01 | 0.03 | 0.82 | 4.13

-- | 4 | 2.77 | | |

Pensamiento

A las palabras de amor
les sienta bien su poquito
de exageración.

Antonio Machado

Productos simultáneos de dos y tres números consecutivos

PorJosé A. Alonso 13-mayo-201925-abril-2021

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ª solución
-- ===========

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

-- 2ª solució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..])

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

-- La comparación es
--    λ> let (n,x) = (3,200) in productos1 n (product [x..x+n-1])
--    [[200,201,202]]
--    (7.38 secs, 7,235,956,440 bytes)
--    λ> let (n,x) = (3,200) in productos2 n (product [x..x+n-1])
--    [[200,201,202]]
--    (0.01 secs, 451,928 bytes)
--
--    λ> productos2 3 1000018000107000210
--    [[1000005,1000006,1000007]]
--    (1.57 secs, 1,560,159,144 bytes)

-- En lo que sigue se usa la 2ª definición
productos :: Integer -> Integer -> [[Integer]]
productos = productos2

-- La propiedad es
prop_productos :: Integer -> Integer -> Property
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.

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

-- El cálculo es
--    λ> take 2 productosDe2y3consecutivos
--    [6,210]

import Test.QuickCheck

-- 1ª solución

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

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

productos1 n x =

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

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

-- 2ª solució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..])

-- Comparación de eficiencia

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

-- La comparación es

-- λ> let (n,x) = (3,200) in productos1 n (product [x..x+n-1])

-- [[200,201,202]]

-- (7.38 secs, 7,235,956,440 bytes)

-- λ> let (n,x) = (3,200) in productos2 n (product [x..x+n-1])

-- [[200,201,202]]

-- (0.01 secs, 451,928 bytes)

-- λ> productos2 3 1000018000107000210

-- [[1000005,1000006,1000007]]

-- (1.57 secs, 1,560,159,144 bytes)

-- En lo que sigue se usa la 2ª definición

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

productos = productos2

-- La propiedad es

prop_productos :: Integer -> Integer -> Property

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.

productosDe2y3consecutivos :: [Integer]

productosDe2y3consecutivos =

[x | x <- [1..]

, let ys = productos 2 x

, not (null ys)

, let zs = productos 3 x

, not (null zs)]

-- El cálculo es

-- λ> take 2 productosDe2y3consecutivos

-- [6,210]

Pensamiento

Mis ojos en el espejo
son ojos ciegos que miran
los ojos con que los veo.

Antonio Machado

Matriz dodecafónica

PorJosé A. Alonso 10-abril-201925-abril-2019

Como se explica en Create a Twelve-Tone Melody With a Twelve-Tone Matrix una matriz dodecafónica es una matriz de 12 filas y 12 columnas construidas siguiendo los siguientes pasos:

Se escribe en la primera fila una permutación de los números del 1 al 12. Por ejemplo,

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( )

Escribir la primera columna de forma que, para todo i (entre 2 y 12), a(i,1) es el número entre 1 y 12 que verifica la siguiente condición

     (a(1,1) - a(i,1)) = (a(1,i) - a(1,1)) (módulo 12)

1	(a(1,1) - a(i,1)) = (a(1,i) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz con la 1ª fila y la 1ª columna es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5                                  )
     (  9                                  )
     (  1                                  )
     (  2                                  )
     ( 12                                  )
     ( 10                                  )
     ( 11                                  )
     (  6                                  )
     (  8                                  )
     (  7                                  )
     (  4                                  )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 )

( 9 )

( 1 )

( 2 )

( 12 )

( 10 )

( 11 )

( 6 )

( 8 )

( 7 )

( 4 )

Escribir la segunda fila de forma que, para todo j (entre 2 y 12), a(j,2) es el número entre 1 y 12 que verifica la siguiente condición

     (a(2,j) - a(1,j)) = (a(2,1) - a(1,1)) (módulo 12)

1	(a(2,j) - a(1,j)) = (a(2,1) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz con la 1ª fila, 1ª columna y 2ª fila es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5  3 11  7  6  8 10  9  2 12  1  4 )
     (  9                                  )
     (  1                                  )
     (  2                                  )
     ( 12                                  )
     ( 10                                  )
     ( 11                                  )
     (  6                                  )
     (  8                                  )
     (  7                                  )
     (  4                                  )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 )

( 1 )

( 2 )

( 12 )

( 10 )

( 11 )

( 6 )

( 8 )

( 7 )

( 4 )

Las restantes filas se completan como la 2ª; es decir, para todo i (entre 3 y 12) y todo j (entre 2 y 12), a(i,j) es el número entre 1 y 12 que verifica la siguiente relación.

     (a(i,j) - a(1,j)) = (a(i,1) - a(1,1)) (módulo 12)

1	(a(i,j) - a(1,j)) = (a(i,1) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz dodecafónica es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5  3 11  7  6  8 10  9  2 12  1  4 )
     (  9  7  3 11 10 12  2  1  6  4  5  8 )
     (  1 11  7  3  2  4  6  5 10  8  9 12 )
     (  2 12  8  4  3  5  7  6 11  9 10  1 )
     ( 12 10  6  2  1  3  5  4  9  7  8 11 )
     ( 10  8  4 12 11  1  3  2  7  5  6  9 )
     ( 11  9  5  1 12  2  4  3  8  6  7 10 )
     (  6  4 12  8  7  9 11 10  3  1  2  5 )
     (  8  6  2 10  9 11  1 12  5  3  4  7 )
     (  7  5  1  9  8 10 12 11  4  2  3  6 )
     (  4  2 10  6  5  7  9  8  1 11 12  3 )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 7 3 11 10 12 2 1 6 4 5 8 )

( 1 11 7 3 2 4 6 5 10 8 9 12 )

( 2 12 8 4 3 5 7 6 11 9 10 1 )

( 12 10 6 2 1 3 5 4 9 7 8 11 )

( 10 8 4 12 11 1 3 2 7 5 6 9 )

( 11 9 5 1 12 2 4 3 8 6 7 10 )

( 6 4 12 8 7 9 11 10 3 1 2 5 )

( 8 6 2 10 9 11 1 12 5 3 4 7 )

( 7 5 1 9 8 10 12 11 4 2 3 6 )

( 4 2 10 6 5 7 9 8 1 11 12 3 )

Definir la función

   matrizDodecafonica :: [Int] -> Matrix Int

1	matrizDodecafonica :: [Int] -> Matrix Int

tal que (matrizDodecafonica xs) es la matriz dodecafónica cuya primera fila es xs (que se supone que es una permutación de los números del 1 al 12). Por ejemplo,

   λ> matrizDodecafonica [3,1,9,5,4,6,8,7,12,10,11,2]
   (  3  1  9  5  4  6  8  7 12 10 11  2 )
   (  5  3 11  7  6  8 10  9  2 12  1  4 )
   (  9  7  3 11 10 12  2  1  6  4  5  8 )
   (  1 11  7  3  2  4  6  5 10  8  9 12 )
   (  2 12  8  4  3  5  7  6 11  9 10  1 )
   ( 12 10  6  2  1  3  5  4  9  7  8 11 )
   ( 10  8  4 12 11  1  3  2  7  5  6  9 )
   ( 11  9  5  1 12  2  4  3  8  6  7 10 )
   (  6  4 12  8  7  9 11 10  3  1  2  5 )
   (  8  6  2 10  9 11  1 12  5  3  4  7 )
   (  7  5  1  9  8 10 12 11  4  2  3  6 )
   (  4  2 10  6  5  7  9  8  1 11 12  3 )

λ> matrizDodecafonica [3,1,9,5,4,6,8,7,12,10,11,2]

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 7 3 11 10 12 2 1 6 4 5 8 )

( 1 11 7 3 2 4 6 5 10 8 9 12 )

( 2 12 8 4 3 5 7 6 11 9 10 1 )

( 12 10 6 2 1 3 5 4 9 7 8 11 )

( 10 8 4 12 11 1 3 2 7 5 6 9 )

( 11 9 5 1 12 2 4 3 8 6 7 10 )

( 6 4 12 8 7 9 11 10 3 1 2 5 )

( 8 6 2 10 9 11 1 12 5 3 4 7 )

( 7 5 1 9 8 10 12 11 4 2 3 6 )

( 4 2 10 6 5 7 9 8 1 11 12 3 )

Comprobar con QuickCheck para toda matriz dodecafónica D se verifican las siguientes propiedades:

todas las filas de D son permutaciones de los números 1 a 12,
todos los elementos de la diagonal de D son iguales y
la suma de todos los elementos de D es 936.

Nota: Este ejercicio ha sido propuesto por Francisco J. Hidalgo.

Soluciones

import Data.List
import Test.QuickCheck
import Data.Matrix

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

matrizDodecafonica :: [Int] -> Matrix Int
matrizDodecafonica xs = matrix 12 12 f
  where f (1,j) = xs !! (j-1)
        f (i,1) = modulo12 (2 * f (1,1) - f (1,i)) 
        f (i,j) = modulo12 (f (1,j) + f (i,1) - f (1,1)) 
        modulo12 0  = 12
        modulo12 12 = 12
        modulo12 x  = x `mod` 12
  
-- 2ª solución
-- ===========

matrizDodecafonica2 :: [Int] -> Matrix Int
matrizDodecafonica2 xs = fromLists (secuencias xs)

secuencias :: [Int] -> [[Int]]
secuencias xs = [secuencia a xs | a <- inversa xs]

inversa :: [Int] -> [Int]
inversa xs = map conv (map (\x -> (-x) + 2* (abs a)) xs)
  where a = head xs
        
secuencia :: Int -> [Int] -> [Int]
secuencia n xs = [conv (a+(n-b)) | a <- xs] 
  where b = head xs

conv :: Int -> Int
conv n | n == 0 = 12
       | n < 0 = conv (n+12)
       | n > 11 = conv (mod n 12)
       | otherwise = n          

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

-- Las propiedades son
prop_dodecafonica :: Int -> Property
prop_dodecafonica n = 
  n >= 0 ==>
  all esPermutacion (toLists d)
  && all (== d!(1,1)) [d!(i,i) | i <- [2..12]]
  && sum d == 936
  where xss = permutations [1..12]
        k   = n `mod` product [1..12]
        d   = matrizDodecafonica (xss !! k)
        esPermutacion ys = sort ys == [1..12]

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

import Data.List

import Test.QuickCheck

import Data.Matrix

-- 1ª solución

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

matrizDodecafonica :: [Int] -> Matrix Int

matrizDodecafonica xs = matrix 12 12 f

where f (1,j) = xs !! (j-1)

f (i,1) = modulo12 (2 * f (1,1) - f (1,i))

f (i,j) = modulo12 (f (1,j) + f (i,1) - f (1,1))

modulo12 0 = 12

modulo12 12 = 12

modulo12 x = x `mod` 12

-- 2ª solución

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

matrizDodecafonica2 :: [Int] -> Matrix Int

matrizDodecafonica2 xs = fromLists (secuencias xs)

secuencias :: [Int] -> [[Int]]

secuencias xs = [secuencia a xs | a <- inversa xs]

inversa :: [Int] -> [Int]

inversa xs = map conv (map (\x -> (-x) + 2* (abs a)) xs)

where a = head xs

secuencia :: Int -> [Int] -> [Int]

secuencia n xs = [conv (a+(n-b)) | a <- xs]

where b = head xs

conv :: Int -> Int

conv n | n == 0 = 12

| n < 0 = conv (n+12)

| n > 11 = conv (mod n 12)

| otherwise = n

-- Propiedades

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

-- Las propiedades son

prop_dodecafonica :: Int -> Property

prop_dodecafonica n =

n >= 0 ==>

all esPermutacion (toLists d)

&& all (== d!(1,1)) [d!(i,i) | i <- [2..12]]

&& sum d == 936

where xss = permutations [1..12]

k = n `mod` product [1..12]

d = matrizDodecafonica (xss !! k)

esPermutacion ys = sort ys == [1..12]

-- La comprobación es

-- λ> quickCheck prop_dodecafonica

-- +++ OK, passed 100 tests.

Pensamiento

Como el olivar,
mucho fruto lleva,
poca sombra da.

Antonio Machado

Números altamente compuestos

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

Un número altamente compuesto es un entero positivo con más divisores que cualquier entero positivo más pequeño. Por ejemplo,

4 es un número altamente compuesto porque es el menor con 3 divisores,
5 no es altamente compuesto porque tiene menos divisores que 4 y
6 es un número altamente compuesto porque es el menor con 4 divisores,

Los primeros números altamente compuestos son

   1, 2, 4, 6, 12, 24, 36, 48, 60, 120, 180, 240, 360, ...

1	1, 2, 4, 6, 12, 24, 36, 48, 60, 120, 180, 240, 360, ...

Definir las funciones

   esAltamenteCompuesto       :: Int -> Bool
   altamenteCompuestos        :: [Int]
   graficaAltamenteCompuestos :: Int -> IO ()

esAltamenteCompuesto :: Int -> Bool

altamenteCompuestos :: [Int]

graficaAltamenteCompuestos :: Int -> IO ()

tales que

(esAltamanteCompuesto x) se verifica si x es altamente compuesto. Por ejemplo,

     esAltamenteCompuesto 4      ==  True
     esAltamenteCompuesto 5      ==  False
     esAltamenteCompuesto 6      ==  True
     esAltamenteCompuesto 1260   ==  True
     esAltamenteCompuesto 2520   ==  True
     esAltamenteCompuesto 27720  ==  True

esAltamenteCompuesto 4 == True

esAltamenteCompuesto 5 == False

esAltamenteCompuesto 6 == True

esAltamenteCompuesto 1260 == True

esAltamenteCompuesto 2520 == True

esAltamenteCompuesto 27720 == True

altamente compuestos es la sucesión de los números altamente compuestos. Por ejemplo,

     λ> take 20 altamenteCompuestos
     [1,2,4,6,12,24,36,48,60,120,180,240,360,720,840,1260,1680,2520,5040,7560]

1 2	λ> take 20 altamenteCompuestos [1,2,4,6,12,24,36,48,60,120,180,240,360,720,840,1260,1680,2520,5040,7560]

(graficaAltamenteCompuestos n) dibuja la gráfica de los n primeros números altamente compuestos. Por ejemplo, (graficaAltamenteCompuestos 25) dibuja

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple

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

esAltamenteCompuesto :: Int -> Bool
esAltamenteCompuesto x =
  and [nDivisores x > nDivisores y | y <- [1..x-1]]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,
--    nDivisores 30  ==  8
nDivisores :: Int -> Int
nDivisores x = length (divisores x)

-- (divisores x) es la lista de los divisores de x. Por ejemplo,
--    divisores 30  ==  [1,2,3,5,6,10,15,30]
divisores :: Int -> [Int]
divisores x =
  [y | y <- [1..x]
     , x `mod` y == 0]

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

esAltamenteCompuesto2 :: Int -> Bool
esAltamenteCompuesto2 x =
  all (nDivisores2 x >) [nDivisores2 y | y <- [1..x-1]]

-- (nDivisores2 x) es el número de divisores de x. Por ejemplo,
--    nDivisores2 30  ==  8
nDivisores2 :: Int -> Int
nDivisores2 = succ . length . divisoresPropios

-- (divisoresPropios x) es la lista de los divisores de x menores que
-- x. Por ejemplo, 
--    divisoresPropios 30  ==  [1,2,3,5,6,10,15]
divisoresPropios :: Int -> [Int]
divisoresPropios x =
  [y | y <- [1..x `div` 2]
     , x `mod` y == 0]

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

esAltamenteCompuesto3 :: Int -> Bool
esAltamenteCompuesto3 x =
  all (nDivisores3 x >) [nDivisores3 y | y <- [1..x-1]]

-- (nDivisores3 x) es el número de divisores de x. Por ejemplo,
--    nDivisores3 30  ==  8
nDivisores3 :: Int -> Int
nDivisores3 x =
  product [1 + length xs | xs <- group (primeFactors x)]

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

esAltamenteCompuesto4 :: Int -> Bool
esAltamenteCompuesto4 x =
  x `pertenece` altamenteCompuestos2

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

altamenteCompuestos :: [Int]
altamenteCompuestos =
  filter esAltamenteCompuesto4 [1..]

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

altamenteCompuestos2 :: [Int]
altamenteCompuestos2 =
  1 : [y | ((x,n),(y,m)) <- zip sucMaxDivisores (tail sucMaxDivisores)
         , m > n]

-- sucMaxDivisores es la sucesión formada por los números enteros
-- positivos y el máximo número de divisores hasta cada número. Por
-- ejemplo,
--    λ> take 12 sucMaxDivisores
--    [(1,1),(2,2),(3,2),(4,3),(5,3),(6,4),(7,4),(8,4),(9,4),(10,4),(11,4),(12,6)]
sucMaxDivisores :: [(Int,Int)]
sucMaxDivisores =
  zip [1..] (scanl1 max (map nDivisores3 [1..]))

pertenece :: Int -> [Int] -> Bool
pertenece x ys =
  x == head (dropWhile (<x) ys)

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

--    λ> esAltamenteCompuesto 1260
--    True
--    (2.99 secs, 499,820,296 bytes)
--    λ> esAltamenteCompuesto2 1260
--    True
--    (0.51 secs, 83,902,744 bytes)
--    λ> esAltamenteCompuesto3 1260
--    True
--    (0.04 secs, 15,294,192 bytes)
--    λ> esAltamenteCompuesto4 1260
--    True
--    (0.04 secs, 15,594,392 bytes)
--    
--    λ> esAltamenteCompuesto2 2520
--    True
--    (2.10 secs, 332,940,168 bytes)
--    λ> esAltamenteCompuesto3 2520
--    True
--    (0.09 secs, 37,896,168 bytes)
--    λ> esAltamenteCompuesto4 2520
--    True
--    (0.06 secs, 23,087,456 bytes)
--
--    λ> esAltamenteCompuesto3 27720
--    True
--    (1.32 secs, 841,010,624 bytes)
--    λ> esAltamenteCompuesto4 27720
--    True
--    (1.33 secs, 810,870,384 bytes)

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

--    λ> altamenteCompuestos !! 25
--    45360
--    (2.84 secs, 1,612,045,976 bytes)
--    λ> altamenteCompuestos2 !! 25
--    45360
--    (0.01 secs, 102,176 bytes)

-- Definición de graficaAltamenteCompuestos
-- ========================================

graficaAltamenteCompuestos :: Int -> IO ()
graficaAltamenteCompuestos n =
  plotList [ Key Nothing
           , PNG ("Numeros_altamente_compuestos.png")
           ]
           (take n altamenteCompuestos2)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple

-- 1ª definición de esAltamenteCompuesto

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

esAltamenteCompuesto :: Int -> Bool

esAltamenteCompuesto x =

and [nDivisores x > nDivisores y | y <- [1..x-1]]

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

-- nDivisores 30 == 8

nDivisores :: Int -> Int

nDivisores x = length (divisores x)

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

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

divisores :: Int -> [Int]

divisores x =

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

, x `mod` y == 0]

-- 2ª definición de esAltamenteCompuesto

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

esAltamenteCompuesto2 :: Int -> Bool

esAltamenteCompuesto2 x =

all (nDivisores2 x >) [nDivisores2 y | y <- [1..x-1]]

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

-- nDivisores2 30 == 8

nDivisores2 :: Int -> Int

nDivisores2 = succ . length . divisoresPropios

-- (divisoresPropios x) es la lista de los divisores de x menores que

-- x. Por ejemplo,

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

divisoresPropios :: Int -> [Int]

divisoresPropios x =

[y | y <- [1..x `div` 2]

, x `mod` y == 0]

-- 3ª definición de esAltamenteCompuesto

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

esAltamenteCompuesto3 :: Int -> Bool

esAltamenteCompuesto3 x =

all (nDivisores3 x >) [nDivisores3 y | y <- [1..x-1]]

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

-- nDivisores3 30 == 8

nDivisores3 :: Int -> Int

nDivisores3 x =

product [1 + length xs | xs <- group (primeFactors x)]

-- 4ª definición de esAltamenteCompuesto

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

esAltamenteCompuesto4 :: Int -> Bool

esAltamenteCompuesto4 x =

x `pertenece` altamenteCompuestos2

-- 1ª definición de altamenteCompuestos

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

altamenteCompuestos :: [Int]

altamenteCompuestos =

filter esAltamenteCompuesto4 [1..]

-- 2ª definición de altamenteCompuestos

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

altamenteCompuestos2 :: [Int]

altamenteCompuestos2 =

1 : [y | ((x,n),(y,m)) <- zip sucMaxDivisores (tail sucMaxDivisores)

, m > n]

-- sucMaxDivisores es la sucesión formada por los números enteros

-- positivos y el máximo número de divisores hasta cada número. Por

-- ejemplo,

-- λ> take 12 sucMaxDivisores

-- [(1,1),(2,2),(3,2),(4,3),(5,3),(6,4),(7,4),(8,4),(9,4),(10,4),(11,4),(12,6)]

sucMaxDivisores :: [(Int,Int)]

sucMaxDivisores =

zip [1..] (scanl1 max (map nDivisores3 [1..]))

pertenece :: Int -> [Int] -> Bool

pertenece x ys =

x == head (dropWhile (<x) ys)

-- Comparación de eficiencia de esAltamenteCompuesto

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

-- λ> esAltamenteCompuesto 1260

-- True

-- (2.99 secs, 499,820,296 bytes)

-- λ> esAltamenteCompuesto2 1260

-- True

-- (0.51 secs, 83,902,744 bytes)

-- λ> esAltamenteCompuesto3 1260

-- True

-- (0.04 secs, 15,294,192 bytes)

-- λ> esAltamenteCompuesto4 1260

-- True

-- (0.04 secs, 15,594,392 bytes)

-- λ> esAltamenteCompuesto2 2520

-- True

-- (2.10 secs, 332,940,168 bytes)

-- λ> esAltamenteCompuesto3 2520

-- True

-- (0.09 secs, 37,896,168 bytes)

-- λ> esAltamenteCompuesto4 2520

-- True

-- (0.06 secs, 23,087,456 bytes)

-- λ> esAltamenteCompuesto3 27720

-- True

-- (1.32 secs, 841,010,624 bytes)

-- λ> esAltamenteCompuesto4 27720

-- True

-- (1.33 secs, 810,870,384 bytes)

-- Comparación de eficiencia de altamenteCompuestos

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

-- λ> altamenteCompuestos !! 25

-- 45360

-- (2.84 secs, 1,612,045,976 bytes)

-- λ> altamenteCompuestos2 !! 25

-- 45360

-- (0.01 secs, 102,176 bytes)

-- Definición de graficaAltamenteCompuestos

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

graficaAltamenteCompuestos :: Int -> IO ()

graficaAltamenteCompuestos n =

plotList [ Key Nothing

, PNG ("Numeros_altamente_compuestos.png")

]

(take n altamenteCompuestos2)

Pensamiento

Nuestras horas son minutos
cuando esperamos saber,
y siglos cuando sabemos
lo que se puede aprender.

Antonio Machado

Árbol de subconjuntos

PorJosé A. Alonso 28-diciembre-20184-enero-2019

Se dice que A es un subconjunto maximal de B si A ⊂ B y no existe ningún C tal que A ⊂ C y C ⊂ B. Por ejemplo, {2,5} es un subconjunto maximal de {2,3,5], pero {3] no lo es.

El árbol de los subconjuntos de un conjunto A es el árbol que tiene como raíz el conjunto A y cada nodo tiene como hijos sus subconjuntos maximales. Por ejemplo, el árbol de subconjuntos de [2,3,5] es

          [2, 3, 5]
          /   |  \
         /    |   \  
        /     |    \   
       /      |     \
      /       |      \
    [5,3]   [2,3]   [2,5]  
    /  \    /  \    /  \  
   [3] [5] [3] [2] [5] [2]
    |   |   |   |   |   | 
   [ ] [ ] [ ] [ ] [ ] [ ]

[2, 3, 5]

/ | \

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

/ \ / \ / \

[3] [5] [3] [2] [5] [2]

| | | | | |

[ ] [ ] [ ] [ ] [ ] [ ]

Usando el tipo de dato

   data Arbol = N [Integer] [Arbol]
     deriving (Eq, Show)

1 2	data Arbol = N [Integer] [Arbol] deriving (Eq, Show)

el árbol anterior se representa por

   N [2,5,3]
     [N [5,3]
        [N [3]
           [N [] []],
         N [5]
           [N [] []]],
      N [2,3]
        [N [3]
           [N [] []],
         N [2]
           [N [] []]],
      N [2,5]
        [N [5]
           [N [] []],
         N [2]
           [N [] []]]]

N [2,5,3]

[N [5,3]

[N [3]

[N [] []],

N [5]

[N [] []]],

N [2,3]

[N [3]

[N [] []],

N [2]

[N [] []]],

N [2,5]

[N [5]

[N [] []],

N [2]

[N [] []]]]

Definir las funciones

   arbolSubconjuntos :: [Int] -> Arbol 
   nOcurrenciasArbolSubconjuntos :: [Int] -> [Int] -> Int

1 2	arbolSubconjuntos :: [Int] -> Arbol nOcurrenciasArbolSubconjuntos :: [Int] -> [Int] -> Int

tales que

(arbolSubconjuntos x) es el árbol de los subconjuntos de xs. Por ejemplo,

     λ> arbolSubconjuntos [2,5,3]
     N [2,5,3] [N [5,3] [N [3] [N [] []],N [5] [N [] []]],
                N [2,3] [N [3] [N [] []],N [2] [N [] []]],
                N [2,5] [N [5] [N [] []],N [2] [N [] []]]]

λ> arbolSubconjuntos [2,5,3]

N [2,5,3] [N [5,3] [N [3] [N [] []],N [5] [N [] []]],

N [2,3] [N [3] [N [] []],N [2] [N [] []]],

N [2,5] [N [5] [N [] []],N [2] [N [] []]]]

(nOcurrenciasArbolSubconjuntos xs ys) es el número de veces que aparece el conjunto xs en el árbol de los subconjuntos de ys. Por ejemplo,

     nOcurrenciasArbolSubconjuntos []      [2,5,3]  ==  6
     nOcurrenciasArbolSubconjuntos [3]     [2,5,3]  ==  2
     nOcurrenciasArbolSubconjuntos [3,5]   [2,5,3]  ==  1
     nOcurrenciasArbolSubconjuntos [3,5,2] [2,5,3]  ==  1

nOcurrenciasArbolSubconjuntos [] [2,5,3] == 6

nOcurrenciasArbolSubconjuntos [3] [2,5,3] == 2

nOcurrenciasArbolSubconjuntos [3,5] [2,5,3] == 1

nOcurrenciasArbolSubconjuntos [3,5,2] [2,5,3] == 1

Comprobar con QuickChek que, para todo entero positivo n, el número de ocurrencia de un subconjunto xs de [1..n] en el árbol de los subconjuntos de [1..n] es el factorial de n-k (donde k es el número de elementos de xs).

Soluciones

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

data Arbol = N [Int] [Arbol]
  deriving (Eq, Show)

arbolSubconjuntos :: [Int] -> Arbol 
arbolSubconjuntos xs =
  N xs (map arbolSubconjuntos (subconjuntosMaximales xs))

-- (subconjuntosMaximales xs) es la lista de los subconjuntos maximales
-- de xs. Por ejemplo,
--    subconjuntosMaximales [2,5,3]  ==  [[5,3],[2,3],[2,5]]
subconjuntosMaximales :: [Int] -> [[Int]]
subconjuntosMaximales xs =
  [delete x xs | x <- xs]

-- Definición de nOcurrenciasArbolSubconjuntos
-- ===========================================

nOcurrenciasArbolSubconjuntos :: [Int] -> [Int] -> Int
nOcurrenciasArbolSubconjuntos xs ys =
  nOcurrencias xs (arbolSubconjuntos ys)

-- (nOcurrencias x a) es el número de veces que aparece x en el árbol
-- a. Por ejemplo,
--    nOcurrencias 3 (arbolSubconjuntos 30)  ==  2
nOcurrencias :: [Int] -> Arbol -> Int
nOcurrencias xs (N ys [])
  | conjunto xs == conjunto ys  = 1
  | otherwise                   = 0
nOcurrencias xs (N ys zs)
  | conjunto xs == conjunto ys = 1 + sum [nOcurrencias xs z | z <- zs]
  | otherwise                  = sum [nOcurrencias xs z | z <- zs]

-- (conjunto xs) es el conjunto ordenado correspondiente a xs. Por
-- ejemplo, 
--    conjunto [3,2,5,2,3,7,2]  ==  [2,3,5,7]
conjunto :: [Int] -> [Int]
conjunto = nub . sort

-- La propiedad es
prop_nOcurrencias :: (Positive Int) -> [Int] -> Bool
prop_nOcurrencias (Positive n) xs =
  nOcurrenciasArbolSubconjuntos ys [1..n] == factorial (n-k)
  where ys = nub [1 + x `mod` n | x <- xs]
        k  = length ys
        factorial m = product [1..m]

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

import Data.List (delete, nub, sort, subsequences)

import Test.QuickCheck

data Arbol = N [Int] [Arbol]

deriving (Eq, Show)

arbolSubconjuntos :: [Int] -> Arbol

arbolSubconjuntos xs =

N xs (map arbolSubconjuntos (subconjuntosMaximales xs))

-- (subconjuntosMaximales xs) es la lista de los subconjuntos maximales

-- de xs. Por ejemplo,

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

subconjuntosMaximales :: [Int] -> [[Int]]

subconjuntosMaximales xs =

[delete x xs | x <- xs]

-- Definición de nOcurrenciasArbolSubconjuntos

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

nOcurrenciasArbolSubconjuntos :: [Int] -> [Int] -> Int

nOcurrenciasArbolSubconjuntos xs ys =

nOcurrencias xs (arbolSubconjuntos ys)

-- (nOcurrencias x a) es el número de veces que aparece x en el árbol

-- a. Por ejemplo,

-- nOcurrencias 3 (arbolSubconjuntos 30) == 2

nOcurrencias :: [Int] -> Arbol -> Int

nOcurrencias xs (N ys [])

| conjunto xs == conjunto ys = 1

| otherwise = 0

nOcurrencias xs (N ys zs)

| conjunto xs == conjunto ys = 1 + sum [nOcurrencias xs z | z <- zs]

| otherwise = sum [nOcurrencias xs z | z <- zs]

-- (conjunto xs) es el conjunto ordenado correspondiente a xs. Por

-- ejemplo,

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

conjunto :: [Int] -> [Int]

conjunto = nub . sort

-- La propiedad es

prop_nOcurrencias :: (Positive Int) -> [Int] -> Bool

prop_nOcurrencias (Positive n) xs =

nOcurrenciasArbolSubconjuntos ys [1..n] == factorial (n-k)

where ys = nub [1 + x `mod` n | x <- xs]

k = length ys

factorial m = product [1..m]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=9}) prop_nOcurrencias

-- +++ OK, passed 100 tests.

Pensamiento

Nunca traces tu frontera,
ni cuides de tu perfil;
todo eso es cosa de fuera.

Antonio Machado

Número de divisores compuestos

PorJosé A. Alonso 25-diciembre-20181-enero-2019

Definir la función

   nDivisoresCompuestos :: Integer -> Integer

1	nDivisoresCompuestos :: Integer -> Integer

tal que (nDivisoresCompuestos x) es el número de divisores de x que son compuestos (es decir, números mayores que 1 que no son primos). Por ejemplo,

   nDivisoresCompuestos 30  ==  4
   nDivisoresCompuestos (product [1..11])  ==  534
   nDivisoresCompuestos (product [1..25])  ==  340022
   length (show (nDivisoresCompuestos (product [1..3*10^4]))) == 1948

nDivisoresCompuestos 30 == 4

nDivisoresCompuestos (product [1..11]) == 534

nDivisoresCompuestos (product [1..25]) == 340022

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

Soluciones

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

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

--    nDivisoresCompuestos 30  ==  4
nDivisoresCompuestos :: Integer -> Integer
nDivisoresCompuestos =
  genericLength . divisoresCompuestos

-- (divisoresCompuestos x) es la lista de los divisores de x que
-- son números compuestos (es decir, números mayores que 1 que no son
-- primos). Por ejemplo,
--    divisoresCompuestos 30  ==  [6,10,15,30]
divisoresCompuestos :: Integer -> [Integer]
divisoresCompuestos =
  sort
  . map product
  . compuestos
  . map concat
  . productoCartesiano
  . map inits
  . group
  . primeFactors
  where compuestos xss = [xs | xs <- xss, length xs > 1]  

-- (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]
  
-- 2ª solución
-- ===========

nDivisoresCompuestos2 :: Integer -> Integer
nDivisoresCompuestos2 x =
  nDivisores x - nDivisoresPrimos x - 1

-- (nDivisores x) es el número de divisores de x. Por ejemplo, 
--    nDivisores 30  ==  8
nDivisores :: Integer -> Integer
nDivisores x =
  product [1 + genericLength xs | xs <- group (primeFactors x)]

-- (nDivisoresPrimos x) es el número de divisores primos de x. Por
-- ejemplo,  
--    nDivisoresPrimos 30  ==  3
nDivisoresPrimos :: Integer -> Integer
nDivisoresPrimos =
  genericLength . group . primeFactors 

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

nDivisoresCompuestos3 :: Integer -> Integer
nDivisoresCompuestos3 x =
  nDivisores - nDivisoresPrimos - 1
  where xss              = group (primeFactors x)
        nDivisores       = product [1 + genericLength xs | xs <-xss]
        nDivisoresPrimos = genericLength xss

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

-- La propiedad es
prop_nDivisoresCompuestos :: (Positive Integer) -> Bool
prop_nDivisoresCompuestos (Positive x) =
  all (== nDivisoresCompuestos x) [f x | f <- [ nDivisoresCompuestos2
                                              , nDivisoresCompuestos3 ]]

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

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

--    λ> nDivisoresCompuestos (product [1..25])
--    340022
--    (2.04 secs, 2,232,146,776 bytes)
--    λ> nDivisoresCompuestos2 (product [1..25])
--    340022
--    (0.00 secs, 220,192 bytes)
--    
--    λ> length (show (nDivisoresCompuestos2 (product [1..3*10^4])))
--    1948
--    (5.22 secs, 8,431,630,288 bytes)
--    λ> length (show (nDivisoresCompuestos3 (product [1..3*10^4])))
--    1948
--    (3.06 secs, 4,662,277,664 bytes)

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª solución

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

-- nDivisoresCompuestos 30 == 4

nDivisoresCompuestos :: Integer -> Integer

nDivisoresCompuestos =

genericLength . divisoresCompuestos

-- (divisoresCompuestos x) es la lista de los divisores de x que

-- son números compuestos (es decir, números mayores que 1 que no son

-- primos). Por ejemplo,

-- divisoresCompuestos 30 == [6,10,15,30]

divisoresCompuestos :: Integer -> [Integer]

divisoresCompuestos =

sort

. map product

. compuestos

. map concat

. productoCartesiano

. map inits

. group

. primeFactors

where compuestos xss = [xs | xs <- xss, length xs > 1]

-- (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]

-- 2ª solución

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

nDivisoresCompuestos2 :: Integer -> Integer

nDivisoresCompuestos2 x =

nDivisores x - nDivisoresPrimos x - 1

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

-- nDivisores 30 == 8

nDivisores :: Integer -> Integer

nDivisores x =

product [1 + genericLength xs | xs <- group (primeFactors x)]

-- (nDivisoresPrimos x) es el número de divisores primos de x. Por

-- ejemplo,

-- nDivisoresPrimos 30 == 3

nDivisoresPrimos :: Integer -> Integer

nDivisoresPrimos =

genericLength . group . primeFactors

-- 3ª solución

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

nDivisoresCompuestos3 :: Integer -> Integer

nDivisoresCompuestos3 x =

nDivisores - nDivisoresPrimos - 1

where xss = group (primeFactors x)

nDivisores = product [1 + genericLength xs | xs <-xss]

nDivisoresPrimos = genericLength xss

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_nDivisoresCompuestos :: (Positive Integer) -> Bool

prop_nDivisoresCompuestos (Positive x) =

all (== nDivisoresCompuestos x) [f x | f <- [ nDivisoresCompuestos2

, nDivisoresCompuestos3 ]]

-- La comprobación es

-- λ> quickCheck prop_nDivisoresCompuestos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> nDivisoresCompuestos (product [1..25])

-- 340022

-- (2.04 secs, 2,232,146,776 bytes)

-- λ> nDivisoresCompuestos2 (product [1..25])

-- 340022

-- (0.00 secs, 220,192 bytes)

-- λ> length (show (nDivisoresCompuestos2 (product [1..3*10^4])))

-- 1948

-- (5.22 secs, 8,431,630,288 bytes)

-- λ> length (show (nDivisoresCompuestos3 (product [1..3*10^4])))

-- 1948

-- (3.06 secs, 4,662,277,664 bytes)

Pensamiento

«Lo corriente en el hombre es la tendencia a creer verdadero cuanto le
reporta alguna utilidad. Por eso hay tantos hombres capaces de comulgar
con ruedas de molino.»

Antonio Machado

Divisores compuestos

PorJosé A. Alonso 24-diciembre-201831-diciembre-2018

Definir la función

   divisoresCompuestos :: Integer -> [Integer]

1	divisoresCompuestos :: Integer -> [Integer]

tal que (divisoresCompuestos x) es la lista de los divisores de x que son números compuestos (es decir, números mayores que 1 que no son primos). Por ejemplo,

   divisoresCompuestos 30  ==  [6,10,15,30]
   length (divisoresCompuestos (product [1..11]))  ==  534
   length (divisoresCompuestos (product [1..14]))  ==  2585
   length (divisoresCompuestos (product [1..16]))  ==  5369
   length (divisoresCompuestos (product [1..25]))  ==  340022

divisoresCompuestos 30 == [6,10,15,30]

length (divisoresCompuestos (product [1..11])) == 534

length (divisoresCompuestos (product [1..14])) == 2585

length (divisoresCompuestos (product [1..16])) == 5369

length (divisoresCompuestos (product [1..25])) == 340022

Soluciones

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

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

divisoresCompuestos :: Integer -> [Integer]
divisoresCompuestos x =
  [y | y <- divisores x
     , y > 1
     , not (isPrime y)]

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

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

divisoresCompuestos2 :: Integer -> [Integer]
divisoresCompuestos2 x =
  [y | y <- divisores2 x
     , y > 1
     , not (isPrime y)]

-- (divisores2 x) es la lista de los divisores de x. Por ejemplo,
--    divisores2 30  ==  [1,2,3,5,6,10,15,30]
divisores2 :: Integer -> [Integer]
divisores2 x =
  [y | y <- [1..x `div` 2], x `mod` y == 0] ++ [x] 

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

divisoresCompuestos3 :: Integer -> [Integer]
divisoresCompuestos3 x =
  [y | y <- divisores2 x
     , y > 1
     , not (isPrime y)]

-- (divisores3 x) es la lista de los divisores de x. Por ejemplo,
--    divisores2 30  ==  [1,2,3,5,6,10,15,30]
divisores3 :: Integer -> [Integer]
divisores3 x =
  nub (sort (ys ++ [x `div` y | y <- ys]))
  where ys = [y | y <- [1..floor (sqrt (fromIntegral x))]
                , x `mod` y == 0]
             
-- 4ª solución
-- ===========

divisoresCompuestos4 :: Integer -> [Integer]
divisoresCompuestos4 x =
  [y | y <- divisores4 x
     , y > 1
     , not (isPrime y)]

-- (divisores4 x) es la lista de los divisores de x. Por ejemplo,
--    divisores4 30  ==  [1,2,3,5,6,10,15,30]
divisores4 :: Integer -> [Integer]
divisores4 =
  nub . sort . map product . subsequences . primeFactors

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

divisoresCompuestos5 :: Integer -> [Integer]
divisoresCompuestos5 x =
  [y | y <- divisores5 x
     , y > 1
     , not (isPrime y)]

-- (divisores5 x) es la lista de los divisores de x. Por ejemplo,
--    divisores5 30  ==  [1,2,3,5,6,10,15,30]
divisores5 :: Integer -> [Integer]
divisores5 =
  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]
  
-- 6ª solución
-- ===========

divisoresCompuestos6 :: Integer -> [Integer]
divisoresCompuestos6 =
  sort
  . map product
  . compuestos
  . map concat
  . productoCartesiano
  . map inits
  . group
  . primeFactors
  where compuestos xss = [xs | xs <- xss, length xs > 1]  

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

-- La propiedad es
prop_divisoresCompuestos :: (Positive Integer) -> Bool
prop_divisoresCompuestos (Positive x) =
  all (== divisoresCompuestos x) [f x | f <- [ divisoresCompuestos2
                                             , divisoresCompuestos3
                                             , divisoresCompuestos4
                                             , divisoresCompuestos5
                                             , divisoresCompuestos6 ]]

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

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

--    λ> length (divisoresCompuestos (product [1..11]))
--    534
--    (14.59 secs, 7,985,108,976 bytes)
--    λ> length (divisoresCompuestos2 (product [1..11]))
--    534
--    (7.36 secs, 3,993,461,168 bytes)
--    λ> length (divisoresCompuestos3 (product [1..11]))
--    534
--    (7.35 secs, 3,993,461,336 bytes)
--    λ> length (divisoresCompuestos4 (product [1..11]))
--    534
--    (0.07 secs, 110,126,392 bytes)
--    λ> length (divisoresCompuestos5 (product [1..11]))
--    534
--    (0.01 secs, 3,332,224 bytes)
--    λ> length (divisoresCompuestos6 (product [1..11]))
--    534
--    (0.01 secs, 1,869,776 bytes)
--    
--    λ> length (divisoresCompuestos4 (product [1..14]))
--    2585
--    (9.11 secs, 9,461,570,720 bytes)
--    λ> length (divisoresCompuestos5 (product [1..14]))
--    2585
--    (0.04 secs, 17,139,872 bytes)
--    λ> length (divisoresCompuestos6 (product [1..14]))
--    2585
--    (0.02 secs, 10,140,744 bytes)
--    
--    λ> length (divisoresCompuestos2 (product [1..16]))
--    5369
--    (1.97 secs, 932,433,176 bytes)
--    λ> length (divisoresCompuestos5 (product [1..16]))
--    5369
--    (0.03 secs, 37,452,088 bytes)
--    λ> length (divisoresCompuestos6 (product [1..16]))
--    5369
--    (0.03 secs, 23,017,480 bytes)
--    
--    λ> length (divisoresCompuestos5 (product [1..25]))
--    340022
--    (2.43 secs, 3,055,140,056 bytes)
--    λ> length (divisoresCompuestos6 (product [1..25]))
--    340022
--    (1.94 secs, 2,145,440,904 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

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

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

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

-- 1ª solución

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

divisoresCompuestos :: Integer -> [Integer]

divisoresCompuestos x =

[y | y <- divisores x

, y > 1

, not (isPrime y)]

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

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

divisores :: Integer -> [Integer]

divisores x =

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

, x `mod` y == 0]

-- 2ª solución

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

divisoresCompuestos2 :: Integer -> [Integer]

divisoresCompuestos2 x =

[y | y <- divisores2 x

, y > 1

, not (isPrime y)]

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

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

divisores2 :: Integer -> [Integer]

divisores2 x =

[y | y <- [1..x `div` 2], x `mod` y == 0] ++ [x]

-- 2ª solución

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

divisoresCompuestos3 :: Integer -> [Integer]

divisoresCompuestos3 x =

[y | y <- divisores2 x

, y > 1

, not (isPrime y)]

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

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

divisores3 :: Integer -> [Integer]

divisores3 x =

nub (sort (ys ++ [x `div` y | y <- ys]))

where ys = [y | y <- [1..floor (sqrt (fromIntegral x))]

, x `mod` y == 0]

-- 4ª solución

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

divisoresCompuestos4 :: Integer -> [Integer]

divisoresCompuestos4 x =

[y | y <- divisores4 x

, y > 1

, not (isPrime y)]

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

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

divisores4 :: Integer -> [Integer]

divisores4 =

nub . sort . map product . subsequences . primeFactors

-- 5ª solución

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

divisoresCompuestos5 :: Integer -> [Integer]

divisoresCompuestos5 x =

[y | y <- divisores5 x

, y > 1

, not (isPrime y)]

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

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

divisores5 :: Integer -> [Integer]

divisores5 =

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]

-- 6ª solución

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

divisoresCompuestos6 :: Integer -> [Integer]

divisoresCompuestos6 =

sort

. map product

. compuestos

. map concat

. productoCartesiano

. map inits

. group

. primeFactors

where compuestos xss = [xs | xs <- xss, length xs > 1]

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_divisoresCompuestos :: (Positive Integer) -> Bool

prop_divisoresCompuestos (Positive x) =

all (== divisoresCompuestos x) [f x | f <- [ divisoresCompuestos2

, divisoresCompuestos3

, divisoresCompuestos4

, divisoresCompuestos5

, divisoresCompuestos6 ]]

-- La comprobación es

-- λ> quickCheck prop_divisoresCompuestos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

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

-- 534

-- (14.59 secs, 7,985,108,976 bytes)

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

-- 534

-- (7.36 secs, 3,993,461,168 bytes)

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

-- 534

-- (7.35 secs, 3,993,461,336 bytes)

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

-- 534

-- (0.07 secs, 110,126,392 bytes)

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

-- 534

-- (0.01 secs, 3,332,224 bytes)

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

-- 534

-- (0.01 secs, 1,869,776 bytes)

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

-- 2585

-- (9.11 secs, 9,461,570,720 bytes)

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

-- 2585

-- (0.04 secs, 17,139,872 bytes)

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

-- 2585

-- (0.02 secs, 10,140,744 bytes)

-- λ> length (divisoresCompuestos2 (product [1..16]))

-- 5369

-- (1.97 secs, 932,433,176 bytes)

-- λ> length (divisoresCompuestos5 (product [1..16]))

-- 5369

-- (0.03 secs, 37,452,088 bytes)

-- λ> length (divisoresCompuestos6 (product [1..16]))

-- 5369

-- (0.03 secs, 23,017,480 bytes)

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

-- 340022

-- (2.43 secs, 3,055,140,056 bytes)

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

-- 340022

-- (1.94 secs, 2,145,440,904 bytes)

Pensamiento

«La verdad del hombre empieza donde acaba su propia tontería, pero la
tontería del hombre es inagotable.»

Antonio Machado

Último dígito no nulo del factorial

PorJosé A. Alonso 13-noviembre-201825-abril-2021

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 Test.QuickCheck

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

ultimoNoNuloFactorial :: Integer -> Integer
ultimoNoNuloFactorial n = ultimoNoNulo (factorial n)

-- (ultimoNoNulo n) es el último dígito no nulo de n. Por ejemplo,
--    ultimoNoNulo 5040  ==  4
ultimoNoNulo :: Integer -> Integer
ultimoNoNulo n
  | m /= 0    = m
  | otherwise = ultimoNoNulo (n `div` 10)
  where m = n `rem` 10
        
-- (factorial n) es el factorial de n. Por ejemplo,
--    factorial 7  ==  5040
factorial :: Integer -> Integer
factorial n = product [1..n]

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

ultimoNoNuloFactorial2 :: Integer -> Integer
ultimoNoNuloFactorial2 n = ultimoNoNulo2 (factorial n)

-- (ultimoNoNulo2 n) es el último dígito no nulo de n. Por ejemplo,
--    ultimoNoNulo 5040  ==  4
ultimoNoNulo2 :: Integer -> Integer
ultimoNoNulo2 n = read [head (dropWhile (=='0') (reverse (show n)))]

-- Comprobación
-- ============

-- 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 Test.QuickCheck

-- 1ª definición

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

ultimoNoNuloFactorial :: Integer -> Integer

ultimoNoNuloFactorial n = ultimoNoNulo (factorial n)

-- (ultimoNoNulo n) es el último dígito no nulo de n. Por ejemplo,

-- ultimoNoNulo 5040 == 4

ultimoNoNulo :: Integer -> Integer

ultimoNoNulo n

| m /= 0 = m

| otherwise = ultimoNoNulo (n `div` 10)

where m = n `rem` 10

-- (factorial n) es el factorial de n. Por ejemplo,

-- factorial 7 == 5040

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª definición

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

ultimoNoNuloFactorial2 :: Integer -> Integer

ultimoNoNuloFactorial2 n = ultimoNoNulo2 (factorial n)

-- (ultimoNoNulo2 n) es el último dígito no nulo de n. Por ejemplo,

-- ultimoNoNulo 5040 == 4

ultimoNoNulo2 :: Integer -> Integer

ultimoNoNulo2 n = read [head (dropWhile (=='0') (reverse (show n)))]

-- Comprobación

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

-- 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.

Pensamiento

Incierto es, lo porvenir. ¿Quién sabe lo que va a pasar? Pero incierto es también lo pretérito. ¿Quién sabe lo que ha pasado? De suerte que ni el porvenir está escrito en ninguna parte, ni el pasado tampoco.

Antonio Machado

Tren de potencias

PorJosé A. Alonso 23-julio-201824-julio-2018

Si n es el número natural cuya expansión decimal es abc… , el tren de potencias de n es a^bc^d… donde el último exponente es 1, si n tiene un número impar de dígitos. Por ejemplo

   trenDePotencias 2453  = 2^4*5^3     = 2000
   trenDePotencias 24536 = 2^4*5^3*6^1 = 12000

1 2	trenDePotencias 2453 = 2^45^3 = 2000 trenDePotencias 24536 = 2^45^3*6^1 = 12000

Definir las funciones

   trenDePotencias            :: Integer -> Integer
   esPuntoFijoTrenDePotencias :: Integer -> Bool
   puntosFijosTrenDePotencias :: [Integer]
   tablaTrenDePotencias       :: Integer -> Integer -> IO ()

trenDePotencias :: Integer -> Integer

esPuntoFijoTrenDePotencias :: Integer -> Bool

puntosFijosTrenDePotencias :: [Integer]

tablaTrenDePotencias :: Integer -> Integer -> IO ()

tales que

(trenDePotencias n) es el tren de potencia de n. Por ejemplo.

     trenDePotencias 20   ==  1
     trenDePotencias 21   ==  2
     trenDePotencias 24   ==  16
     trenDePotencias 39   ==  19683
     trenDePotencias 623  ==  108

trenDePotencias 20 == 1

trenDePotencias 21 == 2

trenDePotencias 24 == 16

trenDePotencias 39 == 19683

trenDePotencias 623 == 108

(esPuntoFijoTrenDePotencias n) se verifica si n es un punto fijo de trenDePotencias; es decir, (trenDePotencias n) es igual a n. Por ejemplo,

     esPuntoFijoTrenDePotencias 2592                        ==  True
     esPuntoFijoTrenDePotencias 24547284284866560000000000  ==  True

1 2	esPuntoFijoTrenDePotencias 2592 == True esPuntoFijoTrenDePotencias 24547284284866560000000000 == True

puntosFijosTrenDePotencias es la lista de los puntso fijos de trenDePotencias. Por ejemplo,

     take 10 puntosFijosTrenDePotencias  ==  [1,2,3,4,5,6,7,8,9,2592]

1	take 10 puntosFijosTrenDePotencias == [1,2,3,4,5,6,7,8,9,2592]

(tablaTrenDePotencias a b) es la tabla de los trenes de potencias de los números entre a y b. Por ejemplo,

     λ> tablaTrenDePotencias 20 39
     | 20 |     1 |
     | 21 |     2 |
     | 22 |     4 |
     | 23 |     8 |
     | 24 |    16 |
     | 25 |    32 |
     | 26 |    64 |
     | 27 |   128 |
     | 28 |   256 |
     | 29 |   512 |
     | 30 |     1 |
     | 31 |     3 |
     | 32 |     9 |
     | 33 |    27 |
     | 34 |    81 |
     | 35 |   243 |
     | 36 |   729 |
     | 37 |  2187 |
     | 38 |  6561 |
     | 39 | 19683 |
     λ> tablaTrenDePotencias 2340 2359
     | 2340 |        8 |
     | 2341 |       32 |
     | 2342 |      128 |
     | 2343 |      512 |
     | 2344 |     2048 |
     | 2345 |     8192 |
     | 2346 |    32768 |
     | 2347 |   131072 |
     | 2348 |   524288 |
     | 2349 |  2097152 |
     | 2350 |        8 |
     | 2351 |       40 |
     | 2352 |      200 |
     | 2353 |     1000 |
     | 2354 |     5000 |
     | 2355 |    25000 |
     | 2356 |   125000 |
     | 2357 |   625000 |
     | 2358 |  3125000 |
     | 2359 | 15625000 |

λ> tablaTrenDePotencias 20 39

| 20 | 1 |

| 21 | 2 |

| 22 | 4 |

| 23 | 8 |

| 24 | 16 |

| 25 | 32 |

| 26 | 64 |

| 27 | 128 |

| 28 | 256 |

| 29 | 512 |

| 30 | 1 |

| 31 | 3 |

| 32 | 9 |

| 33 | 27 |

| 34 | 81 |

| 35 | 243 |

| 36 | 729 |

| 37 | 2187 |

| 38 | 6561 |

| 39 | 19683 |

λ> tablaTrenDePotencias 2340 2359

| 2340 | 8 |

| 2341 | 32 |

| 2342 | 128 |

| 2343 | 512 |

| 2344 | 2048 |

| 2345 | 8192 |

| 2346 | 32768 |

| 2347 | 131072 |

| 2348 | 524288 |

| 2349 | 2097152 |

| 2350 | 8 |

| 2351 | 40 |

| 2352 | 200 |

| 2353 | 1000 |

| 2354 | 5000 |

| 2355 | 25000 |

| 2356 | 125000 |

| 2357 | 625000 |

| 2358 | 3125000 |

| 2359 | 15625000 |

Comprobar con QuickCheck que entre 2593 y 24547284284866559999999999 la función trenDePotencias no tiene puntos fijos.

Soluciones

Puedes escribir tus soluciones en los comentarios o ver las soluciones propuestas pulsando [expand title=»aquí»]

import Test.QuickCheck
import Text.Printf

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

trenDePotencias :: Integer -> Integer
trenDePotencias = trenDePotenciasLN . digitos

-- (digitos n) es la lista de los dígitos del número n. Por ejemplo,
--    digitos 2018  ==   [2,0,1,8]
digitos :: Integer -> [Integer]
digitos n =
  [read [c] | c <- show n]

-- (trenDePotenciasLN xs) es el tren de potencias de la lista de números
-- xs. Por ejemplo,
--    trenDePotenciasLN [2,4,5,3]    ==   2000
--    trenDePotenciasLN [2,4,5,3,6]  ==   12000
trenDePotenciasLN :: [Integer] -> Integer
trenDePotenciasLN []       = 1
trenDePotenciasLN [x]      = x
trenDePotenciasLN (u:v:ws) = u ^ v * (trenDePotenciasLN ws) 

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

trenDePotencias2 :: Integer -> Integer
trenDePotencias2 = trenDePotenciasLN2 . digitos

-- (trenDePotenciasLN2 xs) es el tren de potencias de la lista de números
-- xs. Por ejemplo,
--    trenDePotenciasLN2 [2,4,5,3]    ==   2000
--    trenDePotenciasLN2 [2,4,5,3,6]  ==   12000
trenDePotenciasLN2 :: [Integer] -> Integer
trenDePotenciasLN2 xs =
  product [x^y | (x,y) <- pares xs]

-- (pares xs) es la lista de los pares de elementos en la posiciones
-- pares y sus siguientes; si la longitud de xs es impar, la segunda
-- componente del último par es 1. Por ejemplo,
--    pares [2,4,5,3]    ==   [(2,4),(5,3)]
--    pares [2,4,5,3,6]  ==   [(2,4),(5,3),(6,1)]
pares :: [Integer] -> [(Integer,Integer)]
pares []       = []
pares [x]      = [(x,1)]
pares (x:y:zs) = (x,y) : pares zs

-- Equivalencia
-- ============

-- La propiedad es
prop_equivalencia :: (Positive Integer) -> Bool
prop_equivalencia (Positive n) =
  trenDePotencias n == trenDePotencias2 n

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

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

--    λ> let n = 2*10^5 in trenDePotencias (read (replicate n '2')) == 2^n
--    True
--    (2.11 secs, 2,224,301,136 bytes)
--    λ> let n = 2*10^5 in trenDePotencias2 (read (replicate n '2')) == 2^n
--    True
--    (2.08 secs, 2,237,749,216 bytes)

-- Definición de esPuntoFijoTrenDePotencias
-- ========================================

esPuntoFijoTrenDePotencias :: Integer -> Bool
esPuntoFijoTrenDePotencias n =
  trenDePotencias n == n

-- Definición de puntosFijosTrenDePotencias
-- ========================================

puntosFijosTrenDePotencias :: [Integer]
puntosFijosTrenDePotencias =
  filter esPuntoFijoTrenDePotencias [1..]

-- Definición de tablaTrenDePotencias
-- ==================================

tablaTrenDePotencias :: Integer -> Integer -> IO ()
tablaTrenDePotencias a b =
  sequence_ [printf cabecera x y | (x,y) <- trenes]
  where trenes  = [(n,trenDePotencias n) | n <- [a..b]]
        m1      = show (1 + length (show b))
        m2      = show (length (show (maximum (map snd trenes))))
        cabecera = concat ["|% ",m1,"d | %", m2,"d |\n"]

-- Comprobación
-- ============

-- La propiedad es
prop_puntosFijos :: Positive Integer -> Property
prop_puntosFijos (Positive n) =
  x < 24547284284866560000000000 ==> not (esPuntoFijoTrenDePotencias x)
  where x = 2593 + n 

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

100

101

102

103

104

105

106

107

import Test.QuickCheck

import Text.Printf

-- 1ª definición de trenDePotencias

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

trenDePotencias :: Integer -> Integer

trenDePotencias = trenDePotenciasLN . digitos

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

-- digitos 2018 == [2,0,1,8]

digitos :: Integer -> [Integer]

digitos n =

[read [c] | c <- show n]

-- (trenDePotenciasLN xs) es el tren de potencias de la lista de números

-- xs. Por ejemplo,

-- trenDePotenciasLN [2,4,5,3] == 2000

-- trenDePotenciasLN [2,4,5,3,6] == 12000

trenDePotenciasLN :: [Integer] -> Integer

trenDePotenciasLN [] = 1

trenDePotenciasLN [x] = x

trenDePotenciasLN (u:v:ws) = u ^ v * (trenDePotenciasLN ws)

-- 2ª definición de trenDePotencias

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

trenDePotencias2 :: Integer -> Integer

trenDePotencias2 = trenDePotenciasLN2 . digitos

-- (trenDePotenciasLN2 xs) es el tren de potencias de la lista de números

-- xs. Por ejemplo,

-- trenDePotenciasLN2 [2,4,5,3] == 2000

-- trenDePotenciasLN2 [2,4,5,3,6] == 12000

trenDePotenciasLN2 :: [Integer] -> Integer

trenDePotenciasLN2 xs =

product [x^y | (x,y) <- pares xs]

-- (pares xs) es la lista de los pares de elementos en la posiciones

-- pares y sus siguientes; si la longitud de xs es impar, la segunda

-- componente del último par es 1. Por ejemplo,

-- pares [2,4,5,3] == [(2,4),(5,3)]

-- pares [2,4,5,3,6] == [(2,4),(5,3),(6,1)]

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

pares [] = []

pares [x] = [(x,1)]

pares (x:y:zs) = (x,y) : pares zs

-- Equivalencia

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

-- La propiedad es

prop_equivalencia :: (Positive Integer) -> Bool

prop_equivalencia (Positive n) =

trenDePotencias n == trenDePotencias2 n

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> let n = 2*10^5 in trenDePotencias (read (replicate n '2')) == 2^n

-- True

-- (2.11 secs, 2,224,301,136 bytes)

-- λ> let n = 2*10^5 in trenDePotencias2 (read (replicate n '2')) == 2^n

-- True

-- (2.08 secs, 2,237,749,216 bytes)

-- Definición de esPuntoFijoTrenDePotencias

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

esPuntoFijoTrenDePotencias :: Integer -> Bool

esPuntoFijoTrenDePotencias n =

trenDePotencias n == n

-- Definición de puntosFijosTrenDePotencias

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

puntosFijosTrenDePotencias :: [Integer]

puntosFijosTrenDePotencias =

filter esPuntoFijoTrenDePotencias [1..]

-- Definición de tablaTrenDePotencias

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

tablaTrenDePotencias :: Integer -> Integer -> IO ()

tablaTrenDePotencias a b =

sequence_ [printf cabecera x y | (x,y) <- trenes]

where trenes = [(n,trenDePotencias n) | n <- [a..b]]

m1 = show (1 + length (show b))

m2 = show (length (show (maximum (map snd trenes))))

cabecera = concat ["|% ",m1,"d | %", m2,"d |\n"]

-- Comprobación

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

-- La propiedad es

prop_puntosFijos :: Positive Integer -> Property

prop_puntosFijos (Positive n) =

x < 24547284284866560000000000 ==> not (esPuntoFijoTrenDePotencias x)

where x = 2593 + n

-- La comprobación es

-- λ> quickCheck prop_puntosFijos

-- +++ OK, passed 100 tests.

[/expand]

El problema de las N torres

PorJosé A. Alonso 5-junio-201825-abril-2021

El problema de las N torres consiste en colocar N torres en un tablero con N filas y N columnas de forma que no haya dos torres en la misma fila ni en la misma columna.

Cada solución del problema de puede representar mediante una matriz con ceros y unos donde los unos representan las posiciones ocupadas por las torres y los ceros las posiciones libres. Por ejemplo,

   ( 0 1 0 )
   ( 1 0 0 )
   ( 0 0 1 )

( 0 1 0 )

( 1 0 0 )

( 0 0 1 )

representa una solución del problema de las 3 torres.

Definir las funciones

   torres  :: Int -> [Matrix Int]
   nTorres :: Int -> Integer

1 2	torres :: Int -> [Matrix Int] nTorres :: Int -> Integer

tales que
+ (torres n) es la lista de las soluciones del problema de las n torres. Por ejemplo,

      λ> torres 3
      [( 1 0 0 )
       ( 0 1 0 )
       ( 0 0 1 )
      ,( 1 0 0 )
       ( 0 0 1 )
       ( 0 1 0 )
      ,( 0 1 0 )
       ( 1 0 0 )
       ( 0 0 1 )
      ,( 0 1 0 )
       ( 0 0 1 )
       ( 1 0 0 )
      ,( 0 0 1 )
       ( 1 0 0 )
       ( 0 1 0 )
      ,( 0 0 1 )
       ( 0 1 0 )
       ( 1 0 0 )
      ]

λ> torres 3

[( 1 0 0 )

( 0 1 0 )

( 0 0 1 )

,( 1 0 0 )

( 0 0 1 )

( 0 1 0 )

,( 0 1 0 )

( 1 0 0 )

( 0 0 1 )

,( 0 1 0 )

( 0 0 1 )

( 1 0 0 )

,( 0 0 1 )

( 1 0 0 )

( 0 1 0 )

,( 0 0 1 )

( 0 1 0 )

( 1 0 0 )

]

(nTorres n) es el número de soluciones del problema de las n torres. Por ejemplo,

      λ> nTorres 3
      6
      λ> length (show (nTorres (10^4)))
      35660

λ> nTorres 3

λ> length (show (nTorres (10^4)))

35660

Soluciones

[schedule expon=’2018-06-12′ expat=»06:00″]

Las soluciones se pueden escribir en los comentarios hasta el 17 de abril.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

[/schedule]

[schedule on=’2018-06-12′ at=»06:00″]

import Data.List (genericLength, sort, permutations)
import Data.Matrix 

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

torres1 :: Int -> [Matrix Int]
torres1 n = 
    [permutacionAmatriz n p | p <- sort (permutations [1..n])]

permutacionAmatriz :: Int -> [Int] -> Matrix Int
permutacionAmatriz n p =
    matrix n n f
    where f (i,j) | (i,j) `elem` posiciones = 1
                  | otherwise               = 0
          posiciones = zip [1..n] p    

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

torres2 :: Int -> [Matrix Int]
torres2 = map fromLists . permutations . toLists . identity

-- El cálculo con la definición anterior es:
--    λ> identity 3
--    ( 1 0 0 )
--    ( 0 1 0 )
--    ( 0 0 1 )
--    
--    λ> toLists it
--    [[1,0,0],[0,1,0],[0,0,1]]
--    λ> permutations it
--    [[[1,0,0],[0,1,0],[0,0,1]],
--     [[0,1,0],[1,0,0],[0,0,1]],
--     [[0,0,1],[0,1,0],[1,0,0]],
--     [[0,1,0],[0,0,1],[1,0,0]],
--     [[0,0,1],[1,0,0],[0,1,0]],
--     [[1,0,0],[0,0,1],[0,1,0]]]
--    λ> map fromLists it
--    [( 1 0 0 )
--     ( 0 1 0 )
--     ( 0 0 1 )
--    ,( 0 1 0 )
--     ( 1 0 0 )
--     ( 0 0 1 )
--    ,( 0 0 1 )
--     ( 0 1 0 )
--     ( 1 0 0 )
--    ,( 0 1 0 )
--     ( 0 0 1 )
--     ( 1 0 0 )
--    ,( 0 0 1 )
--     ( 1 0 0 )
--     ( 0 1 0 )
--    ,( 1 0 0 )
--     ( 0 0 1 )
--     ( 0 1 0 )
--    ]

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

nTorres1 :: Int -> Integer
nTorres1 = genericLength . torres1

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

nTorres2 :: Int -> Integer
nTorres2 n = product [1..fromIntegral n]

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

--    λ> nTorres1 9
--    362880
--    (4.22 secs, 693,596,128 bytes)
--    λ> nTorres2 9
--    362880
--    (0.00 secs, 0 bytes)

import Data.List (genericLength, sort, permutations)

import Data.Matrix

-- 1ª definición de torres

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

torres1 :: Int -> [Matrix Int]

torres1 n =

[permutacionAmatriz n p | p <- sort (permutations [1..n])]

permutacionAmatriz :: Int -> [Int] -> Matrix Int

permutacionAmatriz n p =

matrix n n f

where f (i,j) | (i,j) `elem` posiciones = 1

| otherwise = 0

posiciones = zip [1..n] p

-- 2ª definición de torres

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

torres2 :: Int -> [Matrix Int]

torres2 = map fromLists . permutations . toLists . identity

-- El cálculo con la definición anterior es:

-- λ> identity 3

-- ( 1 0 0 )

-- ( 0 1 0 )

-- ( 0 0 1 )

-- λ> toLists it

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

-- λ> permutations it

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

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

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

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

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

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

-- λ> map fromLists it

-- [( 1 0 0 )

-- ( 0 1 0 )

-- ( 0 0 1 )

-- ,( 0 1 0 )

-- ( 1 0 0 )

-- ( 0 0 1 )

-- ,( 0 0 1 )

-- ( 0 1 0 )

-- ( 1 0 0 )

-- ,( 0 1 0 )

-- ( 0 0 1 )

-- ( 1 0 0 )

-- ,( 0 0 1 )

-- ( 1 0 0 )

-- ( 0 1 0 )

-- ,( 1 0 0 )

-- ( 0 0 1 )

-- ( 0 1 0 )

-- ]

-- 1ª definición de nTorres

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

nTorres1 :: Int -> Integer

nTorres1 = genericLength . torres1

-- 2ª definición de nTorres

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

nTorres2 :: Int -> Integer

nTorres2 n = product [1..fromIntegral n]

-- Comparación de eficiencia

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

-- λ> nTorres1 9

-- 362880

-- (4.22 secs, 693,596,128 bytes)

-- λ> nTorres2 9

-- 362880

-- (0.00 secs, 0 bytes)

[/schedule]

Número de triangulaciones de un polígono

PorJosé A. Alonso 7-mayo-201825-abril-2021

Una triangulación de un polígono es una división del área en un conjunto de triángulos, de forma que la unión de todos ellos es igual al polígono original, y cualquier par de triángulos es disjunto o comparte únicamente un vértice o un lado. En el caso de polígonos convexos, la cantidad de triangulaciones posibles depende únicamente del número de vértices del polígono.

Si llamamos T(n) al número de triangulaciones de un polígono de n vértices, se verifica la siguiente relación de recurrencia:

    T(2) = 1
    T(n) = T(2)*T(n-1) + T(3)*T(n-2) + ... + T(n-1)*T(2)

1 2	T(2) = 1 T(n) = T(2)T(n-1) + T(3)T(n-2) + ... + T(n-1)*T(2)

Definir la función

   numeroTriangulaciones :: Integer -> Integer

1	numeroTriangulaciones :: Integer -> Integer

tal que (numeroTriangulaciones n) es el número de triangulaciones de un polígono convexo de n vértices. Por ejemplo,

   numeroTriangulaciones 3  == 1
   numeroTriangulaciones 5  == 5
   numeroTriangulaciones 6  == 14
   numeroTriangulaciones 7  == 42
   numeroTriangulaciones 50 == 131327898242169365477991900
   length (show (numeroTriangulaciones   800)) ==  476
   length (show (numeroTriangulaciones  1000)) ==  597
   length (show (numeroTriangulaciones 10000)) == 6014

numeroTriangulaciones 3 == 1

numeroTriangulaciones 5 == 5

numeroTriangulaciones 6 == 14

numeroTriangulaciones 7 == 42

numeroTriangulaciones 50 == 131327898242169365477991900

length (show (numeroTriangulaciones 800)) == 476

length (show (numeroTriangulaciones 1000)) == 597

length (show (numeroTriangulaciones 10000)) == 6014

Soluciones

import Data.Array (Array, (!), array)
import Data.List  (genericIndex)
import qualified Data.Vector as V

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

numeroTriangulaciones :: Integer -> Integer
numeroTriangulaciones 2 = 1
numeroTriangulaciones n = sum (zipWith (*) ts (reverse ts))
  where ts = [numeroTriangulaciones k | k <- [2..n-1]]

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

--    λ> map numeroTriangulaciones2 [2..15]
--    [1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900]
numeroTriangulaciones2 :: Integer -> Integer
numeroTriangulaciones2 n = 
  head (sucNumeroTriangulacionesInversas `genericIndex` (n-2))

--    λ> mapM_ print (take 10 sucNumeroTriangulacionesInversas)
--    [1]
--    [1,1]
--    [2,1,1]
--    [5,2,1,1]
--    [14,5,2,1,1]
--    [42,14,5,2,1,1]
--    [132,42,14,5,2,1,1]
--    [429,132,42,14,5,2,1,1]
--    [1430,429,132,42,14,5,2,1,1]
--    [4862,1430,429,132,42,14,5,2,1,1]
sucNumeroTriangulacionesInversas :: [[Integer]]        
sucNumeroTriangulacionesInversas = iterate f [1]
  where f ts = sum (zipWith (*) ts (reverse ts)) : ts

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

numeroTriangulaciones3 :: Integer -> Integer
numeroTriangulaciones3 n = vectorTriang n ! n

--    λ> vectorTriang 9
--    array (2,9) [(2,1),(3,1),(4,2),(5,5),(6,14),(7,42),(8,132),(9,429)]
vectorTriang :: Integer -> Array Integer Integer
vectorTriang n = v
  where v = array (2,n) [(i, f i) | i <-[2..n]]
        f 2 = 1
        f i = sum [v!j*v!(i-j+1) | j <-[2..i-1]]

-- 4ª solución (con los números de Catalan)
-- ========================================

-- Al calcular por primeros números de triangulaciones se obtiene
--    λ> map numeroTriangulaciones [2..12]
--    [1,1,2,5,14,42,132,429,1430,4862,16796]
-- Se observa que se corresponden con los números de Catalan
-- http://bit.ly/2FOc1S1
-- 
-- El n-ésimo número de Catalan es
--    (2n)! / (n! * (n+1)!)
-- El número de triangulaciones de un polígono de n lados es el
-- (n-2)-ésimo número de Catalan.

numeroTriangulaciones4 :: Integer -> Integer
numeroTriangulaciones4 n = numeroCatalan (n-2) 

numeroCatalan :: Integer -> Integer
numeroCatalan n =
  factorial (2*n) `div` (factorial (n+1) * factorial n)

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

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

numeroTriangulaciones5 :: Integer -> Integer
numeroTriangulaciones5 n = v V.! (m-2)
   where v = V.constructN m
           (\w -> if   V.null w then 1
                  else V.sum (V.zipWith (*) w (V.reverse w)))
         m = fromIntegral n


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

numeroTriangulaciones6 :: Integer -> Integer
numeroTriangulaciones6 n = nCr (2*n-4, n-2) `div` (n-1)
  where nCr (n,m)   = factorial n `div` (factorial (n-m) * factorial m)
        factorial n = product [1..n]

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

--    λ> numeroTriangulaciones 22
--    6564120420
--    (3.97 secs, 668,070,936 bytes)
--    λ> numeroTriangulaciones2 22
--    6564120420
--    (0.01 secs, 180,064 bytes)
--    λ> numeroTriangulaciones3 22
--    6564120420
--    (0.01 secs, 285,792 bytes)
--    
--    λ> length (show (numeroTriangulaciones2 800))
--    476
--    (0.59 secs, 125,026,824 bytes)
--    λ> length (show (numeroTriangulaciones3 800))
--    476
--    (1.95 secs, 334,652,936 bytes)
--    λ> length (show (numeroTriangulaciones4 800))
--    476
--    (0.01 secs, 2,960,088 bytes)
--    λ> length (show (numeroTriangulaciones5 800))
--    476
--    (0.65 secs, 200,415,640 bytes)
--    λ> length (show (numeroTriangulaciones6 800))
--    476
--    (0.01 secs, 2,960,224 bytes)
--    
--    λ> length (show (numeroTriangulaciones4 (10^4)))
--    6014
--    (1.80 secs, 542,364,320 bytes)
--    λ> length (show (numeroTriangulaciones6 (10^4)))
--    6014
--    (1.87 secs, 542,351,136 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

import Data.Array (Array, (!), array)

import Data.List (genericIndex)

import qualified Data.Vector as V

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

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

numeroTriangulaciones :: Integer -> Integer

numeroTriangulaciones 2 = 1

numeroTriangulaciones n = sum (zipWith (*) ts (reverse ts))

where ts = [numeroTriangulaciones k | k <- [2..n-1]]

-- 2ª solución

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

-- λ> map numeroTriangulaciones2 [2..15]

-- [1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900]

numeroTriangulaciones2 :: Integer -> Integer

numeroTriangulaciones2 n =

head (sucNumeroTriangulacionesInversas `genericIndex` (n-2))

-- λ> mapM_ print (take 10 sucNumeroTriangulacionesInversas)

-- [1]

-- [1,1]

-- [2,1,1]

-- [5,2,1,1]

-- [14,5,2,1,1]

-- [42,14,5,2,1,1]

-- [132,42,14,5,2,1,1]

-- [429,132,42,14,5,2,1,1]

-- [1430,429,132,42,14,5,2,1,1]

-- [4862,1430,429,132,42,14,5,2,1,1]

sucNumeroTriangulacionesInversas :: [[Integer]]

sucNumeroTriangulacionesInversas = iterate f [1]

where f ts = sum (zipWith (*) ts (reverse ts)) : ts

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

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

numeroTriangulaciones3 :: Integer -> Integer

numeroTriangulaciones3 n = vectorTriang n ! n

-- λ> vectorTriang 9

-- array (2,9) [(2,1),(3,1),(4,2),(5,5),(6,14),(7,42),(8,132),(9,429)]

vectorTriang :: Integer -> Array Integer Integer

vectorTriang n = v

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

f 2 = 1

f i = sum [v!j*v!(i-j+1) | j <-[2..i-1]]

-- 4ª solución (con los números de Catalan)

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

-- Al calcular por primeros números de triangulaciones se obtiene

-- λ> map numeroTriangulaciones [2..12]

-- [1,1,2,5,14,42,132,429,1430,4862,16796]

-- Se observa que se corresponden con los números de Catalan

-- http://bit.ly/2FOc1S1

-- El n-ésimo número de Catalan es

-- (2n)! / (n! * (n+1)!)

-- El número de triangulaciones de un polígono de n lados es el

-- (n-2)-ésimo número de Catalan.

numeroTriangulaciones4 :: Integer -> Integer

numeroTriangulaciones4 n = numeroCatalan (n-2)

numeroCatalan :: Integer -> Integer

numeroCatalan n =

factorial (2*n) `div` (factorial (n+1) * factorial n)

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 5ª solución

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

numeroTriangulaciones5 :: Integer -> Integer

numeroTriangulaciones5 n = v V.! (m-2)

where v = V.constructN m

(\w -> if V.null w then 1

else V.sum (V.zipWith (*) w (V.reverse w)))

m = fromIntegral n

-- 6ª solución

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

numeroTriangulaciones6 :: Integer -> Integer

numeroTriangulaciones6 n = nCr (2*n-4, n-2) `div` (n-1)

where nCr (n,m) = factorial n `div` (factorial (n-m) * factorial m)

factorial n = product [1..n]

-- Comparación de eficiencia

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

-- λ> numeroTriangulaciones 22

-- 6564120420

-- (3.97 secs, 668,070,936 bytes)

-- λ> numeroTriangulaciones2 22

-- 6564120420

-- (0.01 secs, 180,064 bytes)

-- λ> numeroTriangulaciones3 22

-- 6564120420

-- (0.01 secs, 285,792 bytes)

-- λ> length (show (numeroTriangulaciones2 800))

-- 476

-- (0.59 secs, 125,026,824 bytes)

-- λ> length (show (numeroTriangulaciones3 800))

-- 476

-- (1.95 secs, 334,652,936 bytes)

-- λ> length (show (numeroTriangulaciones4 800))

-- 476

-- (0.01 secs, 2,960,088 bytes)

-- λ> length (show (numeroTriangulaciones5 800))

-- 476

-- (0.65 secs, 200,415,640 bytes)

-- λ> length (show (numeroTriangulaciones6 800))

-- 476

-- (0.01 secs, 2,960,224 bytes)

-- λ> length (show (numeroTriangulaciones4 (10^4)))

-- 6014

-- (1.80 secs, 542,364,320 bytes)

-- λ> length (show (numeroTriangulaciones6 (10^4)))

-- 6014

-- (1.87 secs, 542,351,136 bytes)

La sucesión de Sylvester

PorJosé A. Alonso 19-marzo-201825-abril-2021

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

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
-- ===========

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)

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

sylvester5 :: Integer -> Integer
sylvester5 0 = 2
sylvester5 n = 1 + (productosSylvester `genericIndex` (n-1))

sucSylvester5 :: [Integer]
sucSylvester5 = map sylvester5 [0..]

productosSylvester :: [Integer]
productosSylvester = scanl1 (*) sucSylvester5

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

--    λ> length (show (sylvester1 20))
--    213441
--    (3.40 secs, 519,249,840 bytes)
--    λ> length (show (sylvester2 20))
--    213441
--    (0.10 secs, 13,716,024 bytes)
--    λ> length (show (sylvester3 20))
--    213441
--    (0.16 secs, 13,646,472 bytes)
--    λ> length (show (sylvester4 20))
--    213441
--    (0.18 secs, 13,654,064 bytes)
--    λ> length (show (sylvester5 20))
--    213441
--    (0.12 secs, 13,497,192 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

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

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)

-- 4ª solución

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

sylvester5 :: Integer -> Integer

sylvester5 0 = 2

sylvester5 n = 1 + (productosSylvester `genericIndex` (n-1))

sucSylvester5 :: [Integer]

sucSylvester5 = map sylvester5 [0..]

productosSylvester :: [Integer]

productosSylvester = scanl1 (*) sucSylvester5

-- Comparación de eficiencia

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

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

-- 213441

-- (3.40 secs, 519,249,840 bytes)

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

-- 213441

-- (0.10 secs, 13,716,024 bytes)

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

-- 213441

-- (0.16 secs, 13,646,472 bytes)

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

-- 213441

-- (0.18 secs, 13,654,064 bytes)

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

-- 213441

-- (0.12 secs, 13,497,192 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]]

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Nuevas soluciones

Soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico