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

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.

Conjunto de primos relativos

PorJosé A. Alonso 30-marzo-20226-abril-2022

Dos números enteros positivos son primos relativos si no tienen ningún factor primo en común; es decit, si 1 es su único divisor común. Por ejemplo, 6 y 35 son primos entre sí, pero 6 y 27 no lo son porque ambos son divisibles por 3.

Definir la función

   primosRelativos :: [Int] -> Bool

1	primosRelativos :: [Int] -> Bool

tal que (primosRelativos xs) se verifica si los elementos de xs son primos relativos dos a dos. Por ejemplo,

   primosRelativos [6,35]         ==  True
   primosRelativos [6,27]         ==  False
   primosRelativos [2,3,4]        ==  False
   primosRelativos [6,35,11]      ==  True
   primosRelativos [6,35,11,221]  ==  True
   primosRelativos [6,35,11,231]  ==  False

primosRelativos [6,35] == True

primosRelativos [6,27] == False

primosRelativos [2,3,4] == False

primosRelativos [6,35,11] == True

primosRelativos [6,35,11,221] == True

primosRelativos [6,35,11,231] == False

Soluciones

import Test.QuickCheck
import Data.List (delete, intersect)
import Data.Numbers.Primes (primeFactors, primes)
import qualified Data.Set as S (disjoint, fromList)

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

primosRelativos1 :: [Int] -> Bool
primosRelativos1 []     = True
primosRelativos1 (x:xs) =
  and [sonPrimosRelativos1 x y | y <- xs] && primosRelativos1 xs

-- (sonPrimosRelativos x y) se verifica si x e y son primos
-- relativos. Por ejemplo,
--    sonPrimosRelativos1 6 35  ==  True
--    sonPrimosRelativos1 6 27  ==  False
sonPrimosRelativos1 :: Int -> Int -> Bool
sonPrimosRelativos1 x y =
  null (divisoresPrimos x `intersect` divisoresPrimos y)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por
-- ejemplo,
--    divisoresPrimos 600  ==  [2,2,2,3,5,5]
divisoresPrimos :: Int -> [Int]
divisoresPrimos 1 = []
divisoresPrimos x =
  y : divisoresPrimos (x `div` y)
  where y = menorDivisorPrimo x

-- (menorDivisorPrimo x) es el menor divisor primo de x. Por ejemplo,
--    menorDivisorPrimo 15  ==  3
--    menorDivisorPrimo 11  ==  11
menorDivisorPrimo :: Int -> Int
menorDivisorPrimo x =
  head [y | y <- [2..], x `mod` y == 0]

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

primosRelativos2 :: [Int] -> Bool
primosRelativos2 []     = True
primosRelativos2 (x:xs) =
  all (sonPrimosRelativos1 x) xs && primosRelativos2 xs

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

primosRelativos3 :: [Int] -> Bool
primosRelativos3 []     = True
primosRelativos3 (x:xs) =
  all (sonPrimosRelativos2 x) xs && primosRelativos3 xs

sonPrimosRelativos2 :: Int -> Int -> Bool
sonPrimosRelativos2 x y =
  null (primeFactors x `intersect` primeFactors y)

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

primosRelativos4 :: [Int] -> Bool
primosRelativos4 []     = True
primosRelativos4 (x:xs) =
  all (sonPrimosRelativos3 x) xs && primosRelativos4 xs

sonPrimosRelativos3 :: Int -> Int -> Bool
sonPrimosRelativos3 x y =
  S.fromList (primeFactors x) `S.disjoint` S.fromList (primeFactors y)

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

primosRelativos5 :: [Int] -> Bool
primosRelativos5 []     = True
primosRelativos5 (x:xs) =
  all (sonPrimosRelativos5 x) xs && primosRelativos5 xs

sonPrimosRelativos5 :: Int -> Int -> Bool
sonPrimosRelativos5 x y =
  gcd x y == 1

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

-- La propiedad es
prop_primosRelativos :: [Positive Int] -> Bool
prop_primosRelativos xs =
  all (== primosRelativos1 ys)
      [primosRelativos2 ys,
       primosRelativos3 ys,
       primosRelativos4 ys,
       primosRelativos5 ys]
  where ys = getPositive <$> xs

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

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

-- La comparación es
--    λ> primosRelativos1 (take 120 primes)
--    True
--    (1.92 secs, 869,909,416 bytes)
--    λ> primosRelativos2 (take 120 primes)
--    True
--    (1.99 secs, 869,045,656 bytes)
--    λ> primosRelativos3 (take 120 primes)
--    True
--    (0.09 secs, 221,183,200 bytes)
--
--    λ> primosRelativos3 (take 600 primes)
--    True
--    (2.62 secs, 11,196,690,856 bytes)
--    λ> primosRelativos4 (take 600 primes)
--    True
--    (2.66 secs, 11,190,940,456 bytes)
--    λ> primosRelativos5 (take 600 primes)
--    True
--    (0.14 secs, 123,673,648 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

import Test.QuickCheck

import Data.List (delete, intersect)

import Data.Numbers.Primes (primeFactors, primes)

import qualified Data.Set as S (disjoint, fromList)

-- 1ª solución

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

primosRelativos1 :: [Int] -> Bool

primosRelativos1 [] = True

primosRelativos1 (x:xs) =

and [sonPrimosRelativos1 x y | y <- xs] && primosRelativos1 xs

-- (sonPrimosRelativos x y) se verifica si x e y son primos

-- relativos. Por ejemplo,

-- sonPrimosRelativos1 6 35 == True

-- sonPrimosRelativos1 6 27 == False

sonPrimosRelativos1 :: Int -> Int -> Bool

sonPrimosRelativos1 x y =

null (divisoresPrimos x `intersect` divisoresPrimos y)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por

-- ejemplo,

-- divisoresPrimos 600 == [2,2,2,3,5,5]

divisoresPrimos :: Int -> [Int]

divisoresPrimos 1 = []

divisoresPrimos x =

y : divisoresPrimos (x `div` y)

where y = menorDivisorPrimo x

-- (menorDivisorPrimo x) es el menor divisor primo de x. Por ejemplo,

-- menorDivisorPrimo 15 == 3

-- menorDivisorPrimo 11 == 11

menorDivisorPrimo :: Int -> Int

menorDivisorPrimo x =

head [y | y <- [2..], x `mod` y == 0]

-- 2ª solución

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

primosRelativos2 :: [Int] -> Bool

primosRelativos2 [] = True

primosRelativos2 (x:xs) =

all (sonPrimosRelativos1 x) xs && primosRelativos2 xs

-- 3ª solución

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

primosRelativos3 :: [Int] -> Bool

primosRelativos3 [] = True

primosRelativos3 (x:xs) =

all (sonPrimosRelativos2 x) xs && primosRelativos3 xs

sonPrimosRelativos2 :: Int -> Int -> Bool

sonPrimosRelativos2 x y =

null (primeFactors x `intersect` primeFactors y)

-- 4ª solución

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

primosRelativos4 :: [Int] -> Bool

primosRelativos4 [] = True

primosRelativos4 (x:xs) =

all (sonPrimosRelativos3 x) xs && primosRelativos4 xs

sonPrimosRelativos3 :: Int -> Int -> Bool

sonPrimosRelativos3 x y =

S.fromList (primeFactors x) `S.disjoint` S.fromList (primeFactors y)

-- 5ª solución

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

primosRelativos5 :: [Int] -> Bool

primosRelativos5 [] = True

primosRelativos5 (x:xs) =

all (sonPrimosRelativos5 x) xs && primosRelativos5 xs

sonPrimosRelativos5 :: Int -> Int -> Bool

sonPrimosRelativos5 x y =

gcd x y == 1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosRelativos :: [Positive Int] -> Bool

prop_primosRelativos xs =

all (== primosRelativos1 ys)

[primosRelativos2 ys,

primosRelativos3 ys,

primosRelativos4 ys,

primosRelativos5 ys]

where ys = getPositive <$> xs

-- La comprobación es

-- λ> quickCheck prop_primosRelativos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosRelativos1 (take 120 primes)

-- True

-- (1.92 secs, 869,909,416 bytes)

-- λ> primosRelativos2 (take 120 primes)

-- True

-- (1.99 secs, 869,045,656 bytes)

-- λ> primosRelativos3 (take 120 primes)

-- True

-- (0.09 secs, 221,183,200 bytes)

-- λ> primosRelativos3 (take 600 primes)

-- True

-- (2.62 secs, 11,196,690,856 bytes)

-- λ> primosRelativos4 (take 600 primes)

-- True

-- (2.66 secs, 11,190,940,456 bytes)

-- λ> primosRelativos5 (take 600 primes)

-- True

-- (0.14 secs, 123,673,648 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>

Números de Perrin

PorJosé A. Alonso 19-mayo-202026-mayo-2020

Los números de Perrin se definen por la elación de recurrencia

   P(n) = P(n - 2) + P(n - 3) si n > 2,

1	P(n) = P(n - 2) + P(n - 3) si n > 2,

con los valores iniciales

   P(0) = 3, P(1) = 0 y P(2) = 2.

1	P(0) = 3, P(1) = 0 y P(2) = 2.

Definir la sucesión

   sucPerrin :: [Integer]

1	sucPerrin :: [Integer]

cuyos elementos son los números de Perrin. Por ejemplo,

   λ> take 15 sucPerrin
   [3,0,2,3,2,5,5,7,10,12,17,22,29,39,51]
   λ> length (show (sucPerrin !! (2*10^5)))
   24425

λ> take 15 sucPerrin

[3,0,2,3,2,5,5,7,10,12,17,22,29,39,51]

λ> length (show (sucPerrin !! (2*10^5)))

24425

Comprobar con QuickCheck si se verifica la siguiente propiedad: para todo entero n > 1, el n-ésimo término de la sucesión de Perrin es divisible por n si y sólo si n es primo.

Soluciones

import Data.List (genericIndex, unfoldr)
import Data.Numbers.Primes (isPrime)
import Test.QuickCheck

-- 1ª solución
sucPerrin1 :: [Integer]
sucPerrin1 = 3 : 0 : 2 : zipWith (+) sucPerrin1 (tail sucPerrin1)

-- 2ª solución
sucPerrin2 :: [Integer]
sucPerrin2 = [x | (x,_,_) <- iterate op (3,0,2)]
  where op (a,b,c) = (b,c,a+b)
 
-- 3ª solución
sucPerrin3 :: [Integer]
sucPerrin3 =
  unfoldr (\(a, (b,c)) -> Just (a, (b,(c,a+b)))) (3,(0,2))

-- 4ª solución
sucPerrin4 :: [Integer]
sucPerrin4 = [vectorPerrin n ! n | n <- [0..]]

vectorPerrin :: Integer -> Array Integer Integer
vectorPerrin n = v where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 3
  f 1 = 0
  f 2 = 2
  f i = v ! (i-2) + v ! (i-3)

-- Comparación de eficiencia
--    λ> length (show (sucPerrin1 !! (3*10^5)))
--    36638
--    (1.62 secs, 2,366,238,984 bytes)
--    λ> length (show (sucPerrin2 !! (3*10^5)))
--    36638
--    (1.40 secs, 2,428,701,384 bytes)
--    λ> length (show (sucPerrin3 !! (3*10^5)))
--    36638
--    (1.14 secs, 2,409,504,864 bytes)
--    λ> length (show (sucPerrin4 !! (3*10^5)))
--    36638
--    (1.78 secs, 2,585,400,776 bytes)


-- Usaremos la 3ª
sucPerrin :: [Integer]
sucPerrin = sucPerrin3

-- La propiedad es  
conjeturaPerrin :: Integer -> Property
conjeturaPerrin n =
  n > 1 ==>
  (perrin n `mod` n == 0) == isPrime n

-- (perrin n) es el n-ésimo término de la sucesión de Perrin. Por
-- ejemplo,
--    perrin 4  ==  2
--    perrin 5  ==  5
--    perrin 6  ==  5
perrin :: Integer -> Integer
perrin n = sucPerrin `genericIndex` n

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

-- Nota: Aunque QuickCheck no haya encontrado contraejemplos, la
-- propiedad no es cierta. Sólo lo es una de las implicaciones: si n es
-- primo, entonces el  n-ésimo término de la sucesión de Perrin es
-- divisible por n. La otra es falsa y los primeros contraejemplos son
--    271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291

import Data.List (genericIndex, unfoldr)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª solución

sucPerrin1 :: [Integer]

sucPerrin1 = 3 : 0 : 2 : zipWith (+) sucPerrin1 (tail sucPerrin1)

-- 2ª solución

sucPerrin2 :: [Integer]

sucPerrin2 = [x | (x,_,_) <- iterate op (3,0,2)]

where op (a,b,c) = (b,c,a+b)

-- 3ª solución

sucPerrin3 :: [Integer]

sucPerrin3 =

unfoldr (\(a, (b,c)) -> Just (a, (b,(c,a+b)))) (3,(0,2))

-- 4ª solución

sucPerrin4 :: [Integer]

sucPerrin4 = [vectorPerrin n ! n | n <- [0..]]

vectorPerrin :: Integer -> Array Integer Integer

vectorPerrin n = v where

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

f 0 = 3

f 1 = 0

f 2 = 2

f i = v ! (i-2) + v ! (i-3)

-- Comparación de eficiencia

-- λ> length (show (sucPerrin1 !! (3*10^5)))

-- 36638

-- (1.62 secs, 2,366,238,984 bytes)

-- λ> length (show (sucPerrin2 !! (3*10^5)))

-- 36638

-- (1.40 secs, 2,428,701,384 bytes)

-- λ> length (show (sucPerrin3 !! (3*10^5)))

-- 36638

-- (1.14 secs, 2,409,504,864 bytes)

-- λ> length (show (sucPerrin4 !! (3*10^5)))

-- 36638

-- (1.78 secs, 2,585,400,776 bytes)

-- Usaremos la 3ª

sucPerrin :: [Integer]

sucPerrin = sucPerrin3

-- La propiedad es

conjeturaPerrin :: Integer -> Property

conjeturaPerrin n =

n > 1 ==>

(perrin n `mod` n == 0) == isPrime n

-- (perrin n) es el n-ésimo término de la sucesión de Perrin. Por

-- ejemplo,

-- perrin 4 == 2

-- perrin 5 == 5

-- perrin 6 == 5

perrin :: Integer -> Integer

perrin n = sucPerrin `genericIndex` n

-- La comprobación es

-- λ> quickCheck conjeturaPerrin

-- +++ OK, passed 100 tests.

-- Nota: Aunque QuickCheck no haya encontrado contraejemplos, la

-- propiedad no es cierta. Sólo lo es una de las implicaciones: si n es

-- primo, entonces el n-ésimo término de la sucesión de Perrin es

-- divisible por n. La otra es falsa y los primeros contraejemplos son

-- 271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291

Cadenas de divisores

PorJosé A. Alonso 14-mayo-202021-mayo-2020

Una cadena de divisores de un número n es una lista donde cada elemento es un divisor de su siguiente elemento en la lista. Por ejemplo, las cadenas de divisores de 12 son [2,4,12], [2,6,12], [2,12], [3,6,12], [3,12], [4,12], [6,12] y [12].

Definir la función

   cadenasDivisores :: Int -> [[Int]]

1	cadenasDivisores :: Int -> [[Int]]

tal que (cadenasDivisores n) es la lista de las cadenas de divisores de n. Por ejemplo,

   λ> cadenasDivisores 12
   [[2,4,12],[2,6,12],[2,12],[3,6,12],[3,12],[4,12],[6,12],[12]]
   λ> length (cadenaDivisores 48)
   48
   λ> length (cadenaDivisores 120)
   132

λ> cadenasDivisores 12

[[2,4,12],[2,6,12],[2,12],[3,6,12],[3,12],[4,12],[6,12],[12]]

λ> length (cadenaDivisores 48)

λ> length (cadenaDivisores 120)

132

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (isPrime)

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

cadenasDivisores :: Int -> [[Int]]
cadenasDivisores n = sort (extiendeLista [[n]])
    where extiendeLista []           = []
          extiendeLista ((1:xs):yss) = xs : extiendeLista yss
          extiendeLista ((x:xs):yss) =
              extiendeLista ([y:x:xs | y <- divisores x] ++ yss)

-- (divisores x) es la lista decreciente de los divisores de x distintos
-- de x. Por ejemplo,
--    divisores 12  ==  [6,4,3,2,1]
divisores :: Int -> [Int]
divisores x = 
    [y | y <- [a,a-1..1], x `mod` y == 0]
    where a = x `div` 2

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

cadenasDivisores2 :: Int -> [[Int]]
cadenasDivisores2 = sort . aux
    where aux 1 = [[]]
          aux n = [xs ++ [n] | xs <- concatMap aux (divisores n)]

-- 3ª definición
-- =============

cadenasDivisores3 :: Int -> [[Int]]
cadenasDivisores3 = sort . map reverse . aux
    where aux 1 = [[]]
          aux n = map (n:) (concatMap aux (divisores3 n))

-- (divisores3 x) es la lista creciente de los divisores de x distintos
-- de x. Por ejemplo,
--    divisores3 12  ==  [1,2,3,4,6]
divisores3 :: Int -> [Int]
divisores3 x = 
    [y | y <- [1..a], x `mod` y == 0]
    where a = x `div` 2

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

nCadenasDivisores1 :: Int -> Int
nCadenasDivisores1 = length . cadenasDivisores

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

nCadenasDivisores2 :: Int -> Int
nCadenasDivisores2 1 = 1
nCadenasDivisores2 n = 
    sum [nCadenasDivisores2 x | x <- divisores n]

import Data.List (sort)

import Data.Numbers.Primes (isPrime)

-- 1ª definición

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

cadenasDivisores :: Int -> [[Int]]

cadenasDivisores n = sort (extiendeLista [[n]])

where extiendeLista [] = []

extiendeLista ((1:xs):yss) = xs : extiendeLista yss

extiendeLista ((x:xs):yss) =

extiendeLista ([y:x:xs | y <- divisores x] ++ yss)

-- (divisores x) es la lista decreciente de los divisores de x distintos

-- de x. Por ejemplo,

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

divisores :: Int -> [Int]

divisores x =

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

where a = x `div` 2

-- 2ª definición

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

cadenasDivisores2 :: Int -> [[Int]]

cadenasDivisores2 = sort . aux

where aux 1 = [[]]

aux n = [xs ++ [n] | xs <- concatMap aux (divisores n)]

-- 3ª definición

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

cadenasDivisores3 :: Int -> [[Int]]

cadenasDivisores3 = sort . map reverse . aux

where aux 1 = [[]]

aux n = map (n:) (concatMap aux (divisores3 n))

-- (divisores3 x) es la lista creciente de los divisores de x distintos

-- de x. Por ejemplo,

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

divisores3 :: Int -> [Int]

divisores3 x =

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

where a = x `div` 2

-- 1ª definición de nCadenasDivisores

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

nCadenasDivisores1 :: Int -> Int

nCadenasDivisores1 = length . cadenasDivisores

-- 2ª definición de nCadenasDivisores

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

nCadenasDivisores2 :: Int -> Int

nCadenasDivisores2 1 = 1

nCadenasDivisores2 n =

sum [nCadenasDivisores2 x | x <- divisores n]

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

Período de una lista

PorJosé A. Alonso 30-abril-20206-mayo-2020

El período de una lista xs es la lista más corta ys tal que xs se puede obtener concatenando varias veces la lista ys. Por ejemplo, el período «abababab» es «ab» ya que «abababab» se obtiene repitiendo tres veces la lista «ab».

Definir la función

   periodo :: Eq a => [a] -> [a]

1	periodo :: Eq a => [a] -> [a]

tal que (periodo xs) es el período de xs. Por ejemplo,

   periodo "ababab"      ==  "ab"
   periodo "buenobueno"  ==  "bueno"
   periodo "oooooo"      ==  "o"
   periodo "sevilla"     ==  "sevilla"

periodo "ababab" == "ab"

periodo "buenobueno" == "bueno"

periodo "oooooo" == "o"

periodo "sevilla" == "sevilla"

Soluciones

import Data.List (isPrefixOf, inits)

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

periodo1 :: Eq a => [a] -> [a]
periodo1 xs = take n xs
    where l = length xs
          n = head [m | m <- divisores l, 
                        concat (replicate (l `div` m) (take m xs)) == xs]

-- (divisores n) es la lista de los divisores de n. Por ejemplo,
--    divisores 96  ==  [1,2,3,4,6,8,12,16,24,32,48,96]
divisores :: Int -> [Int]
divisores n = [x | x <- [1..n], n `mod` x == 0]

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

periodo2 :: Eq a => [a] -> [a]
periodo2 xs = take n xs
    where l = length xs
          n = head [m | m <- divisores l, 
                        xs `isPrefixOf` cycle (take m xs)]

import Data.List (isPrefixOf, inits)

-- 1ª solución

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

periodo1 :: Eq a => [a] -> [a]

periodo1 xs = take n xs

where l = length xs

n = head [m | m <- divisores l,

concat (replicate (l `div` m) (take m xs)) == xs]

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

-- divisores 96 == [1,2,3,4,6,8,12,16,24,32,48,96]

divisores :: Int -> [Int]

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

-- 2ª solución

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

periodo2 :: Eq a => [a] -> [a]

periodo2 xs = take n xs

where l = length xs

n = head [m | m <- divisores l,

xs `isPrefixOf` cycle (take m xs)]

Mayor capicúa producto de dos números de n cifras

PorJosé A. Alonso 29-abril-202012-mayo-2020

Un capicúa es un número que es igual leído de izquierda a derecha que de derecha a izquierda.

Definir la función

   mayorCapicuaP :: Integer -> Integer

1	mayorCapicuaP :: Integer -> Integer

tal que (mayorCapicuaP n) es el mayor capicúa que es el producto de dos números de n cifras. Por ejemplo,

   mayorCapicuaP 2  ==  9009
   mayorCapicuaP 3  ==  906609
   mayorCapicuaP 4  ==  99000099
   mayorCapicuaP 5  ==  9966006699
   mayorCapicuaP 6  ==  999000000999
   mayorCapicuaP 7  ==  99956644665999

mayorCapicuaP 2 == 9009

mayorCapicuaP 3 == 906609

mayorCapicuaP 4 == 99000099

mayorCapicuaP 5 == 9966006699

mayorCapicuaP 6 == 999000000999

mayorCapicuaP 7 == 99956644665999

Soluciones

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

mayorCapicuaP1 :: Integer -> Integer
mayorCapicuaP1 n = head (capicuasP n)

-- (capicuasP n) es la lista de las capicúas de 2*n cifras que
-- pueden escribirse como productos de dos números de n cifras. Por
-- ejemplo, Por ejemplo,
--    ghci> capicuasP 2
--    [9009,8448,8118,8008,7227,7007,6776,6336,6006,5775,5445,5335,
--     5225,5115,5005,4884,4774,4664,4554,4224,4004,3773,3663,3003,
--     2992,2772,2552,2442,2332,2112,2002,1881,1771,1551,1221,1001]
capicuasP n = [x | x <- capicuas n,
                        not (null (productosDosNumerosCifras n x))]

-- (capicuas n) es la lista de las capicúas de 2*n cifras de mayor a
-- menor. Por ejemplo, 
--    capicuas 1           ==  [99,88,77,66,55,44,33,22,11]
--    take 7 (capicuas 2)  ==  [9999,9889,9779,9669,9559,9449,9339]
capicuas :: Integer -> [Integer]
capicuas n = [capicua x | x <- numerosCifras n]

-- (numerosCifras n) es la lista de los números de n cifras de mayor a
-- menor. Por ejemplo,
--    numerosCifras 1           ==  [9,8,7,6,5,4,3,2,1]
--    take 7 (numerosCifras 2)  ==  [99,98,97,96,95,94,93]
--    take 7 (numerosCifras 3)  ==  [999,998,997,996,995,994,993]
numerosCifras :: Integer -> [Integer]
numerosCifras n = [a,a-1..b]
  where a = 10^n-1
        b = 10^(n-1) 

-- (capicua n) es la capicúa formada añadiendo el inverso de n a
--  continuación de n. Por ejemplo,
--    capicua 93  ==  9339
capicua :: Integer -> Integer
capicua n = read (xs ++ (reverse xs))
  where xs = show n

-- (productosDosNumerosCifras n x) es la lista de los números y de n
-- cifras tales que existe un z de n cifras y x es el producto de y por
-- z. Por ejemplo, 
--    productosDosNumerosCifras 2 9009  ==  [99,91]
productosDosNumerosCifras n x = [y | y <- numeros,
                                     mod x y == 0,
                                     div x y `elem` numeros]
  where numeros = numerosCifras n

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

mayorCapicuaP2 :: Integer -> Integer
mayorCapicuaP2 n = maximum [x*y | x <- [a,a-1..b],
                                  y <- [a,a-1..b],
                                  esCapicua (x*y)] 
  where a = 10^n-1
        b = 10^(n-1)

-- (esCapicua x) se verifica si x es capicúa. Por ejemplo,
--    esCapicua 353  ==  True
--    esCapicua 357  ==  False
esCapicua :: Integer -> Bool
esCapicua n = xs == reverse xs
  where xs = show n

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

mayorCapicuaP3 :: Integer -> Integer
mayorCapicuaP3 n = maximum [x*y | (x,y) <- pares a b, 
                                  esCapicua (x*y)] 
  where a = 10^n-1
        b = 10^(n-1)

-- (pares a b) es la lista de los pares de números entre a y b de forma
-- que su suma es decreciente. Por ejemplo,
--    pares 9 7  ==  [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]
pares a b = [(x,z-x) | z <- [a1,a1-1..b1],
                       x <- [a,a-1..b],
                       x <= z-x, z-x <= a]
  where a1 = 2*a
        b1 = 2*b

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

mayorCapicuaP4 :: Integer -> Integer
mayorCapicuaP4 n = maximum [x | y <- [a..b],
                                z <- [y..b],
                                let x = y * z,
                                let s = show x,
                                s == reverse s]
  where a = 10^(n-1)
        b = 10^n-1
     
-- 5ª solución
-- ===========

mayorCapicuaP5 :: Integer -> Integer
mayorCapicuaP5 n = maximum [x*y | (x,y) <- pares2 b a, esCapicua (x*y)]
  where a = 10^(n-1)
        b = 10^n-1
                       
-- (pares2 a b) es la lista de los pares de números entre a y b de forma
-- que su suma es decreciente. Por ejemplo,
--    pares2 9 7  ==  [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]
pares2 a b = [(x,y) | x <- [a,a-1..b], y <- [a,a-1..x]]

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

mayorCapicuaP6 :: Integer -> Integer
mayorCapicuaP6 n = maximum [x*y | x <- [a..b], 
                                  y <- [x..b] , 
                                  esCapicua (x*y)]
  where a = 10^(n-1)
        b = 10^n-1
 
-- (cifras n) es la lista de las cifras de n en orden inverso. Por
-- ejemplo,  
--    cifras 325  == [5,2,3]
cifras :: Integer -> [Integer]
cifras n 
    | n < 10    = [n]
    | otherwise = (ultima n) : (cifras (quitarUltima n))

-- (ultima n) es la última cifra de n. Por ejemplo,
--    ultima 325  ==  5
ultima  :: Integer -> Integer
ultima n =  n - (n `div` 10)*10

-- (quitarUltima n) es el número obtenido al quitarle a n su última
-- cifra. Por ejemplo,
--    quitarUltima 325  =>  32 
quitarUltima :: Integer -> Integer
quitarUltima n = (n - (ultima n)) `div` 10

-- 7ª solución
-- ===========

mayorCapicuaP7 :: Integer -> Integer
mayorCapicuaP7 n = head [x | x <- capicuas n, esFactorizable x n]

-- (esFactorizable x n) se verifica si x se puede escribir como producto
-- de dos números de n dígitos. Por ejemplo,
--    esFactorizable 1219 2  ==  True
--    esFactorizable 1217 2  ==  False
esFactorizable x n = aux i x
  where b = 10^n-1
        i = floor (sqrt (fromIntegral x))
        aux i x | i > b          = False
                | x `mod` i == 0 = x `div` i < b 
                | otherwise      = aux (i+1) x

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

--    λ> mayorCapicuaP1 3
--    906609
--    (0.07 secs, 18,248,224 bytes)
--    λ> mayorCapicuaP2 3
--    906609
--    (0.51 secs, 555,695,720 bytes)
--    λ> mayorCapicuaP3 3
--    906609
--    (0.96 secs, 780,794,768 bytes)
--    λ> mayorCapicuaP4 3
--    906609
--    (0.24 secs, 255,445,448 bytes)
--    λ> mayorCapicuaP5 3
--    906609
--    (0.33 secs, 317,304,080 bytes)
--    λ> mayorCapicuaP6 3
--    906609
--    (0.26 secs, 274,987,472 bytes)
--    λ> mayorCapicuaP7 3
--    906609
--    (0.02 secs, 1,807,720 bytes)
--    
--    λ> mayorCapicuaP1 5
--    9966006699
--    (9.90 secs, 6,349,454,544 bytes)
--    λ> mayorCapicuaP7 5
--    9966006699
--    (0.06 secs, 15,958,616 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

176

177

178

179

180

181

182

183

184

185

186

-- 1ª solución

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

mayorCapicuaP1 :: Integer -> Integer

mayorCapicuaP1 n = head (capicuasP n)

-- (capicuasP n) es la lista de las capicúas de 2*n cifras que

-- pueden escribirse como productos de dos números de n cifras. Por

-- ejemplo, Por ejemplo,

-- ghci> capicuasP 2

-- [9009,8448,8118,8008,7227,7007,6776,6336,6006,5775,5445,5335,

-- 5225,5115,5005,4884,4774,4664,4554,4224,4004,3773,3663,3003,

-- 2992,2772,2552,2442,2332,2112,2002,1881,1771,1551,1221,1001]

capicuasP n = [x | x <- capicuas n,

not (null (productosDosNumerosCifras n x))]

-- (capicuas n) es la lista de las capicúas de 2*n cifras de mayor a

-- menor. Por ejemplo,

-- capicuas 1 == [99,88,77,66,55,44,33,22,11]

-- take 7 (capicuas 2) == [9999,9889,9779,9669,9559,9449,9339]

capicuas :: Integer -> [Integer]

capicuas n = [capicua x | x <- numerosCifras n]

-- (numerosCifras n) es la lista de los números de n cifras de mayor a

-- menor. Por ejemplo,

-- numerosCifras 1 == [9,8,7,6,5,4,3,2,1]

-- take 7 (numerosCifras 2) == [99,98,97,96,95,94,93]

-- take 7 (numerosCifras 3) == [999,998,997,996,995,994,993]

numerosCifras :: Integer -> [Integer]

numerosCifras n = [a,a-1..b]

where a = 10^n-1

b = 10^(n-1)

-- (capicua n) es la capicúa formada añadiendo el inverso de n a

-- continuación de n. Por ejemplo,

-- capicua 93 == 9339

capicua :: Integer -> Integer

capicua n = read (xs ++ (reverse xs))

where xs = show n

-- (productosDosNumerosCifras n x) es la lista de los números y de n

-- cifras tales que existe un z de n cifras y x es el producto de y por

-- z. Por ejemplo,

-- productosDosNumerosCifras 2 9009 == [99,91]

productosDosNumerosCifras n x = [y | y <- numeros,

mod x y == 0,

div x y `elem` numeros]

where numeros = numerosCifras n

-- 2ª solución

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

mayorCapicuaP2 :: Integer -> Integer

mayorCapicuaP2 n = maximum [x*y | x <- [a,a-1..b],

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

esCapicua (x*y)]

where a = 10^n-1

b = 10^(n-1)

-- (esCapicua x) se verifica si x es capicúa. Por ejemplo,

-- esCapicua 353 == True

-- esCapicua 357 == False

esCapicua :: Integer -> Bool

esCapicua n = xs == reverse xs

where xs = show n

-- 3ª solución

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

mayorCapicuaP3 :: Integer -> Integer

mayorCapicuaP3 n = maximum [x*y | (x,y) <- pares a b,

esCapicua (x*y)]

where a = 10^n-1

b = 10^(n-1)

-- (pares a b) es la lista de los pares de números entre a y b de forma

-- que su suma es decreciente. Por ejemplo,

-- pares 9 7 == [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]

pares a b = [(x,z-x) | z <- [a1,a1-1..b1],

x <- [a,a-1..b],

x <= z-x, z-x <= a]

where a1 = 2*a

b1 = 2*b

-- 4ª solución

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

mayorCapicuaP4 :: Integer -> Integer

mayorCapicuaP4 n = maximum [x | y <- [a..b],

z <- [y..b],

let x = y * z,

let s = show x,

s == reverse s]

where a = 10^(n-1)

b = 10^n-1

-- 5ª solución

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

mayorCapicuaP5 :: Integer -> Integer

mayorCapicuaP5 n = maximum [x*y | (x,y) <- pares2 b a, esCapicua (x*y)]

where a = 10^(n-1)

b = 10^n-1

-- (pares2 a b) es la lista de los pares de números entre a y b de forma

-- que su suma es decreciente. Por ejemplo,

-- pares2 9 7 == [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]

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

-- 6ª solución

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

mayorCapicuaP6 :: Integer -> Integer

mayorCapicuaP6 n = maximum [x*y | x <- [a..b],

y <- [x..b] ,

esCapicua (x*y)]

where a = 10^(n-1)

b = 10^n-1

-- (cifras n) es la lista de las cifras de n en orden inverso. Por

-- ejemplo,

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

cifras :: Integer -> [Integer]

cifras n

| n < 10 = [n]

| otherwise = (ultima n) : (cifras (quitarUltima n))

-- (ultima n) es la última cifra de n. Por ejemplo,

-- ultima 325 == 5

ultima :: Integer -> Integer

ultima n = n - (n `div` 10)*10

-- (quitarUltima n) es el número obtenido al quitarle a n su última

-- cifra. Por ejemplo,

-- quitarUltima 325 => 32

quitarUltima :: Integer -> Integer

quitarUltima n = (n - (ultima n)) `div` 10

-- 7ª solución

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

mayorCapicuaP7 :: Integer -> Integer

mayorCapicuaP7 n = head [x | x <- capicuas n, esFactorizable x n]

-- (esFactorizable x n) se verifica si x se puede escribir como producto

-- de dos números de n dígitos. Por ejemplo,

-- esFactorizable 1219 2 == True

-- esFactorizable 1217 2 == False

esFactorizable x n = aux i x

where b = 10^n-1

i = floor (sqrt (fromIntegral x))

aux i x | i > b = False

| x `mod` i == 0 = x `div` i < b

| otherwise = aux (i+1) x

-- Comparación de eficiencia

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

-- λ> mayorCapicuaP1 3

-- 906609

-- (0.07 secs, 18,248,224 bytes)

-- λ> mayorCapicuaP2 3

-- 906609

-- (0.51 secs, 555,695,720 bytes)

-- λ> mayorCapicuaP3 3

-- 906609

-- (0.96 secs, 780,794,768 bytes)

-- λ> mayorCapicuaP4 3

-- 906609

-- (0.24 secs, 255,445,448 bytes)

-- λ> mayorCapicuaP5 3

-- 906609

-- (0.33 secs, 317,304,080 bytes)

-- λ> mayorCapicuaP6 3

-- 906609

-- (0.26 secs, 274,987,472 bytes)

-- λ> mayorCapicuaP7 3

-- 906609

-- (0.02 secs, 1,807,720 bytes)

-- λ> mayorCapicuaP1 5

-- 9966006699

-- (9.90 secs, 6,349,454,544 bytes)

-- λ> mayorCapicuaP7 5

-- 9966006699

-- (0.06 secs, 15,958,616 bytes)

Combinaciones divisibles

PorJosé A. Alonso 21-abril-202029-abril-2020

Definir la función

   tieneCombinacionDivisible :: [Int] -> Int -> Bool

1	tieneCombinacionDivisible :: [Int] -> Int -> Bool

tal que (tieneCombinacionDivisible xs m) se verifica si existe alguna forma de combinar todos los elementos de la lista (con las operaciones suma o resta) de forma que el resultado sea divisible por m. Por ejemplo,

   tieneCombinacionDivisible [1,3,4,6] 4  ==  True
   tieneCombinacionDivisible [1,3,9]   2  ==  False

1 2	tieneCombinacionDivisible [1,3,4,6] 4 == True tieneCombinacionDivisible [1,3,9] 2 == False

En el primer ejemplo, 1 – 2 + 3 + 4 + 6 = 12 es una combinación divisible por 4. En el segundo ejemplo, las combinaciones de [1,3,9] son

   1 + 3 + 9 =  13
  -1 + 3 + 9 =  11
   1 - 3 + 9 =   7
  -1 - 3 + 9 =   5
   1 + 3 - 9 =  -5
  -1 + 3 - 9 =  -7
   1 - 3 - 9 = -11
  -1 - 3 - 9 = -13

1 + 3 + 9 = 13

-1 + 3 + 9 = 11

1 - 3 + 9 = 7

-1 - 3 + 9 = 5

1 + 3 - 9 = -5

-1 + 3 - 9 = -7

1 - 3 - 9 = -11

-1 - 3 - 9 = -13

y ninguna de las 4 es divisible por 2.

Soluciones

import Test.QuickCheck

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

tieneCombinacionDivisible :: [Int] -> Int -> Bool
tieneCombinacionDivisible xs m =
  any esDivisible (valoresCombinaciones xs)
  where esDivisible x = x `mod` m == 0

-- (valoresCombinaciones xs) es la lista de los valores de todas las
-- combinaciones de todos los elementos de la lista con las operaciones
-- suma o resta. Por ejemplo,
--    λ> valoresCombinaciones [1,3,4,6]
--    [14,12,8,6,6,4,0,-2,2,0,-4,-6,-6,-8,-12,-14]
--    λ> valoresCombinaciones [1,3,-4,6]
--    [6,4,0,-2,14,12,8,6,-6,-8,-12,-14,2,0,-4,-6]
valoresCombinaciones :: [Int] -> [Int]
valoresCombinaciones []     = []
valoresCombinaciones [x]    = [x,-x]
valoresCombinaciones (x:xs) = concat [[y + x, y - x] | y <- ys]
  where ys = valoresCombinaciones xs

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

tieneCombinacionDivisible2 :: [Int] -> Int -> Bool
tieneCombinacionDivisible2 xs m =
  tieneCombinacionCongruente xs m 0

-- (tieneCombinacionCongruente xs m a) se verifica si existe alguna
-- forma de combinar todos los elementos de la lista xs (con las
-- operaciones suma o resta) de forma que el resultado sea congruente
-- con a módulo m. Por ejemplo,
--    tieneCombinacionCongruente [1,3,4,6] 4 0  ==  True
--    tieneCombinacionCongruente [1,3,4,6] 4 1  ==  False
--    tieneCombinacionCongruente [1,3,9] 2 0    ==  False
--    tieneCombinacionCongruente [1,3,9] 2 1    ==  True
tieneCombinacionCongruente :: [Int] -> Int -> Int -> Bool
tieneCombinacionCongruente []  _  _ = False
tieneCombinacionCongruente [x] m  a = (x - a) `mod` m == 0
tieneCombinacionCongruente (x:xs) m a =
  tieneCombinacionCongruente xs m (a-x) ||
  tieneCombinacionCongruente xs m (a+x)

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

-- La propiedad es
prop_tieneCombinacionDivisible :: [Int] -> Positive Int -> Bool
prop_tieneCombinacionDivisible xs (Positive m) =
  tieneCombinacionDivisible xs m == tieneCombinacionDivisible2 xs m

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

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

--    λ> (n,xs,m) = (200,[-n..n],sum [1..n]) 
--    (0.00 secs, 0 bytes)
--    λ> and [tieneCombinacionDivisible xs a | a <- [1..m]]
--    True
--    (4.74 secs, 6,536,494,976 bytes)
--    λ> and [tieneCombinacionDivisible2 xs a | a <- [1..m]]
--    True
--    (2.97 secs, 3,381,932,664 bytes)

import Test.QuickCheck

-- 1ª solución

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

tieneCombinacionDivisible :: [Int] -> Int -> Bool

tieneCombinacionDivisible xs m =

any esDivisible (valoresCombinaciones xs)

where esDivisible x = x `mod` m == 0

-- (valoresCombinaciones xs) es la lista de los valores de todas las

-- combinaciones de todos los elementos de la lista con las operaciones

-- suma o resta. Por ejemplo,

-- λ> valoresCombinaciones [1,3,4,6]

-- [14,12,8,6,6,4,0,-2,2,0,-4,-6,-6,-8,-12,-14]

-- λ> valoresCombinaciones [1,3,-4,6]

-- [6,4,0,-2,14,12,8,6,-6,-8,-12,-14,2,0,-4,-6]

valoresCombinaciones :: [Int] -> [Int]

valoresCombinaciones [] = []

valoresCombinaciones [x] = [x,-x]

valoresCombinaciones (x:xs) = concat [[y + x, y - x] | y <- ys]

where ys = valoresCombinaciones xs

-- 2ª solución

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

tieneCombinacionDivisible2 :: [Int] -> Int -> Bool

tieneCombinacionDivisible2 xs m =

tieneCombinacionCongruente xs m 0

-- (tieneCombinacionCongruente xs m a) se verifica si existe alguna

-- forma de combinar todos los elementos de la lista xs (con las

-- operaciones suma o resta) de forma que el resultado sea congruente

-- con a módulo m. Por ejemplo,

-- tieneCombinacionCongruente [1,3,4,6] 4 0 == True

-- tieneCombinacionCongruente [1,3,4,6] 4 1 == False

-- tieneCombinacionCongruente [1,3,9] 2 0 == False

-- tieneCombinacionCongruente [1,3,9] 2 1 == True

tieneCombinacionCongruente :: [Int] -> Int -> Int -> Bool

tieneCombinacionCongruente [] _ _ = False

tieneCombinacionCongruente [x] m a = (x - a) `mod` m == 0

tieneCombinacionCongruente (x:xs) m a =

tieneCombinacionCongruente xs m (a-x) ||

tieneCombinacionCongruente xs m (a+x)

-- Equivalencia

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

-- La propiedad es

prop_tieneCombinacionDivisible :: [Int] -> Positive Int -> Bool

prop_tieneCombinacionDivisible xs (Positive m) =

tieneCombinacionDivisible xs m == tieneCombinacionDivisible2 xs m

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=25}) prop_tieneCombinacionDivisible

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> (n,xs,m) = (200,[-n..n],sum [1..n])

-- (0.00 secs, 0 bytes)

-- λ> and [tieneCombinacionDivisible xs a | a <- [1..m]]

-- True

-- (4.74 secs, 6,536,494,976 bytes)

-- λ> and [tieneCombinacionDivisible2 xs a | a <- [1..m]]

-- True

-- (2.97 secs, 3,381,932,664 bytes)

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>

Orden de divisibilidad

PorJosé A. Alonso 14-abril-202021-abril-2020

El orden de divisibilidad de un número x es el mayor n tal que para todo i menor o igual que n, los i primeros dígitos de n es divisible por i. Por ejemplo, el orden de divisibilidad de 74156 es 3 porque

   7       es divisible por 1
   74      es divisible por 2
   741     es divisible por 3
   7415 no es divisible por 4

7 es divisible por 1

74 es divisible por 2

741 es divisible por 3

7415 no es divisible por 4

Definir la función

   ordenDeDivisibilidad :: Integer -> Int

1	ordenDeDivisibilidad :: Integer -> Int

tal que (ordenDeDivisibilidad x) es el orden de divisibilidad de x. Por ejemplo,

   ordenDeDivisibilidad 74156                      ==  3
   ordenDeDivisibilidad 12                         ==  2
   ordenDeDivisibilidad 7                          ==  1
   ordenDeDivisibilidad 3608528850368400786036725  ==  25

ordenDeDivisibilidad 74156 == 3

ordenDeDivisibilidad 12 == 2

ordenDeDivisibilidad 7 == 1

ordenDeDivisibilidad 3608528850368400786036725 == 25

Soluciones

import Data.List (inits)

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

ordenDeDivisibilidad :: Integer -> Int
ordenDeDivisibilidad n = 
  length (takeWhile (\(x,k) -> x `mod` k == 0) (zip (sucDigitos n) [1..]))

-- (sucDigitos x) es la sucesión de los dígitos de x. Por ejemplo,
--    sucDigitos 325    ==  [3,32,325]
--    sucDigitos 32050  ==  [3,32,320,3205,32050]
sucDigitos :: Integer -> [Integer]
sucDigitos n = 
    [n `div` (10^i) | i <- [k-1,k-2..0]]
    where k = length (show n)

-- 2ª definición de sucDigitos
sucDigitos2 :: Integer -> [Integer]
sucDigitos2 n = [read xs | xs <- aux (show n)]
  where aux []     = []
        aux (d:ds) = [d] : map (d:) (aux ds)

-- 3ª definición de sucDigitos
sucDigitos3 :: Integer -> [Integer]
sucDigitos3 n = 
  [read (take k ds) | k <- [1..length ds]]
  where ds = show n

-- 4ª definición de sucDigitos
sucDigitos4 :: Integer -> [Integer]
sucDigitos4 n = [read xs | xs <- tail (inits (show n))]

-- 5ª definición de sucDigitos
sucDigitos5 :: Integer -> [Integer]
sucDigitos5 n = map read (tail (inits (show n)))

-- 6ª definición de sucDigitos
sucDigitos6 :: Integer -> [Integer]
sucDigitos6 = map read . (tail . inits . show)

-- Eficiencia de las definiciones de sucDigitos
--    ghci> length (sucDigitos (10^5000))
--    5001
--    (0.01 secs, 1550688 bytes)
--    ghci> length (sucDigitos2 (10^5000))
--    5001
--    (1.25 secs, 729411872 bytes)
--    ghci> length (sucDigitos3 (10^5000))
--    5001
--    (0.02 secs, 2265120 bytes)
--    ghci> length (sucDigitos4 (10^5000))
--    5001
--    (1.10 secs, 728366872 bytes)
--    ghci> length (sucDigitos5 (10^5000))
--    5001
--    (1.12 secs, 728393864 bytes)
--    ghci> length (sucDigitos6 (10^5000))
--    5001
--    (1.20 secs, 728403052 bytes)
-- 
--    ghci> length (sucDigitos (10^3000000))
--    3000001
--    (2.73 secs, 820042696 bytes)
--    ghci> length (sucDigitos3 (10^3000000))
--    3000001
--    (3.69 secs, 820043688 bytes)

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

ordenDeDivisibilidad2 :: Integer -> Int
ordenDeDivisibilidad2 x =
  length
  $ takeWhile (==0)
  $ zipWith (mod . read) (tail $ inits $ show x) [1..]

import Data.List (inits)

-- 1ª definición de ordenDeDivisibilidad

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

ordenDeDivisibilidad :: Integer -> Int

ordenDeDivisibilidad n =

length (takeWhile (\(x,k) -> x `mod` k == 0) (zip (sucDigitos n) [1..]))

-- (sucDigitos x) es la sucesión de los dígitos de x. Por ejemplo,

-- sucDigitos 325 == [3,32,325]

-- sucDigitos 32050 == [3,32,320,3205,32050]

sucDigitos :: Integer -> [Integer]

sucDigitos n =

[n `div` (10^i) | i <- [k-1,k-2..0]]

where k = length (show n)

-- 2ª definición de sucDigitos

sucDigitos2 :: Integer -> [Integer]

sucDigitos2 n = [read xs | xs <- aux (show n)]

where aux [] = []

aux (d:ds) = [d] : map (d:) (aux ds)

-- 3ª definición de sucDigitos

sucDigitos3 :: Integer -> [Integer]

sucDigitos3 n =

[read (take k ds) | k <- [1..length ds]]

where ds = show n

-- 4ª definición de sucDigitos

sucDigitos4 :: Integer -> [Integer]

sucDigitos4 n = [read xs | xs <- tail (inits (show n))]

-- 5ª definición de sucDigitos

sucDigitos5 :: Integer -> [Integer]

sucDigitos5 n = map read (tail (inits (show n)))

-- 6ª definición de sucDigitos

sucDigitos6 :: Integer -> [Integer]

sucDigitos6 = map read . (tail . inits . show)

-- Eficiencia de las definiciones de sucDigitos

-- ghci> length (sucDigitos (10^5000))

-- 5001

-- (0.01 secs, 1550688 bytes)

-- ghci> length (sucDigitos2 (10^5000))

-- 5001

-- (1.25 secs, 729411872 bytes)

-- ghci> length (sucDigitos3 (10^5000))

-- 5001

-- (0.02 secs, 2265120 bytes)

-- ghci> length (sucDigitos4 (10^5000))

-- 5001

-- (1.10 secs, 728366872 bytes)

-- ghci> length (sucDigitos5 (10^5000))

-- 5001

-- (1.12 secs, 728393864 bytes)

-- ghci> length (sucDigitos6 (10^5000))

-- 5001

-- (1.20 secs, 728403052 bytes)

-- ghci> length (sucDigitos (10^3000000))

-- 3000001

-- (2.73 secs, 820042696 bytes)

-- ghci> length (sucDigitos3 (10^3000000))

-- 3000001

-- (3.69 secs, 820043688 bytes)

-- 2ª definición de ordenDeDivisibilidad

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

ordenDeDivisibilidad2 :: Integer -> Int

ordenDeDivisibilidad2 x =

length

$ takeWhile (==0)

$ zipWith (mod . read) (tail $ inits $ show x) [1..]

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>

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>

El sesgo de Chebyshev

PorJosé A. Alonso 23-marzo-202026-marzo-2022

Un número primo distinto de 2 tiene la forma 4k + 1 o 4k + 3. Chebyshev notó en 1853 que la mayoría de las veces hay más números primos de la forma 4k + 3 que números primos de la forma 4k + 1 menores que un número dado. Esto se llama el sesgo de Chebyshev.

Definir las funciones

   distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]
   empatesRestosModulo4 :: [Integer]
   mayoria1RestosModulo4 :: [Integer]
   grafica_Chebyshev :: Int -> IO ()

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]

empatesRestosModulo4 :: [Integer]

mayoria1RestosModulo4 :: [Integer]

grafica_Chebyshev :: Int -> IO ()

tales que

distribucionPrimosModulo4 es la lista de las ternas (p,a,b) tales que p es un números primo, a es la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 y b es la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 7 distribucionPrimosModulo4
     [(2,0,0),(3,0,1),(5,1,1),(7,1,2),(11,1,3),(13,2,3),(17,3,3)]
     λ> distribucionPrimosModulo4 !! (5*10^5)
     (7368791,249888,250112)

λ> take 7 distribucionPrimosModulo4

[(2,0,0),(3,0,1),(5,1,1),(7,1,2),(11,1,3),(13,2,3),(17,3,3)]

λ> distribucionPrimosModulo4 !! (5*10^5)

(7368791,249888,250112)

empatesRestosModulo4 es la lista de los primos p tales que la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 es igual a la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 10 empatesRestosModulo4
     [2,5,17,41,461,26833,26849,26863,26881,26893]
     λ> length (takeWhile (<= 10^6) empatesRestosModulo4)
     112

λ> take 10 empatesRestosModulo4

[2,5,17,41,461,26833,26849,26863,26881,26893]

λ> length (takeWhile (<= 10^6) empatesRestosModulo4)

112

mayoria1RestosModulo4 es la lista de los primos p tales que la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 es mayor que la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 10 mayoria1RestosModulo4
     [26861,616841,616849,616877,616897,616909,616933,616943,616951,616961]
     λ> length (takeWhile (<= 10^6) mayoria1RestosModulo4)
     239

λ> take 10 mayoria1RestosModulo4

[26861,616841,616849,616877,616897,616909,616933,616943,616951,616961]

λ> length (takeWhile (<= 10^6) mayoria1RestosModulo4)

239

(graficaChebyshev n) dibuja la gráfica de los puntos (p,b-a) donde p es uno de los n primeros primos impares, a es la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 y b es la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo, (graficaChebyshev 5000) dibuja la figura

Soluciones

[schedule expon=’2020-03-30′ expat=»06:00″]

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

Pensamiento

«El valor de un problema no es tanto el de encontrar la respuesta como el de las ideas e intentos que obliga su resolución.»

Israel Nathan Herstein.

[/schedule]

[schedule on=’2020-03-30′ at=»06:00″]

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

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]
distribucionPrimosModulo4 = (2,0,0) : aux (tail primes) (0, 0)
  where aux (p:ps) (a,b)
          | p `mod` 4 == 1 = (p,a+1,b) : aux ps (a+1,b)
          | otherwise      = (p,a,b+1) : aux ps (a,b+1)

empatesRestosModulo4 :: [Integer]
empatesRestosModulo4 =
  [p | (p,a,b) <- distribucionPrimosModulo4
     , a == b]
  
mayoria1RestosModulo4 :: [Integer]
mayoria1RestosModulo4 =
  [p | (p,a,b) <- distribucionPrimosModulo4
     , a > b]
  
grafica :: Int -> IO ()
grafica n = 
  plotLists [Key Nothing]
            [ [(p, a) | (p,a,b) <- xs]
            , [(p, b) | (p,a,b) <- xs]
            , [(p, b-a) | (p,a,b) <- xs]]
  where xs = take n (tail (distribucionPrimosModulo4))

graficaChebyshev :: Int -> IO ()
graficaChebyshev n = 
  plotList [ Key Nothing
           , PNG "El_sesgo_de_Chebyshev.png"
           ]
           [(p, b-a) | (p,a,b) <- xs]
  where xs = take n (tail (distribucionPrimosModulo4))

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]

distribucionPrimosModulo4 = (2,0,0) : aux (tail primes) (0, 0)

where aux (p:ps) (a,b)

| p `mod` 4 == 1 = (p,a+1,b) : aux ps (a+1,b)

| otherwise = (p,a,b+1) : aux ps (a,b+1)

empatesRestosModulo4 :: [Integer]

empatesRestosModulo4 =

[p | (p,a,b) <- distribucionPrimosModulo4

, a == b]

mayoria1RestosModulo4 :: [Integer]

mayoria1RestosModulo4 =

[p | (p,a,b) <- distribucionPrimosModulo4

, a > b]

grafica :: Int -> IO ()

grafica n =

plotLists [Key Nothing]

[ [(p, a) | (p,a,b) <- xs]

, [(p, b) | (p,a,b) <- xs]

, [(p, b-a) | (p,a,b) <- xs]]

where xs = take n (tail (distribucionPrimosModulo4))

graficaChebyshev :: Int -> IO ()

graficaChebyshev n =

plotList [ Key Nothing

, PNG "El_sesgo_de_Chebyshev.png"

]

[(p, b-a) | (p,a,b) <- xs]

where xs = take n (tail (distribucionPrimosModulo4))

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>

[/schedule]

La menos conocida de las conjeturas de Goldbach

PorJosé A. Alonso 7-febrero-202017-febrero-2020

Goldbach, el de la famosa conjetura, hizo por lo menos otra conjetura que finalmente resultó ser falsa.

Esta última decía que todo número compuesto impar puede expresarse como la suma de un número primo más dos veces la suma de un cuadrado. Así por ejemplo,

    9 =  7 + 2×1^2
   15 =  7 + 2×2^2
   21 =  3 + 2×3^2
   25 =  7 + 2×3^2
   27 = 19 + 2×2^2
   33 = 31 + 2×1^2

9 = 7 + 2×1^2

15 = 7 + 2×2^2

21 = 3 + 2×3^2

25 = 7 + 2×3^2

27 = 19 + 2×2^2

33 = 31 + 2×1^2

Definir las sucesiones

   imparesCompuestos :: [Integer]
   descomposiciones :: Integer -> [(Integer,Integer)]
   contraejemplosGoldbach :: [Integer]

imparesCompuestos :: [Integer]

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

contraejemplosGoldbach :: [Integer]

tales que

imparesCompuestos es la lista de los números impares compuestos. Por ejemplo,

     take 9 imparesCompuestos  ==  [9,15,21,25,27,33,35,39,45]

1	take 9 imparesCompuestos == [9,15,21,25,27,33,35,39,45]

(descomposiciones n) es la lista de las descomposiciones de n de n como la suma de un número primo más dos veces la suma de un cuadrado. Por ejemplo,

     descomposiciones 9     ==  [(7,1)]
     descomposiciones 21    ==  [(3,9),(13,4),(19,1)]
     descomposiciones 5777  ==  []

descomposiciones 9 == [(7,1)]

descomposiciones 21 == [(3,9),(13,4),(19,1)]

descomposiciones 5777 == []

Las 3 descomposiciones de 21 son

     21 =  3 + 2*9 = 21 + 2*3^2
     21 = 13 + 2*4 = 13 + 2*3^2
     21 = 19 + 2*1 = 19 + 2*1^2

21 = 3 + 2*9 = 21 + 2*3^2

21 = 13 + 2*4 = 13 + 2*3^2

21 = 19 + 2*1 = 19 + 2*1^2

contraejemplosGoldbach es la lista de los contraejemplos de la anterior conjetura de Goldbach; es decir, los números impares compuestos que no pueden expresarse como la suma de un número primo más dos veces la suma de un cuadrado. Por ejemplo,

   take 2 contraejemplosGoldbach  ==  [5777,5993]

1	take 2 contraejemplosGoldbach == [5777,5993]

Comprobar con QuickCheck que la conjetura de Golbach se verifica a partir de 5993; es decir, todo número compuesto impar mayor que 5993 puede expresarse como la suma de un número primo más dos veces la suma de un cuadrado.

Nota: Basado en el artículo La menos conocida de las conjeturas de Goldbach de Claudio Meller en el blog Números y algo más.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

imparesCompuestos :: [Integer]
imparesCompuestos = filter esCompuesto [3,5..]

-- (esCompuesto x) se verifica si x es un número compuesto. Por ejemplo,
--    esCompuesto 6  ==  True
--    esCompuesto 7  ==  False
esCompuesto :: Integer -> Bool
esCompuesto = not . isPrime

contraejemplosGoldbach :: [Integer]
contraejemplosGoldbach = filter esContraejemplo imparesCompuestos

-- (esContraejemplo x) es verifica si el número impar compuesto x es un
-- contraejemplo de la conjetura de Goldbach. Por ejemplo,
--    esContraejemplo 5777  ==  True
--    esContraejemplo 15    ==  False
esContraejemplo :: Integer -> Bool
esContraejemplo = null . descomposiciones

descomposiciones :: Integer -> [(Integer,Integer)]
descomposiciones n =
  [(p,x) | p <- takeWhile (<=n) primes
         , (n - p) `mod` 2 == 0
         , let x = (n - p) `div` 2
         , esCuadrado x]

-- (esCuadrado x) es verifica si x es un cuadrado perfecto. Por ejemplo, 
--    esCuadrado 16  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = y^2 == x
  where y = ceiling (sqrt (fromIntegral x))

-- La propiedad es
prop_conjetura :: Int -> Property
prop_conjetura n =
  n >= 0 ==> not (esContraejemplo (imparesCompuestosMayore5993 !! n))
  where imparesCompuestosMayore5993 = dropWhile (<=5993) imparesCompuestos

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

import Data.Numbers.Primes

import Test.QuickCheck

imparesCompuestos :: [Integer]

imparesCompuestos = filter esCompuesto [3,5..]

-- (esCompuesto x) se verifica si x es un número compuesto. Por ejemplo,

-- esCompuesto 6 == True

-- esCompuesto 7 == False

esCompuesto :: Integer -> Bool

esCompuesto = not . isPrime

contraejemplosGoldbach :: [Integer]

contraejemplosGoldbach = filter esContraejemplo imparesCompuestos

-- (esContraejemplo x) es verifica si el número impar compuesto x es un

-- contraejemplo de la conjetura de Goldbach. Por ejemplo,

-- esContraejemplo 5777 == True

-- esContraejemplo 15 == False

esContraejemplo :: Integer -> Bool

esContraejemplo = null . descomposiciones

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

descomposiciones n =

[(p,x) | p <- takeWhile (<=n) primes

, (n - p) `mod` 2 == 0

, let x = (n - p) `div` 2

, esCuadrado x]

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

-- esCuadrado 16 == True

-- esCuadrado 27 == False

esCuadrado :: Integer -> Bool

esCuadrado x = y^2 == x

where y = ceiling (sqrt (fromIntegral x))

-- La propiedad es

prop_conjetura :: Int -> Property

prop_conjetura n =

n >= 0 ==> not (esContraejemplo (imparesCompuestosMayore5993 !! n))

where imparesCompuestosMayore5993 = dropWhile (<=5993) imparesCompuestos

-- La comprobación es

-- λ> quickCheck prop_conjetura

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

«Obvio es la palabra más peligrosa de las matemáticas.»

Eric Temple Bell

Longitud de la parte periódica

PorJosé A. Alonso 30-enero-202026-marzo-2022

La propiedad de la longitud de la parte periódica afirma que

Si p es un número primo distinto de 2 y de 5, entonces la longitud del período de 1/p es el menor entero positivo n tal que p divide a $10^n - 1$ .

El objetivo de este ejercicio es la verificación de dicha propiedad.

Las fracciones se representan por un par de enteros. Por ejemplo, el número 2/3 se representa por (2,3). Su tipo es

   type Fraccion = (Integer,Integer)

1	type Fraccion = (Integer,Integer)

Los números decimales se representan por ternas, donde el primer elemento es la parte entera, el segundo es el anteperíodo y el tercero es el período. Por ejemplo,

 6/2  = 3                  se representa por (3,[],[])
 1/2  = 0.5                se representa por (0,[5],[])
 1/3  = 0.333333...        se representa por (0,[],[3])  
23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

6/2 = 3 se representa por (3,[],[])

1/2 = 0.5 se representa por (0,[5],[])

1/3 = 0.333333... se representa por (0,[],[3])

23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

Su tipo es

   type Decimal = (Integer,[Integer],[Integer])

1	type Decimal = (Integer,[Integer],[Integer])

Definir, usando las funciones cocientesRestos y primerRepetido de los ejercicios anteriores, las funciones

   decimal :: Fraccion -> Decimal
   longitudPeriodo :: Fraccion -> Integer

1 2	decimal :: Fraccion -> Decimal longitudPeriodo :: Fraccion -> Integer

tales que

(decimal f) es la representación decimal de la fracción f. Por ejemplo,

     decimal (6,2)          ==  (3,[],[])
     decimal (3,4)          ==  (0,[7,5],[])
     decimal (1,3)          ==  (0,[],[3])
     decimal (23,14)        ==  (1,[6],[4,2,8,5,7,1])
     decimal (247813,19980) ==  (12,[4,0],[3,0,5])
     decimal (1,101)        ==  (0,[],[0,0,9,9])

decimal (6,2) == (3,[],[])

decimal (3,4) == (0,[7,5],[])

decimal (1,3) == (0,[],[3])

decimal (23,14) == (1,[6],[4,2,8,5,7,1])

decimal (247813,19980) == (12,[4,0],[3,0,5])

decimal (1,101) == (0,[],[0,0,9,9])

(longitudPeriodo f) es la longitud de la parte periódica de la representación decimal de la fracción f. Por ejemplo,

     longitudPeriodo (6,2)           ==  0
     longitudPeriodo (3,4)           ==  0
     longitudPeriodo (1,3)           ==  1
     longitudPeriodo (23,14)         ==  6
     longitudPeriodo (247813,19980)  ==  3
     longitudPeriodo (1,101)         ==  4
     longitudPeriodo (1,1229)        ==  1228

longitudPeriodo (6,2) == 0

longitudPeriodo (3,4) == 0

longitudPeriodo (1,3) == 1

longitudPeriodo (23,14) == 6

longitudPeriodo (247813,19980) == 3

longitudPeriodo (1,101) == 4

longitudPeriodo (1,1229) == 1228

Comprobar con QuickCheck la propiedad de la longitud de la parte periódica; es decir, k es un número natural distinto de 0 y 2 y p es el primo k-ésimo, entonces la longitud del período de 1/p es el menor entero positivo n tal que p divide a $10^n - 1$ ..

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

type Fraccion = (Integer,Integer)
type Decimal = (Integer,[Integer],[Integer])

decimal :: Fraccion -> Decimal
decimal (n,d) 
  | snd y == 0 = (fst x, map fst xs, [])
  | otherwise  = (fst x, map fst xs, map fst (y:zs))
  where
    qrs         = cocientesRestos (n,d)
    Just (q,r)  = primerRepetido qrs
    (x:xs,y:ys) = break (==(q,r)) qrs
    zs          = takeWhile (/=(q,r)) ys

cocientesRestos :: Fraccion -> [(Integer,Integer)]
cocientesRestos (n,d) =
  (q,r) : cocientesRestos (10*r, d)
  where (q,r) = quotRem n d

primerRepetido :: Eq a => [a] -> Maybe a
primerRepetido xs = aux xs []
  where
    aux [] _                     = Nothing
    aux (x:xs') ys | x `elem` ys = Just x
                   | otherwise   = aux xs' (x:ys) 

longitudPeriodo :: Fraccion -> Int
longitudPeriodo (n,d) = length xs
  where (_,_,xs) = decimal (n,d)

-- La propiedad es
prop_LongitudPeriodo :: Int -> Property
prop_LongitudPeriodo k =
  k > 0 && k /= 2 
  ==>
  longitudPeriodo (1,p) ==
  head [n | n <- [1..], (10^n-1) `mod` p == 0]
  where p = primes !! k

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

import Data.Numbers.Primes

import Test.QuickCheck

type Fraccion = (Integer,Integer)

type Decimal = (Integer,[Integer],[Integer])

decimal :: Fraccion -> Decimal

decimal (n,d)

| snd y == 0 = (fst x, map fst xs, [])

| otherwise = (fst x, map fst xs, map fst (y:zs))

where

qrs = cocientesRestos (n,d)

Just (q,r) = primerRepetido qrs

(x:xs,y:ys) = break (==(q,r)) qrs

zs = takeWhile (/=(q,r)) ys

cocientesRestos :: Fraccion -> [(Integer,Integer)]

cocientesRestos (n,d) =

(q,r) : cocientesRestos (10*r, d)

where (q,r) = quotRem n d

primerRepetido :: Eq a => [a] -> Maybe a

primerRepetido xs = aux xs []

where

aux [] _ = Nothing

aux (x:xs') ys | x `elem` ys = Just x

| otherwise = aux xs' (x:ys)

longitudPeriodo :: Fraccion -> Int

longitudPeriodo (n,d) = length xs

where (_,_,xs) = decimal (n,d)

-- La propiedad es

prop_LongitudPeriodo :: Int -> Property

prop_LongitudPeriodo k =

k > 0 && k /= 2

==>

longitudPeriodo (1,p) ==

head [n | n <- [1..], (10^n-1) `mod` p == 0]

where p = primes !! k

-- La comprobación es

-- λ> quickCheck prop_LongitudPeriodo

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

«En el desarrollo de la comprensión de los fenómenos complejos, la herramienta más poderosa de que dispone el intelecto humano es la abstracción. La abstracción surge del reconocimiento de las similitudes entre ciertos objetos, situaciones o procesos en el mundo real y de la decisión de concentrarse en estas similitudes e ignorar, por el momento, sus diferencias.»

Tony Hoare

Elementos múltiplos de la longitud de la lista

PorJosé A. Alonso 16-enero-202023-enero-2020

Definir las funciones

   multiplosDeLaLongitud :: [Int] -> [Int]
   multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

1 2	multiplosDeLaLongitud :: [Int] -> [Int] multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

tales que

(multiplosDeLaLongitud xs) es la lista de los elementos de xs que son múltiplos de la longitud de xs. Por ejemplo,

     multiplosDeLaLongitud [2,4,6,8] == [4,8]

1	multiplosDeLaLongitud [2,4,6,8] == [4,8]

(multiplosDeLaLongitudDeConsecutivos n m) es la lista de elementos de [n..n+m-1] que son múltiplos de n. Por ejemplo,

     multiplosDeLaLongitudDeConsecutivos 9 2  ==  [10]
     multiplosDeLaLongitudDeConsecutivos 9 12 ==  [12]

1 2	multiplosDeLaLongitudDeConsecutivos 9 2 == [10] multiplosDeLaLongitudDeConsecutivos 9 12 == [12]

Comprobar con QuickCheck si se verifican las siguientes propiedades

En cualquier conjunto de m elementos consecutivos, m divide exactamente a uno de dichos elementos. En otras palabras, si n y m son enteros positivos, entonces (multiplosDeLaLongitudDeConsecutivos n m) tiene exactamente un elemento.
Si n es un entero positivo y m >= n, entonces (multiplosDeLaLongitudDeConsecutivos n m) es igual a [m]
Si n y n son enteros positivos y m < n, entonces (multiplosDeLaLongitudDeConsecutivos n m) es igual a [m * ceiling (n’ / m’)] donde n’ y m’ son las formas decimales de n y m respectivamente.

Soluciones

import Test.QuickCheck

multiplosDeLaLongitud :: [Int] -> [Int]
multiplosDeLaLongitud xs =
  [x | x <- xs
     , x `mod`  n == 0]
  where n = length xs

multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]
multiplosDeLaLongitudDeConsecutivos n m =
  multiplosDeLaLongitud [n..n+m-1]

-- La 1ª propiedad es
prop_multiplosDeLaLongitud :: Int -> Int -> Property
prop_multiplosDeLaLongitud n m =
  n > 0 && m > 0
  ==>
  length (multiplosDeLaLongitudDeConsecutivos n m) == 1

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

-- La 2ª propiedad es
prop_multiplosDeLaLongitud2 :: Int -> Int -> Property
prop_multiplosDeLaLongitud2 n m =
  n > 0 && m >= n
  ==>
  multiplosDeLaLongitudDeConsecutivos n m == [m]

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

-- La 3ª propiedad es
prop_multiplosDeLaLongitud3 :: Int -> Int -> Property
prop_multiplosDeLaLongitud3 n m =
  n > 0 && 0 < m && m < n
  ==>
  multiplosDeLaLongitudDeConsecutivos n m == [m * ceiling (n' / m')]
  where n' = fromIntegral n
        m' = fromIntegral m 
  
-- La comprobación es
--    λ> quickCheck prop_multiplosDeLaLongitud3
--    +++ OK, passed 100 tests.

-- Con las propiedades anteriores se puede redefinir multiplos
multiplosDeLaLongitudDeConsecutivos2 :: Int -> Int -> [Int]
multiplosDeLaLongitudDeConsecutivos2 n m
  | m < n     = [m * ceiling (n' / m')]
  | otherwise = [m]
  where n' = fromIntegral n
        m' = fromIntegral m 

-- Comprobación de la equivalencia
prop_equiv :: Int -> Int -> Property
prop_equiv n m =
  n > 0 && m > 0
  ==>
  multiplosDeLaLongitudDeConsecutivos  n m ==
  multiplosDeLaLongitudDeConsecutivos2 n m

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

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

--    λ> xs1 = [multiplosDeLaLongitudDeConsecutivos n (10^3) | n <- [1..10^4]]
--    (0.01 secs, 0 bytes)
--    λ> xs2 = [multiplosDeLaLongitudDeConsecutivos2 n (10^3) | n <- [1..10^4]]
--    (0.01 secs, 0 bytes)
--    λ> maximum xs1
--    [10000]
--    (2.51 secs, 2,093,241,168 bytes)
--    λ> maximum xs2
--    [10000]
--    (0.02 secs, 16,742,408 bytes)

import Test.QuickCheck

multiplosDeLaLongitud :: [Int] -> [Int]

multiplosDeLaLongitud xs =

[x | x <- xs

, x `mod` n == 0]

where n = length xs

multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

multiplosDeLaLongitudDeConsecutivos n m =

multiplosDeLaLongitud [n..n+m-1]

-- La 1ª propiedad es

prop_multiplosDeLaLongitud :: Int -> Int -> Property

prop_multiplosDeLaLongitud n m =

n > 0 && m > 0

==>

length (multiplosDeLaLongitudDeConsecutivos n m) == 1

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud

-- +++ OK, passed 100 tests.

-- La 2ª propiedad es

prop_multiplosDeLaLongitud2 :: Int -> Int -> Property

prop_multiplosDeLaLongitud2 n m =

n > 0 && m >= n

==>

multiplosDeLaLongitudDeConsecutivos n m == [m]

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud2

-- +++ OK, passed 100 tests.

-- La 3ª propiedad es

prop_multiplosDeLaLongitud3 :: Int -> Int -> Property

prop_multiplosDeLaLongitud3 n m =

n > 0 && 0 < m && m < n

==>

multiplosDeLaLongitudDeConsecutivos n m == [m * ceiling (n' / m')]

where n' = fromIntegral n

m' = fromIntegral m

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud3

-- +++ OK, passed 100 tests.

-- Con las propiedades anteriores se puede redefinir multiplos

multiplosDeLaLongitudDeConsecutivos2 :: Int -> Int -> [Int]

multiplosDeLaLongitudDeConsecutivos2 n m

| m < n = [m * ceiling (n' / m')]

| otherwise = [m]

where n' = fromIntegral n

m' = fromIntegral m

-- Comprobación de la equivalencia

prop_equiv :: Int -> Int -> Property

prop_equiv n m =

n > 0 && m > 0

==>

multiplosDeLaLongitudDeConsecutivos n m ==

multiplosDeLaLongitudDeConsecutivos2 n m

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> xs1 = [multiplosDeLaLongitudDeConsecutivos n (10^3) | n <- [1..10^4]]

-- (0.01 secs, 0 bytes)

-- λ> xs2 = [multiplosDeLaLongitudDeConsecutivos2 n (10^3) | n <- [1..10^4]]

-- (0.01 secs, 0 bytes)

-- λ> maximum xs1

-- [10000]

-- (2.51 secs, 2,093,241,168 bytes)

-- λ> maximum xs2

-- [10000]

-- (0.02 secs, 16,742,408 bytes)

Referencia

Set of successive numbers contains unique multiple en ProofWiki.

Pensamiento

Pensando que no veía
porque Dios no le miraba,
dijo Abel cuando moría:
Se acabó lo que se daba.

Antonio Machado

Teorema de existencia de divisores

PorJosé A. Alonso 10-enero-20206-junio-2022

El teorema de existencia de divisores afirma que

En cualquier subconjunto de {1, 2, …, 2m} con al menos m+1 elementos existen números distintos a, b tales que a divide a b.

Un conjunto de números naturales xs es mayoritario si existe un m tal que la lista de xs es un subconjunto de {1,2,…,2m} con al menos m+1 elementos. Por ejemplo, {2,3,5,6} porque es un subconjunto de {1,2,…,6} con más de 3 elementos.

Definir las funciones

   divisoresMultiplos :: [Integer] -> [(Integer,Integer)]
   esMayoritario :: [Integer] -> Bool

1 2	divisoresMultiplos :: [Integer] -> [(Integer,Integer)] esMayoritario :: [Integer] -> Bool

tales que

(divisores xs) es la lista de pares de elementos distintos de (a,b) tales que a divide a b. Por ejemplo,

     divisoresMultiplos [2,3,5,6]  ==  [(2,6),(3,6)]
     divisoresMultiplos [2,3,5]    ==  []
     divisoresMultiplos [4..8]     ==  [(4,8)]

divisoresMultiplos [2,3,5,6] == [(2,6),(3,6)]

divisoresMultiplos [2,3,5] == []

divisoresMultiplos [4..8] == [(4,8)]

(esMayoritario xs) se verifica xs es mayoritario. Por ejemplo,

     esMayoritario [2,3,5,6]  ==  True
     esMayoritario [2,3,5]    ==  False

1 2	esMayoritario [2,3,5,6] == True esMayoritario [2,3,5] == False

Comprobar con QuickCheck el teorema de existencia de divisores; es decir, en cualquier conjunto mayoritario existen números distintos a, b tales que a divide a b. Para la comprobación se puede usar el siguiente generador de conjuntos mayoritarios

   mayoritario :: Gen [Integer]
   mayoritario = do
     m' <- arbitrary
     let m = 1 + abs m'
     xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)
     return xs'

mayoritario :: Gen [Integer]

mayoritario = do

m' <- arbitrary

let m = 1 + abs m'

xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)

return xs'

con lo que la propiedad que hay que comprobar con QuickCheck es

   teorema_de_existencia_de_divisores :: Property
   teorema_de_existencia_de_divisores =
     forAll mayoritario (not . null . divisoresMultiplos)

teorema_de_existencia_de_divisores :: Property

teorema_de_existencia_de_divisores =

forAll mayoritario (not . null . divisoresMultiplos)

Soluciones

import Data.List (genericLength)
import Test.QuickCheck

divisoresMultiplos :: [Integer] -> [(Integer,Integer)]
divisoresMultiplos xs =
  [(x,y) | x <- xs
         , y <- xs
         , y /= x
         , y `mod` x == 0]

esMayoritario :: [Integer] -> Bool
esMayoritario xs =
  not (null xs) && length xs > ceiling (n / 2) 
  where n = fromIntegral (maximum xs)

-- Comprobación del teorema
-- ========================

-- La propiedad es
teorema_de_existencia_de_divisores :: Property
teorema_de_existencia_de_divisores =
  forAll mayoritario (not . null . divisoresMultiplos)

-- mayoritario es un generador de conjuntos mayoritarios. Por ejemplo, 
--    λ> sample mayoritario
--    [1,2]
--    [2,5,7,8]
--    [1,2,8,10,14]
--    [3,8,11,12,13,15,18,19,22,23,25,26]
--    [1,3,4,6]
--    [3,6,9,11,12,14,17,19]
mayoritario :: Gen [Integer]
mayoritario = do
  m' <- arbitrary
  let m = 1 + abs m'
  xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)
  return xs'

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

import Data.List (genericLength)

import Test.QuickCheck

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

divisoresMultiplos xs =

[(x,y) | x <- xs

, y <- xs

, y /= x

, y `mod` x == 0]

esMayoritario :: [Integer] -> Bool

esMayoritario xs =

not (null xs) && length xs > ceiling (n / 2)

where n = fromIntegral (maximum xs)

-- Comprobación del teorema

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

-- La propiedad es

teorema_de_existencia_de_divisores :: Property

teorema_de_existencia_de_divisores =

forAll mayoritario (not . null . divisoresMultiplos)

-- mayoritario es un generador de conjuntos mayoritarios. Por ejemplo,

-- λ> sample mayoritario

-- [1,2]

-- [2,5,7,8]

-- [1,2,8,10,14]

-- [3,8,11,12,13,15,18,19,22,23,25,26]

-- [1,3,4,6]

-- [3,6,9,11,12,14,17,19]

mayoritario :: Gen [Integer]

mayoritario = do

m' <- arbitrary

let m = 1 + abs m'

xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)

return xs'

-- La comprobación es

-- λ> quickCheck teorema_de_existencia_de_divisores

-- +++ OK, passed 100 tests.

Pensamiento

Guiomar, Guiomar,
mírame en ti castigado:
reo de haberte creado,
ya no te puedo olvidar.

Antonio Machado

La mayor potencia de dos no es divisor

PorJosé A. Alonso 7-enero-202014-enero-2020

Para cada número entero positivo n, se define el conjunto

   S(n) = {1, 2, 3, ..., n}

1	S(n) = {1, 2, 3, ..., n}

de los números desde el 1 hasta n.

Definir la función

   mayorPotenciaDeDosEnS :: Integer -> Integer

1	mayorPotenciaDeDosEnS :: Integer -> Integer

tal que (mayorPotenciaDeDosEnS n) es la mayor potencia de 2 en S(n). Por ejemplo,

   mayorPotenciaDeDosEnS 3  ==  2
   mayorPotenciaDeDosEnS 4  ==  4

1 2	mayorPotenciaDeDosEnS 3 == 2 mayorPotenciaDeDosEnS 4 == 4

Comprobar con QuickCheck que la mayor potencia de 2 en S(n) no divide a ningún otro elemento de S(n).

Soluciones

import Data.List (delete)
import Test.QuickCheck

mayorPotenciaDeDosEnS :: Integer -> Integer
mayorPotenciaDeDosEnS n =
  last (takeWhile (<=n) potenciasDeDos)

-- potenciasDeDos es la lista de las potencias de 2. Por ejemplo,
--    take 5 potenciasDeDos  ==  [1,2,4,8,16]
potenciasDeDos :: [Integer]
potenciasDeDos = iterate (*2) 1

-- La propiedad es
prop_mayorPotenciaDeDosEnS :: Integer -> Property
prop_mayorPotenciaDeDosEnS n =
  n > 0 ==> all (x `noDivide`) (delete x [1..n])
  where x = mayorPotenciaDeDosEnS n
        x `noDivide` y = y `mod` x /= 0

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

import Data.List (delete)

import Test.QuickCheck

mayorPotenciaDeDosEnS :: Integer -> Integer

mayorPotenciaDeDosEnS n =

last (takeWhile (<=n) potenciasDeDos)

-- potenciasDeDos es la lista de las potencias de 2. Por ejemplo,

-- take 5 potenciasDeDos == [1,2,4,8,16]

potenciasDeDos :: [Integer]

potenciasDeDos = iterate (*2) 1

-- La propiedad es

prop_mayorPotenciaDeDosEnS :: Integer -> Property

prop_mayorPotenciaDeDosEnS n =

n > 0 ==> all (x `noDivide`) (delete x [1..n])

where x = mayorPotenciaDeDosEnS n

x `noDivide` y = y `mod` x /= 0

-- La comprobación es

-- λ> quickCheck prop_mayorPotenciaDeDosEnS

-- +++ OK, passed 100 tests.

Referencia

Basado en Greatest power of two not divisor de ProofWiki.

Pensamiento

¡Sólo tu figura,
como una centella blanca,
en mi noche oscura.

Antonio Machado

Factorizaciones de 4n+1

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

Sea S el conjunto

   S = {1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, ...}

1	S = {1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, ...}

de los enteros positivos congruentes con 1 módulo 4; es decir,

   S = {4n+1 | n ∈ N}

1	S = {4n+1 \| n ∈ N}

Un elemento n de S es irreducible si sólo es divisible por dos elementos de S: 1 y n. Por ejemplo, 9 es irreducible; pero 45 no lo es (ya que es el proctos de 5 y 9 que son elementos de S).

Definir las funciones

   esIrreducible :: Integer -> Bool
   factorizaciones :: Integer -> [[Integer]]
   conFactorizacionNoUnica :: [Integer]

esIrreducible :: Integer -> Bool

factorizaciones :: Integer -> [[Integer]]

conFactorizacionNoUnica :: [Integer]

tales que

(esIrreducible n) se verifica si n es irreducible. Por ejemplo,

     esIrreducible 9   ==  True
     esIrreducible 45  ==  False

1 2	esIrreducible 9 == True esIrreducible 45 == False

(factorizaciones n) es la lista de conjuntos de elementos irreducibles de S cuyo producto es n. Por ejemplo,

     factorizaciones 9     ==  [[9]]
     factorizaciones 693   ==  [[9,77],[21,33]]
     factorizaciones 4389  ==  [[21,209],[33,133],[57,77]]
     factorizaciones 2205  ==  [[5,9,49],[5,21,21]]

factorizaciones 9 == [[9]]

factorizaciones 693 == [[9,77],[21,33]]

factorizaciones 4389 == [[21,209],[33,133],[57,77]]

factorizaciones 2205 == [[5,9,49],[5,21,21]]

conFactorizacionNoUnica es la lista de elementos de S cuya factorización no es única. Por ejemplo,

     λ> take 10 conFactorizacionNoUnica
     [441,693,1089,1197,1449,1617,1881,1953,2205,2277]

1 2	λ> take 10 conFactorizacionNoUnica [441,693,1089,1197,1449,1617,1881,1953,2205,2277]

Soluciones

esIrreducible :: Integer -> Bool
esIrreducible n = divisores n == [1,n]

-- (divisores n) es el conjunto de elementos de S que dividen a n. Por
-- ejemplo, 
--    divisores 9    ==  [9]
--    divisores 693  ==  [9,21,33,77,693]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1,5..n], n `mod` x == 0] 
  
factorizaciones :: Integer -> [[Integer]]
factorizaciones 1 = [[1]]
factorizaciones n
  | esIrreducible n = [[n]]
  | otherwise       = [d:ds | d <- divisores n
                            , esIrreducible d
                            , ds@(e:_) <- factorizaciones (n `div` d)
                            , d <= e]

conFactorizacionNoUnica :: [Integer]
conFactorizacionNoUnica =
  [n | n <- [1,5..], length (factorizaciones n) > 1]

esIrreducible :: Integer -> Bool

esIrreducible n = divisores n == [1,n]

-- (divisores n) es el conjunto de elementos de S que dividen a n. Por

-- ejemplo,

-- divisores 9 == [9]

-- divisores 693 == [9,21,33,77,693]

divisores :: Integer -> [Integer]

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

factorizaciones :: Integer -> [[Integer]]

factorizaciones 1 = [[1]]

factorizaciones n

| esIrreducible n = [[n]]

| otherwise = [d:ds | d <- divisores n

, esIrreducible d

, ds@(e:_) <- factorizaciones (n `div` d)

, d <= e]

conFactorizacionNoUnica :: [Integer]

conFactorizacionNoUnica =

[n | n <- [1,5..], length (factorizaciones n) > 1]

Pensamiento

¡Qué bien los nombres ponía
quien puso Sierra Morena
a esta serranía!

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

Teorema de Carmichael

PorJosé A. Alonso 1-enero-20208-enero-2020

La sucesión de Fibonacci, F(n), es la siguiente sucesión infinita de números naturales:

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

La sucesión comieanza con los números 0 y 1. A partir de estos, cada término es la suma de los dos anteriores.

El teorema de Carmichael establece que para todo n mayor que 12, el n-ésimo número de Fibonacci F(n) tiene al menos un factor primo que no es factor de ninguno de los términos anteriores de la sucesión.

Si un número primo p es un factor de F(n) y no es factor de ningún F(m) con m < n, entonces se dice que p es un factor característico o un divisor primitivo de F(n).

Definir la función

   factoresCaracteristicos :: Int -> [Integer]

1	factoresCaracteristicos :: Int -> [Integer]

tal que (factoresCaracteristicos n) es la lista de los factores característicos de F(n). Por ejemplo,

   factoresCaracteristicos  4  ==  [3]
   factoresCaracteristicos  6  ==  []
   factoresCaracteristicos 19  ==  [37,113]
   factoresCaracteristicos 20  ==  [41]
   factoresCaracteristicos 37  ==  [73,149,2221]

factoresCaracteristicos 4 == [3]

factoresCaracteristicos 6 == []

factoresCaracteristicos 19 == [37,113]

factoresCaracteristicos 20 == [41]

factoresCaracteristicos 37 == [73,149,2221]

Comprobar con QuickCheck el teorema de Carmichael; es decir, para todo número entero (factoresCaracteristicos (13 + abs n)) es una lista no vacía.

Soluciones

import Data.List (nub)
import Data.Numbers.Primes
import Test.QuickCheck

factoresCaracteristicos :: Int -> [Integer]
factoresCaracteristicos n =
  [x | x <- factoresPrimos (fib n)
     , and [fib m `mod` x /= 0 | m <- [1..n-1]]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por
-- ejemplo,
--    fib 6  ==  8
fib :: Int -> Integer
fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,
--    λ> take 20 fibs
--    [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

-- (factoresPrimos n) es la lista de los factores primos de n. Por
-- ejemplo, 
--    factoresPrimos 600  ==  [2,3,5]
factoresPrimos :: Integer -> [Integer]
factoresPrimos 0 = []
factoresPrimos n = nub (primeFactors n)

-- Teorema
-- =======

-- El teorema es
teorema_de_Carmichael :: Int -> Bool
teorema_de_Carmichael n =
  not (null (factoresCaracteristicos n'))
  where n' = 13 + abs n

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

import Data.List (nub)

import Data.Numbers.Primes

import Test.QuickCheck

factoresCaracteristicos :: Int -> [Integer]

factoresCaracteristicos n =

[x | x <- factoresPrimos (fib n)

, and [fib m `mod` x /= 0 | m <- [1..n-1]]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por

-- ejemplo,

-- fib 6 == 8

fib :: Int -> Integer

fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,

-- λ> take 20 fibs

-- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]

fibs :: [Integer]

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

-- (factoresPrimos n) es la lista de los factores primos de n. Por

-- ejemplo,

-- factoresPrimos 600 == [2,3,5]

factoresPrimos :: Integer -> [Integer]

factoresPrimos 0 = []

factoresPrimos n = nub (primeFactors n)

-- Teorema

-- =======

-- El teorema es

teorema_de_Carmichael :: Int -> Bool

teorema_de_Carmichael n =

not (null (factoresCaracteristicos n'))

where n' = 13 + abs n

-- La comprobación es

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

-- +++ OK, passed 100 tests.

Pensamiento

No puede ser
amor de tanta fortuna:
dos soledades en una.

Antonio Machado

El teorema de Navidad de Fermat

PorJosé A. Alonso 20-diciembre-201927-diciembre-2019

El 25 de diciembre de 1640, en una carta a Mersenne, Fermat demostró la conjetura de Girard: todo primo de la forma 4n+1 puede expresarse de manera única como suma de dos cuadrados. Por eso es conocido como el Teorema de Navidad de Fermat.

Definir las funciones

   representaciones :: Integer -> [(Integer,Integer)]
   primosImparesConRepresentacionUnica :: [Integer]
   primos4nM1 :: [Integer]

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

(representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2 con x <= y. Por ejemplo,

     representaciones  20           ==  [(2,4)]
     representaciones  25           ==  [(0,5),(3,4)]
     representaciones 325           ==  [(1,18),(6,17),(10,15)]
     representaciones 100000147984  ==  [(0,316228)]
     length (representaciones (10^10))    ==  6
     length (representaciones (4*10^12))  ==  7

representaciones 20 == [(2,4)]

representaciones 25 == [(0,5),(3,4)]

representaciones 325 == [(1,18),(6,17),(10,15)]

representaciones 100000147984 == [(0,316228)]

length (representaciones (10^10)) == 6

length (representaciones (4*10^12)) == 7

primosImparesConRepresentacionUnica es la lista de los números primos impares que se pueden escribir exactamente de una manera como suma de cuadrados de pares de números naturales (x,y) con x <= y. Por ejemplo,

     λ> take 20 primosImparesConRepresentacionUnica
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primosImparesConRepresentacionUnica [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

primos4nM1 es la lista de los números primos que se pueden escribir como uno más un múltiplo de 4 (es decir, que son congruentes con 1 módulo 4). Por ejemplo,

     λ> take 20 primos4nM1
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primos4nM1 [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

El teorema de Navidad de Fermat afirma que un número primo impar p se puede escribir exactamente de una manera como suma de dos cuadrados de números naturales p = x² + y^2 (con x <= y) si, y sólo si, p se puede escribir como uno más un múltiplo de 4 (es decir, que es congruente con 1 módulo 4).

Comprobar con QuickCheck el teorema de Navidad de Fermat; es decir, que para todo número n, los n-ésimos elementos de primosImparesConRepresentacionUnica y de primos4nM1 son iguales.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck

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

representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

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

representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)] 
              , let z = n - x*x
              , esCuadrado z]

-- (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 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2

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

representaciones3 :: Integer -> [(Integer,Integer)]
representaciones3 n =
  [(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)] 
               , let z = n - x*x
               , esCuadrado3 z]

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

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz3 25  ==  5
--    raiz3 24  ==  4
--    raiz3 26  ==  5
raiz3 :: Integer -> Integer
raiz3 x = floor (sqrt (fromIntegral x))
          
-- 4ª definición de representaciones
-- =================================
 
representaciones4 :: Integer -> [(Integer, Integer)]
representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))
  where aux x y
          | x > y     = [] 
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)

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

-- La propiedad es
prop_representaciones_equiv :: (Positive Integer) -> Bool
prop_representaciones_equiv (Positive n) =
  representaciones  n == representaciones2 n &&
  representaciones2 n == representaciones3 n &&
  representaciones3 n == representaciones4 n

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

-- Comparación de eficiencia de las definiciones de representaciones
-- =================================================================

--    λ> representaciones 3025
--    [(0,55),(33,44)]
--    (2.86 secs, 1,393,133,528 bytes)
--    λ> representaciones2 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 867,944 bytes)
--    λ> representaciones3 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 173,512 bytes)
--    λ> representaciones4 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 423,424 bytes)
--    
--    λ> length (representaciones2 (10^10))
--    6
--    (3.38 secs, 2,188,903,544 bytes)
--    λ> length (representaciones3 (10^10))
--    6
--    (0.10 secs, 62,349,048 bytes)
--    λ> length (representaciones4 (10^10))
--    6
--    (0.11 secs, 48,052,360 bytes)
--
--    λ> length (representaciones3 (4*10^12))
--    7
--    (1.85 secs, 1,222,007,176 bytes)
--    λ> length (representaciones4 (4*10^12))
--    7
--    (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica
-- =================================================

primosImparesConRepresentacionUnica :: [Integer]
primosImparesConRepresentacionUnica =
  [x | x <- tail primes
     , length (representaciones4 x) == 1]

-- Definición de primos4nM1
-- ========================

primos4nM1 :: [Integer]
primos4nM1 = [x | x <- primes
                , x `mod` 4 == 1]

-- Teorema de Navidad de Fermat
-- ============================

-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
             
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ 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

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª definición de representaciones

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

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

representaciones n =

[(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

-- 2ª definición de representaciones

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

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

representaciones2 n =

[(x,raiz z) | x <- [0..raiz (n `div` 2)]

, let z = n - x*x

, esCuadrado z]

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

raiz 1 = 1

raiz x = aux (0,x)

where aux (a,b) | d == x = c

| c == a = a

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^2

-- 3ª definición de representaciones

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

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

representaciones3 n =

[(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)]

, let z = n - x*x

, esCuadrado3 z]

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

-- ejemplo,

-- esCuadrado3 25 == True

-- esCuadrado3 26 == False

esCuadrado3 :: Integer -> Bool

esCuadrado3 x = x == y * y

where y = raiz3 x

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

-- raiz3 25 == 5

-- raiz3 24 == 4

-- raiz3 26 == 5

raiz3 :: Integer -> Integer

raiz3 x = floor (sqrt (fromIntegral x))

-- 4ª definición de representaciones

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

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

representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))

where aux x y

| x > y = []

| otherwise = case compare (x*x + y*y) n of

LT -> aux (x + 1) y

EQ -> (x, y) : aux (x + 1) (y - 1)

GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones

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

-- La propiedad es

prop_representaciones_equiv :: (Positive Integer) -> Bool

prop_representaciones_equiv (Positive n) =

representaciones n == representaciones2 n &&

representaciones2 n == representaciones3 n &&

representaciones3 n == representaciones4 n

-- La comprobación es

-- λ> quickCheck prop_representaciones_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de las definiciones de representaciones

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

-- λ> representaciones 3025

-- [(0,55),(33,44)]

-- (2.86 secs, 1,393,133,528 bytes)

-- λ> representaciones2 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 867,944 bytes)

-- λ> representaciones3 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 173,512 bytes)

-- λ> representaciones4 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 423,424 bytes)

-- λ> length (representaciones2 (10^10))

-- 6

-- (3.38 secs, 2,188,903,544 bytes)

-- λ> length (representaciones3 (10^10))

-- 6

-- (0.10 secs, 62,349,048 bytes)

-- λ> length (representaciones4 (10^10))

-- 6

-- (0.11 secs, 48,052,360 bytes)

-- λ> length (representaciones3 (4*10^12))

-- 7

-- (1.85 secs, 1,222,007,176 bytes)

-- λ> length (representaciones4 (4*10^12))

-- 7

-- (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica

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

primosImparesConRepresentacionUnica :: [Integer]

primosImparesConRepresentacionUnica =

[x | x <- tail primes

, length (representaciones4 x) == 1]

-- Definición de primos4nM1

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

primos4nM1 :: [Integer]

primos4nM1 = [x | x <- primes

, x `mod` 4 == 1]

-- Teorema de Navidad de Fermat

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

-- La propiedad es

prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool

prop_teoremaDeNavidadDeFermat (Positive n) =

primosImparesConRepresentacionUnica !! n == primos4nM1 !! n

-- La comprobación es

-- λ> quickCheck prop_teoremaDeNavidadDeFermat

-- +++ OK, passed 100 tests.

Pensamiento

Dijo Dios: brote la nada
Y alzó su mano derecha,
hasta ocultar su mirada.
Y quedó la nada hecha.

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

Suma de primos menores

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

La suma de los primos menores que 10 es 2 + 3 + 5 + 7 = 17.

Definir la función

   sumaPrimosMenores :: Integer -> Integer

1	sumaPrimosMenores :: Integer -> Integer

tal que (sumaPrimosMenores n) es la suma de los primos menores que n. Por ejemplo,

   sumaPrimosMenores 10        ==  17
   sumaPrimosMenores (5*10^5)  ==  9914236195

1 2	sumaPrimosMenores 10 == 17 sumaPrimosMenores (5*10^5) == 9914236195

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

Soluciones

import Data.Numbers.Primes (primes)

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

sumaPrimosMenores :: Integer -> Integer
sumaPrimosMenores n = sum (takeWhile (<n) primos)

-- primos es la lista de los números primos. Por ejemplo,
--    λ> take 12 primos
--    [2,3,5,7,11,13,17,19,23,29,31,37]
primos :: [Integer]
primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool
esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]
                    , n `mod` x == 0]

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

sumaPrimosMenores2 :: Integer -> Integer
sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]
primos2 = 2 : filter esPrimo2 [3,5..]
 
esPrimo2 :: Integer -> Bool
esPrimo2 x =
  all ((/= 0) . mod x)
  (takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

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

sumaPrimosMenores3 :: Integer -> Integer
sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

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

--    λ> sumaPrimosMenores (2*10^5)
--    1709600813
--    (2.56 secs, 1,522,015,240 bytes)
--    λ> sumaPrimosMenores2 (2*10^5)
--    1709600813
--    (0.56 secs, 376,951,456 bytes)
--    λ> sumaPrimosMenores3 (2*10^5)
--    1709600813
--    (0.07 secs, 62,321,888 bytes)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

sumaPrimosMenores :: Integer -> Integer

sumaPrimosMenores n = sum (takeWhile (<n) primos)

-- primos es la lista de los números primos. Por ejemplo,

-- λ> take 12 primos

-- [2,3,5,7,11,13,17,19,23,29,31,37]

primos :: [Integer]

primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool

esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]

, n `mod` x == 0]

-- 2ª solución

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

sumaPrimosMenores2 :: Integer -> Integer

sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]

primos2 = 2 : filter esPrimo2 [3,5..]

esPrimo2 :: Integer -> Bool

esPrimo2 x =

all ((/= 0) . mod x)

(takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

-- 3ª solución

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

sumaPrimosMenores3 :: Integer -> Integer

sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

-- Comparación de eficiencia

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

-- λ> sumaPrimosMenores (2*10^5)

-- 1709600813

-- (2.56 secs, 1,522,015,240 bytes)

-- λ> sumaPrimosMenores2 (2*10^5)

-- 1709600813

-- (0.56 secs, 376,951,456 bytes)

-- λ> sumaPrimosMenores3 (2*10^5)

-- 1709600813

-- (0.07 secs, 62,321,888 bytes)

Pensamiento

El movimiento no es nada esencial. La fuerza puede ser inmóvil (lo es en su estado de pureza); mas no por ello deja de ser activa.

Antonio Machado

Menor divisible por todos

PorJosé A. Alonso 28-noviembre-20195-diciembre-2019

Definir la función

   menorDivisible :: Integer -> Integer -> Integer

1	menorDivisible :: Integer -> Integer -> Integer

tal que (menorDivisible a b) es el menor número divisible por todos los números desde a hasta b, ambos inclusive. Por ejemplo,

   menorDivisible 1 10                        ==  2520
   length (show (menorDivisible 1 (3*10^5)))  ==  130141

1 2	menorDivisible 1 10 == 2520 length (show (menorDivisible 1 (3*10^5))) == 130141

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

Soluciones

import Data.List (foldl')

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

menorDivisible :: Integer -> Integer -> Integer
menorDivisible a b =
  head [x | x <- [b..]
          , and [x `mod` y == 0 | y <- [a..b]]]

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

menorDivisible2 :: Integer -> Integer -> Integer
menorDivisible2 a b  
  | a == b    = a
  | otherwise = lcm a (menorDivisible (a+1) b)

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

menorDivisible3 :: Integer -> Integer -> Integer
menorDivisible3 a b = foldl lcm 1 [a..b] 

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

menorDivisible4 :: Integer -> Integer -> Integer
menorDivisible4 a b = foldl1 lcm [a..b] 

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

menorDivisible5 :: Integer -> Integer -> Integer
menorDivisible5 a b = foldl' lcm 1 [a..b] 

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

menorDivisible6 :: Integer -> Integer -> Integer
menorDivisible6 a b = foldr1 lcm [a..b] 

-- 7ª solución
-- ===========

menorDivisible7 :: Integer -> Integer -> Integer
menorDivisible7 a = foldr1 lcm . enumFromTo a

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

--   λ> menorDivisible 1 17
--   12252240
--   (18.63 secs, 15,789,475,488 bytes)
--   λ> menorDivisible2 1 17
--   12252240
--   (13.29 secs, 11,868,764,272 bytes)
--   λ> menorDivisible3 1 17
--   12252240
--   (0.00 secs, 114,688 bytes)
--   λ> menorDivisible4 1 17
--   12252240
--   (0.01 secs, 114,752 bytes)
--   λ> menorDivisible5 1 17
--   12252240
--   (0.01 secs, 110,640 bytes)
--   λ> menorDivisible6 1 17
--   12252240
--   (0.01 secs, 114,752 bytes)
--   λ> menorDivisible7 1 17
--   12252240
--   (0.00 secs, 110,912 bytes)
--   
--   λ> length (show (menorDivisible3 1 (10^5)))
--   43452
--   (1.54 secs, 2,021,887,000 bytes)
--   λ> length (show (menorDivisible4 1 (10^5)))
--   43452
--   (1.47 secs, 2,021,886,616 bytes)
--   λ> length (show (menorDivisible5 1 (10^5)))
--   43452
--   (0.65 secs, 2,009,595,568 bytes)
--   λ> length (show (menorDivisible6 1 (10^5)))
--   43452
--   (0.30 secs, 172,986,840 bytes)
--   λ> length (show (menorDivisible7 1 (10^5)))
--   43452
--   (0.30 secs, 172,986,920 bytes)
--   
--   λ> length (show (menorDivisible5 1 (2*10^5)))
--   86871
--   (2.47 secs, 7,989,147,304 bytes)
--   λ> length (show (menorDivisible6 1 (2*10^5)))
--   86871
--   (0.89 secs, 533,876,496 bytes)
--   λ> length (show (menorDivisible7 1 (2*10^5)))
--   86871
--   (0.88 secs, 533,875,608 bytes)

import Data.List (foldl')

-- 1ª solución

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

menorDivisible :: Integer -> Integer -> Integer

menorDivisible a b =

head [x | x <- [b..]

, and [x `mod` y == 0 | y <- [a..b]]]

-- 2ª solución

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

menorDivisible2 :: Integer -> Integer -> Integer

menorDivisible2 a b

| a == b = a

| otherwise = lcm a (menorDivisible (a+1) b)

-- 3ª solución

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

menorDivisible3 :: Integer -> Integer -> Integer

menorDivisible3 a b = foldl lcm 1 [a..b]

-- 4ª solución

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

menorDivisible4 :: Integer -> Integer -> Integer

menorDivisible4 a b = foldl1 lcm [a..b]

-- 5ª solución

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

menorDivisible5 :: Integer -> Integer -> Integer

menorDivisible5 a b = foldl' lcm 1 [a..b]

-- 6ª solución

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

menorDivisible6 :: Integer -> Integer -> Integer

menorDivisible6 a b = foldr1 lcm [a..b]

-- 7ª solución

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

menorDivisible7 :: Integer -> Integer -> Integer

menorDivisible7 a = foldr1 lcm . enumFromTo a

-- Comparación de eficiencia

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

-- λ> menorDivisible 1 17

-- 12252240

-- (18.63 secs, 15,789,475,488 bytes)

-- λ> menorDivisible2 1 17

-- 12252240

-- (13.29 secs, 11,868,764,272 bytes)

-- λ> menorDivisible3 1 17

-- 12252240

-- (0.00 secs, 114,688 bytes)

-- λ> menorDivisible4 1 17

-- 12252240

-- (0.01 secs, 114,752 bytes)

-- λ> menorDivisible5 1 17

-- 12252240

-- (0.01 secs, 110,640 bytes)

-- λ> menorDivisible6 1 17

-- 12252240

-- (0.01 secs, 114,752 bytes)

-- λ> menorDivisible7 1 17

-- 12252240

-- (0.00 secs, 110,912 bytes)

-- λ> length (show (menorDivisible3 1 (10^5)))

-- 43452

-- (1.54 secs, 2,021,887,000 bytes)

-- λ> length (show (menorDivisible4 1 (10^5)))

-- 43452

-- (1.47 secs, 2,021,886,616 bytes)

-- λ> length (show (menorDivisible5 1 (10^5)))

-- 43452

-- (0.65 secs, 2,009,595,568 bytes)

-- λ> length (show (menorDivisible6 1 (10^5)))

-- 43452

-- (0.30 secs, 172,986,840 bytes)

-- λ> length (show (menorDivisible7 1 (10^5)))

-- 43452

-- (0.30 secs, 172,986,920 bytes)

-- λ> length (show (menorDivisible5 1 (2*10^5)))

-- 86871

-- (2.47 secs, 7,989,147,304 bytes)

-- λ> length (show (menorDivisible6 1 (2*10^5)))

-- 86871

-- (0.89 secs, 533,876,496 bytes)

-- λ> length (show (menorDivisible7 1 (2*10^5)))

-- 86871

-- (0.88 secs, 533,875,608 bytes)

Pensamiento

Será el peor de los malos
bribón que olvide
su vocación de diablo.

Antonio Machado

Mayor divisor primo

PorJosé A. Alonso 26-noviembre-20193-diciembre-2019

Los divisores primos de 13195 son 5, 7, 13 y 29. Por tanto, el mayor divisor primo de 13195 es 29.

Definir la función

   mayorDivisorPrimo :: Integer -> Integer

1	mayorDivisorPrimo :: Integer -> Integer

tal que (mayorDivisorPrimo n) es el mayor divisor primo de n. Por ejemplo,

   mayorDivisorPrimo 13195            ==  29
   mayorDivisorPrimo 152416333181401  ==  12345701

1 2	mayorDivisorPrimo 13195 == 29 mayorDivisorPrimo 152416333181401 == 12345701

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

Soluciones

import Data.Numbers.Primes (primeFactors)

-- 1ª solución  (sin librerías auxiliares)
-- =======================================

mayorDivisorPrimo :: Integer -> Integer
mayorDivisorPrimo = last . divisoresPrimos 

-- (divisoresPrimos n) es la lista de los divisores primos de n. Por
-- ejemplo, 
--    divisoresPrimos 13195  ==  [5,7,13,29]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos 0 = []
divisoresPrimos 1 = []
divisoresPrimos n = m : divisoresPrimos (n `div` m)
  where m = menorDivisorPrimo n 

-- (menorDivisorPrimo n) es el menor divisor primo de n. Por ejemplo, 
--    menorDivisorPrimo 24  ==  2
--    menorDivisorPrimo 25  ==  5
--    menorDivisorPrimo 29  ==  29
menorDivisorPrimo :: Integer -> Integer
menorDivisorPrimo x =
  head [y | y <- 2 : [3,5..(ceiling . sqrt . fromIntegral) x] ++ [x]
          , x `mod` y == 0]

-- 2ª solución (con la librería Data.Numbers.Primes)
-- =================================================

mayorDivisorPrimo2 :: Integer -> Integer
mayorDivisorPrimo2 = last . primeFactors

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

--   λ> mayorDivisorPrimo 152416333181401
--   12345701
--   (1.96 secs, 1,630,201,856 bytes)
--   λ> mayorDivisorPrimo2 152416333181401
--   12345701
--   (2.01 secs, 5,445,284,432 bytes)

import Data.Numbers.Primes (primeFactors)

-- 1ª solución (sin librerías auxiliares)

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

mayorDivisorPrimo :: Integer -> Integer

mayorDivisorPrimo = last . divisoresPrimos

-- (divisoresPrimos n) es la lista de los divisores primos de n. Por

-- ejemplo,

-- divisoresPrimos 13195 == [5,7,13,29]

divisoresPrimos :: Integer -> [Integer]

divisoresPrimos 0 = []

divisoresPrimos 1 = []

divisoresPrimos n = m : divisoresPrimos (n `div` m)

where m = menorDivisorPrimo n

-- (menorDivisorPrimo n) es el menor divisor primo de n. Por ejemplo,

-- menorDivisorPrimo 24 == 2

-- menorDivisorPrimo 25 == 5

-- menorDivisorPrimo 29 == 29

menorDivisorPrimo :: Integer -> Integer

menorDivisorPrimo x =

head [y | y <- 2 : [3,5..(ceiling . sqrt . fromIntegral) x] ++ [x]

, x `mod` y == 0]

-- 2ª solución (con la librería Data.Numbers.Primes)

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

mayorDivisorPrimo2 :: Integer -> Integer

mayorDivisorPrimo2 = last . primeFactors

-- Comparación de eficiencia

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

-- λ> mayorDivisorPrimo 152416333181401

-- 12345701

-- (1.96 secs, 1,630,201,856 bytes)

-- λ> mayorDivisorPrimo2 152416333181401

-- 12345701

-- (2.01 secs, 5,445,284,432 bytes)

Pensamiento

«Un programa de ordenador es una demostración.» ~ Igor Rivin

Factorización prima

PorJosé A. Alonso 21-noviembre-201926-marzo-2022

La descomposición prima de 600 es

   600 = 2³ * 3 * 5²

1	600 = 2³ * 3 * 5²

Definir la función

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

1	factorizacion :: Integer -> [(Integer,Integer)]

tal que (factorizacion x) ses la lista de las bases y exponentes de la descomposición prima de x. Por ejemplo,

   factorizacion 600  ==  [(2,3),(3,1),(5,2)]
   length (factorizacion (product [1..3*10^4]))  ==  3245

1 2	factorizacion 600 == [(2,3),(3,1),(5,2)] length (factorizacion (product [1..3*10^4])) == 3245

Soluciones

import Data.List (genericLength, group, inits, nub, sort, subsequences)
import Data.Numbers.Primes (primes, primeFactors)

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

factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(x,nOcurrencias x xs) | x <- elementos xs]
  where xs = factoresPrimos n

-- (factores primos n) es la lista de los factores primos de n. Por
-- ejemplo, 
--   factoresPrimos 600  ==  [2,2,2,3,5,5]
factoresPrimos :: Integer -> [Integer]
factoresPrimos 1 = []
factoresPrimos n = x : factoresPrimos (n `div` x)
  where x = menorFactor n

-- (menorFactor n) es el menor factor primo de n. Por ejemplo,
--   menorFactor 10  ==  2
--   menorFactor 11  ==  11
menorFactor :: Integer -> Integer
menorFactor n = head [x | x <- [2..n], n `mod` x == 0]

-- (elementos xs) es la lista de los elementos, sin repeticiones, de
-- xs. Por ejemplo,
--   elementos [3,2,3,5,2]  ==  [3,2,5]
elementos :: Eq a => [a] -> [a]
elementos [] = []
elementos (x:xs) = x : elementos (filter (/=x) xs)

-- (nOcurrencias x ys) es el número de ocurrencias de x en ys. Por
-- ejemplo, 
--   nOcurrencias 'a' "Salamanca"  ==  4
nOcurrencias :: Eq a => a -> [a] -> Integer
nOcurrencias x ys = genericLength (filter (==x) ys) 

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

factorizacion2 :: Integer -> [(Integer,Integer)]
factorizacion2 n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

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

factorizacion3 :: Integer -> [(Integer,Integer)]
factorizacion3 = 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)

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

--   λ> length (factorizacion (product [1..10^4]))
--   1229
--   (4.84 secs, 2,583,331,768 bytes)
--   λ> length (factorizacion2 (product [1..10^4]))
--   1229
--   (0.24 secs, 452,543,360 bytes)
--   λ> length (factorizacion3 (product [1..10^4]))
--   1229
--   (0.23 secs, 452,433,504 bytes)
--   
--   λ> length (factorizacion (product (take (2*10^3) primes)))
--   2000
--   (6.58 secs, 3,415,098,552 bytes)
--   λ> length (factorizacion2 (product (take (2*10^3) primes)))
--   2000
--   (0.02 secs, 23,060,512 bytes)
--   λ> length (factorizacion3 (product (take (2*10^3) primes)))
--   2000
--   (0.02 secs, 22,882,080 bytes)

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

import Data.Numbers.Primes (primes, primeFactors)

-- 1ª solución

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

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

factorizacion n =

[(x,nOcurrencias x xs) | x <- elementos xs]

where xs = factoresPrimos n

-- (factores primos n) es la lista de los factores primos de n. Por

-- ejemplo,

-- factoresPrimos 600 == [2,2,2,3,5,5]

factoresPrimos :: Integer -> [Integer]

factoresPrimos 1 = []

factoresPrimos n = x : factoresPrimos (n `div` x)

where x = menorFactor n

-- (menorFactor n) es el menor factor primo de n. Por ejemplo,

-- menorFactor 10 == 2

-- menorFactor 11 == 11

menorFactor :: Integer -> Integer

menorFactor n = head [x | x <- [2..n], n `mod` x == 0]

-- (elementos xs) es la lista de los elementos, sin repeticiones, de

-- xs. Por ejemplo,

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

elementos :: Eq a => [a] -> [a]

elementos [] = []

elementos (x:xs) = x : elementos (filter (/=x) xs)

-- (nOcurrencias x ys) es el número de ocurrencias de x en ys. Por

-- ejemplo,

-- nOcurrencias 'a' "Salamanca" == 4

nOcurrencias :: Eq a => a -> [a] -> Integer

nOcurrencias x ys = genericLength (filter (==x) ys)

-- 2ª solución

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

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

factorizacion2 n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 3ª solución

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

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

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

-- Comparación de eficiencia de sumaDivisores

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

-- λ> length (factorizacion (product [1..10^4]))

-- 1229

-- (4.84 secs, 2,583,331,768 bytes)

-- λ> length (factorizacion2 (product [1..10^4]))

-- 1229

-- (0.24 secs, 452,543,360 bytes)

-- λ> length (factorizacion3 (product [1..10^4]))

-- 1229

-- (0.23 secs, 452,433,504 bytes)

-- λ> length (factorizacion (product (take (2*10^3) primes)))

-- 2000

-- (6.58 secs, 3,415,098,552 bytes)

-- λ> length (factorizacion2 (product (take (2*10^3) primes)))

-- 2000

-- (0.02 secs, 23,060,512 bytes)

-- λ> length (factorizacion3 (product (take (2*10^3) primes)))

-- 2000

-- (0.02 secs, 22,882,080 bytes)

Pensamiento

¿Todo para los demás?
Mancebo, llena tu jarro,
que ya te lo beberán.

Antonio Machado

Período de una lista

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

Definir la función

   periodo :: Eq a => [a] -> [a]

1	periodo :: Eq a => [a] -> [a]

tal que (periodo xs) es el período de xs. Por ejemplo,

   periodo "ababab"      ==  "ab"
   periodo "buenobueno"  ==  "bueno"
   periodo "oooooo"      ==  "o"
   periodo "sevilla"     ==  "sevilla"

periodo "ababab" == "ab"

periodo "buenobueno" == "bueno"

periodo "oooooo" == "o"

periodo "sevilla" == "sevilla"

Soluciones

import Data.List (isPrefixOf, inits)

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

periodo1 :: Eq a => [a] -> [a]
periodo1 xs = take n xs
    where l = length xs
          n = head [m | m <- divisores l, 
                        concat (replicate (l `div` m) (take m xs)) == xs]

-- (divisores n) es la lista de los divisores de n. Por ejemplo,
--    divisores 96  ==  [1,2,3,4,6,8,12,16,24,32,48,96]
divisores :: Int -> [Int]
divisores n = [x | x <- [1..n], n `mod` x == 0]

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

periodo2 :: Eq a => [a] -> [a]
periodo2 xs = take n xs
    where l = length xs
          n = head [m | m <- divisores l, 
                        xs `isPrefixOf` cycle (take m xs)]

import Data.List (isPrefixOf, inits)

-- 1ª solución

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

periodo1 :: Eq a => [a] -> [a]

periodo1 xs = take n xs

where l = length xs

n = head [m | m <- divisores l,

concat (replicate (l `div` m) (take m xs)) == xs]

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

-- divisores 96 == [1,2,3,4,6,8,12,16,24,32,48,96]

divisores :: Int -> [Int]

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

-- 2ª solución

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

periodo2 :: Eq a => [a] -> [a]

periodo2 xs = take n xs

where l = length xs

n = head [m | m <- divisores l,

xs `isPrefixOf` cycle (take m xs)]

Pensamiento

Esta luz de Sevilla … Es el palacio
donde nací, con su rumor de fuente.

Antonio Machado

Mayor capicúa producto de dos números de n cifras

PorJosé A. Alonso 25-abril-201925-abril-2021

Un capicúa es un número que es igual leído de izquierda a derecha que de derecha a izquierda.

Definir la función

   mayorCapicuaP :: Integer -> Integer

1	mayorCapicuaP :: Integer -> Integer

tal que (mayorCapicuaP n) es el mayor capicúa que es el producto de dos números de n cifras. Por ejemplo,

   mayorCapicuaP 2  ==  9009
   mayorCapicuaP 3  ==  906609
   mayorCapicuaP 4  ==  99000099
   mayorCapicuaP 5  ==  9966006699

mayorCapicuaP 2 == 9009

mayorCapicuaP 3 == 906609

mayorCapicuaP 4 == 99000099

mayorCapicuaP 5 == 9966006699

Soluciones

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

mayorCapicuaP1 :: Integer -> Integer
mayorCapicuaP1 n = head (capicuasP n)

-- (capicuasP n) es la lista de las capicúas de 2*n cifras que
-- pueden escribirse como productos de dos números de n cifras. Por
-- ejemplo, Por ejemplo,
--    ghci> capicuasP 2
--    [9009,8448,8118,8008,7227,7007,6776,6336,6006,5775,5445,5335,
--     5225,5115,5005,4884,4774,4664,4554,4224,4004,3773,3663,3003,
--     2992,2772,2552,2442,2332,2112,2002,1881,1771,1551,1221,1001]
capicuasP n = [x | x <- capicuas n,
                        not (null (productosDosNumerosCifras n x))]

-- (capicuas n) es la lista de las capicúas de 2*n cifras de mayor a
-- menor. Por ejemplo, 
--    capicuas 1           ==  [99,88,77,66,55,44,33,22,11]
--    take 7 (capicuas 2)  ==  [9999,9889,9779,9669,9559,9449,9339]
capicuas :: Integer -> [Integer]
capicuas n = [capicua x | x <- numerosCifras n]

-- (numerosCifras n) es la lista de los números de n cifras de mayor a
-- menor. Por ejemplo,
--    numerosCifras 1           ==  [9,8,7,6,5,4,3,2,1]
--    take 7 (numerosCifras 2)  ==  [99,98,97,96,95,94,93]
--    take 7 (numerosCifras 3)  ==  [999,998,997,996,995,994,993]
numerosCifras :: Integer -> [Integer]
numerosCifras n = [a,a-1..b]
  where a = 10^n-1
        b = 10^(n-1) 

-- (capicua n) es la capicúa formada añadiendo el inverso de n a
--  continuación de n. Por ejemplo,
--    capicua 93  ==  9339
capicua :: Integer -> Integer
capicua n = read (xs ++ (reverse xs))
  where xs = show n

-- (productosDosNumerosCifras n x) es la lista de los números y de n
-- cifras tales que existe un z de n cifras y x es el producto de y por
-- z. Por ejemplo, 
--    productosDosNumerosCifras 2 9009  ==  [99,91]
productosDosNumerosCifras n x = [y | y <- numeros,
                                     mod x y == 0,
                                     div x y `elem` numeros]
  where numeros = numerosCifras n

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

mayorCapicuaP2 :: Integer -> Integer
mayorCapicuaP2 n = maximum [x*y | x <- [a,a-1..b],
                                  y <- [a,a-1..b],
                                  esCapicua (x*y)] 
  where a = 10^n-1
        b = 10^(n-1)

-- (esCapicua x) se verifica si x es capicúa. Por ejemplo,
--    esCapicua 353  ==  True
--    esCapicua 357  ==  False
esCapicua :: Integer -> Bool
esCapicua n = xs == reverse xs
  where xs = show n

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

mayorCapicuaP3 :: Integer -> Integer
mayorCapicuaP3 n = maximum [x*y | (x,y) <- pares a b, 
                                  esCapicua (x*y)] 
  where a = 10^n-1
        b = 10^(n-1)

-- (pares a b) es la lista de los pares de números entre a y b de forma
-- que su suma es decreciente. Por ejemplo,
--    pares 9 7  ==  [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]
pares a b = [(x,z-x) | z <- [a1,a1-1..b1],
                       x <- [a,a-1..b],
                       x <= z-x, z-x <= a]
  where a1 = 2*a
        b1 = 2*b

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

mayorCapicuaP4 :: Integer -> Integer
mayorCapicuaP4 n = maximum [x | y <- [a..b],
                                z <- [y..b],
                                let x = y * z,
                                let s = show x,
                                s == reverse s]
  where a = 10^(n-1)
        b = 10^n-1
     
-- 5ª solución
-- ===========

mayorCapicuaP5 :: Integer -> Integer
mayorCapicuaP5 n = maximum [x*y | (x,y) <- pares2 b a, esCapicua (x*y)]
  where a = 10^(n-1)
        b = 10^n-1
                       
-- (pares2 a b) es la lista de los pares de números entre a y b de forma
-- que su suma es decreciente. Por ejemplo,
--    pares2 9 7  ==  [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]
pares2 a b = [(x,y) | x <- [a,a-1..b], y <- [a,a-1..x]]

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

mayorCapicuaP6 :: Integer -> Integer
mayorCapicuaP6 n = maximum [x*y | x <- [a..b], 
                                  y <- [x..b] , 
                                  esCapicua (x*y)]
  where a = 10^(n-1)
        b = 10^n-1
 
-- (cifras n) es la lista de las cifras de n en orden inverso. Por
-- ejemplo,  
--    cifras 325  == [5,2,3]
cifras :: Integer -> [Integer]
cifras n 
    | n < 10    = [n]
    | otherwise = (ultima n) : (cifras (quitarUltima n))

-- (ultima n) es la última cifra de n. Por ejemplo,
--    ultima 325  ==  5
ultima  :: Integer -> Integer
ultima n =  n - (n `div` 10)*10

-- (quitarUltima n) es el número obtenido al quitarle a n su última
-- cifra. Por ejemplo,
--    quitarUltima 325  =>  32 
quitarUltima :: Integer -> Integer
quitarUltima n = (n - (ultima n)) `div` 10

-- 7ª solución
-- ===========

mayorCapicuaP7 :: Integer -> Integer
mayorCapicuaP7 n = head [x | x <- capicuas n, esFactorizable x n]

-- (esFactorizable x n) se verifica si x se puede escribir como producto
-- de dos números de n dígitos. Por ejemplo,
--    esFactorizable 1219 2  ==  True
--    esFactorizable 1217 2  ==  False
esFactorizable x n = aux i x
  where b = 10^n-1
        i = floor (sqrt (fromIntegral x))
        aux i x | i > b          = False
                | x `mod` i == 0 = x `div` i < b 
                | otherwise      = aux (i+1) x

-- ---------------------------------------------------------------------
-- Comparación de soluciones                                          --
-- ---------------------------------------------------------------------

-- El tiempo de cálculo de (mayorCapicuaP n) para las 7 definiciones es
--    +------+------+------+------+
--    | Def. | 2    | 3    | 4    |
--    |------+------+------+------|
--    |    1 | 0.01 | 0.13 | 1.39 |
--    |    2 | 0.03 | 2.07 |      |
--    |    3 | 0.05 | 3.86 |      |
--    |    4 | 0.01 | 0.89 |      |
--    |    5 | 0.03 | 1.23 |      |
--    |    6 | 0.02 | 1.03 |      |
--    |    7 | 0.01 | 0.02 | 0.02 |
--    +------+------+------+------+

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

-- 1ª solución

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

mayorCapicuaP1 :: Integer -> Integer

mayorCapicuaP1 n = head (capicuasP n)

-- (capicuasP n) es la lista de las capicúas de 2*n cifras que

-- pueden escribirse como productos de dos números de n cifras. Por

-- ejemplo, Por ejemplo,

-- ghci> capicuasP 2

-- [9009,8448,8118,8008,7227,7007,6776,6336,6006,5775,5445,5335,

-- 5225,5115,5005,4884,4774,4664,4554,4224,4004,3773,3663,3003,

-- 2992,2772,2552,2442,2332,2112,2002,1881,1771,1551,1221,1001]

capicuasP n = [x | x <- capicuas n,

not (null (productosDosNumerosCifras n x))]

-- (capicuas n) es la lista de las capicúas de 2*n cifras de mayor a

-- menor. Por ejemplo,

-- capicuas 1 == [99,88,77,66,55,44,33,22,11]

-- take 7 (capicuas 2) == [9999,9889,9779,9669,9559,9449,9339]

capicuas :: Integer -> [Integer]

capicuas n = [capicua x | x <- numerosCifras n]

-- (numerosCifras n) es la lista de los números de n cifras de mayor a

-- menor. Por ejemplo,

-- numerosCifras 1 == [9,8,7,6,5,4,3,2,1]

-- take 7 (numerosCifras 2) == [99,98,97,96,95,94,93]

-- take 7 (numerosCifras 3) == [999,998,997,996,995,994,993]

numerosCifras :: Integer -> [Integer]

numerosCifras n = [a,a-1..b]

where a = 10^n-1

b = 10^(n-1)

-- (capicua n) es la capicúa formada añadiendo el inverso de n a

-- continuación de n. Por ejemplo,

-- capicua 93 == 9339

capicua :: Integer -> Integer

capicua n = read (xs ++ (reverse xs))

where xs = show n

-- (productosDosNumerosCifras n x) es la lista de los números y de n

-- cifras tales que existe un z de n cifras y x es el producto de y por

-- z. Por ejemplo,

-- productosDosNumerosCifras 2 9009 == [99,91]

productosDosNumerosCifras n x = [y | y <- numeros,

mod x y == 0,

div x y `elem` numeros]

where numeros = numerosCifras n

-- 2ª solución

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

mayorCapicuaP2 :: Integer -> Integer

mayorCapicuaP2 n = maximum [x*y | x <- [a,a-1..b],

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

esCapicua (x*y)]

where a = 10^n-1

b = 10^(n-1)

-- (esCapicua x) se verifica si x es capicúa. Por ejemplo,

-- esCapicua 353 == True

-- esCapicua 357 == False

esCapicua :: Integer -> Bool

esCapicua n = xs == reverse xs

where xs = show n

-- 3ª solución

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

mayorCapicuaP3 :: Integer -> Integer

mayorCapicuaP3 n = maximum [x*y | (x,y) <- pares a b,

esCapicua (x*y)]

where a = 10^n-1

b = 10^(n-1)

-- (pares a b) es la lista de los pares de números entre a y b de forma

-- que su suma es decreciente. Por ejemplo,

-- pares 9 7 == [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]

pares a b = [(x,z-x) | z <- [a1,a1-1..b1],

x <- [a,a-1..b],

x <= z-x, z-x <= a]

where a1 = 2*a

b1 = 2*b

-- 4ª solución

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

mayorCapicuaP4 :: Integer -> Integer

mayorCapicuaP4 n = maximum [x | y <- [a..b],

z <- [y..b],

let x = y * z,

let s = show x,

s == reverse s]

where a = 10^(n-1)

b = 10^n-1

-- 5ª solución

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

mayorCapicuaP5 :: Integer -> Integer

mayorCapicuaP5 n = maximum [x*y | (x,y) <- pares2 b a, esCapicua (x*y)]

where a = 10^(n-1)

b = 10^n-1

-- (pares2 a b) es la lista de los pares de números entre a y b de forma

-- que su suma es decreciente. Por ejemplo,

-- pares2 9 7 == [(9,9),(8,9),(8,8),(7,9),(7,8),(7,7)]

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

-- 6ª solución

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

mayorCapicuaP6 :: Integer -> Integer

mayorCapicuaP6 n = maximum [x*y | x <- [a..b],

y <- [x..b] ,

esCapicua (x*y)]

where a = 10^(n-1)

b = 10^n-1

-- (cifras n) es la lista de las cifras de n en orden inverso. Por

-- ejemplo,

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

cifras :: Integer -> [Integer]

cifras n

| n < 10 = [n]

| otherwise = (ultima n) : (cifras (quitarUltima n))

-- (ultima n) es la última cifra de n. Por ejemplo,

-- ultima 325 == 5

ultima :: Integer -> Integer

ultima n = n - (n `div` 10)*10

-- (quitarUltima n) es el número obtenido al quitarle a n su última

-- cifra. Por ejemplo,

-- quitarUltima 325 => 32

quitarUltima :: Integer -> Integer

quitarUltima n = (n - (ultima n)) `div` 10

-- 7ª solución

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

mayorCapicuaP7 :: Integer -> Integer

mayorCapicuaP7 n = head [x | x <- capicuas n, esFactorizable x n]

-- (esFactorizable x n) se verifica si x se puede escribir como producto

-- de dos números de n dígitos. Por ejemplo,

-- esFactorizable 1219 2 == True

-- esFactorizable 1217 2 == False

esFactorizable x n = aux i x

where b = 10^n-1

i = floor (sqrt (fromIntegral x))

aux i x | i > b = False

| x `mod` i == 0 = x `div` i < b

| otherwise = aux (i+1) x

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

-- Comparación de soluciones --

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

-- El tiempo de cálculo de (mayorCapicuaP n) para las 7 definiciones es

-- +------+------+------+------+

-- | Def. | 2 | 3 | 4 |

-- |------+------+------+------|

-- | 1 | 0.01 | 0.13 | 1.39 |

-- | 2 | 0.03 | 2.07 | |

-- | 3 | 0.05 | 3.86 | |

-- | 4 | 0.01 | 0.89 | |

-- | 5 | 0.03 | 1.23 | |

-- | 6 | 0.02 | 1.03 | |

-- | 7 | 0.01 | 0.02 | 0.02 |

-- +------+------+------+------+

Pensamiento

Mi corazón está donde ha nacido,
no a la vida, al amor, cerca del Duero.

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

Soluciones

Soluciones

Soluciones

Soluciones

Nuevas soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Pensamiento

Otras soluciones

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Referencia

Pensamiento

Soluciones

Pensamiento

Soluciones

Referencia

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

Entradas recientes

Buscar

Correo electrónico