Teoría de números – Exercitium

Números de Pentanacci

PorJosé A. Alonso 5-agosto-202214-diciembre-2022

Los números de Fibonacci se definen mediante las ecuaciones

   F(0) = 0
   F(1) = 1
   F(n) = F(n-1) + F(n-2), si n > 1

F(0) = 0

F(1) = 1

F(n) = F(n-1) + F(n-2), si n > 1

Los primeros números de Fibonacci son

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

Una generalización de los anteriores son los números de Pentanacci definidos por las siguientes ecuaciones

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

P(0) = 0

P(1) = 1

P(2) = 1

P(3) = 2

P(4) = 4

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

Los primeros números de Pentanacci son

  0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525, ...

1	0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525, ...

Definir la sucesión

   pentanacci :: [Integer]

1	pentanacci :: [Integer]

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

   λ> take 15 pentanacci
   [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]
   λ> (pentanacci !! (10^5)) `mod` (10^30) 
   482929150584077921552549215816
   231437922897686901289110700696
   λ> length (show (pentanacci !! (10^5)))
   29357

λ> take 15 pentanacci

[0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]

λ> (pentanacci !! (10^5)) `mod` (10^30)

482929150584077921552549215816

231437922897686901289110700696

λ> length (show (pentanacci !! (10^5)))

29357

Soluciones

import Data.List (zipWith5)
import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

pentanacci1 :: [Integer]
pentanacci1 = [pent n | n <- [0..]]

pent :: Integer -> Integer
pent 0 = 0
pent 1 = 1
pent 2 = 1
pent 3 = 2
pent 4 = 4
pent n = pent (n-1) + pent (n-2) + pent (n-3) + pent (n-4) + pent (n-5)

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

pentanacci2 :: [Integer]
pentanacci2 = 
  0 : 1 : 1 : 2 : 4 : zipWith5 f (r 0) (r 1) (r 2) (r 3) (r 4)
  where f a b c d e = a+b+c+d+e
        r n         = drop n pentanacci2

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

pentanacci3 :: [Integer]
pentanacci3 = p (0, 1, 1, 2, 4)
  where p (a, b, c, d, e) = a : p (b, c, d, e, a + b + c + d + e)

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

pentanacci4 :: [Integer]
pentanacci4 = 0: 1: 1: 2: 4: p pentanacci4
  where p (a:b:c:d:e:xs) = (a+b+c+d+e): p (b:c:d:e:xs)

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

-- La propiedad es
prop_pentanacci :: NonNegative Int -> Bool
prop_pentanacci (NonNegative n) =
  all (== pentanacci1 !! n)
      [pentanacci1 !! n,
       pentanacci2 !! n,
       pentanacci3 !! n]

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

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

-- La comparación es
--    λ> pentanacci1 !! 25
--    5976577
--    (3.18 secs, 1,025,263,896 bytes)
--    λ> pentanacci2 !! 25
--    5976577
--    (0.00 secs, 562,360 bytes)
--    
--    λ> length (show (pentanacci2 !! (10^5)))
--    29357
--    (1.04 secs, 2,531,259,408 bytes)
--    λ> length (show (pentanacci3 !! (10^5)))
--    29357
--    (1.00 secs, 2,548,868,384 bytes)
--    λ> length (show (pentanacci4 !! (10^5)))
--    29357
--    (0.96 secs, 2,580,065,520 bytes)

import Data.List (zipWith5)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

pentanacci1 :: [Integer]

pentanacci1 = [pent n | n <- [0..]]

pent :: Integer -> Integer

pent 0 = 0

pent 1 = 1

pent 2 = 1

pent 3 = 2

pent 4 = 4

pent n = pent (n-1) + pent (n-2) + pent (n-3) + pent (n-4) + pent (n-5)

-- 2ª solución

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

pentanacci2 :: [Integer]

pentanacci2 =

0 : 1 : 1 : 2 : 4 : zipWith5 f (r 0) (r 1) (r 2) (r 3) (r 4)

where f a b c d e = a+b+c+d+e

r n = drop n pentanacci2

-- 3ª solución

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

pentanacci3 :: [Integer]

pentanacci3 = p (0, 1, 1, 2, 4)

where p (a, b, c, d, e) = a : p (b, c, d, e, a + b + c + d + e)

-- 4ª solución

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

pentanacci4 :: [Integer]

pentanacci4 = 0: 1: 1: 2: 4: p pentanacci4

where p (a:b:c:d:e:xs) = (a+b+c+d+e): p (b:c:d:e:xs)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_pentanacci :: NonNegative Int -> Bool

prop_pentanacci (NonNegative n) =

all (== pentanacci1 !! n)

[pentanacci1 !! n,

pentanacci2 !! n,

pentanacci3 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> pentanacci1 !! 25

-- 5976577

-- (3.18 secs, 1,025,263,896 bytes)

-- λ> pentanacci2 !! 25

-- 5976577

-- (0.00 secs, 562,360 bytes)

-- λ> length (show (pentanacci2 !! (10^5)))

-- 29357

-- (1.04 secs, 2,531,259,408 bytes)

-- λ> length (show (pentanacci3 !! (10^5)))

-- 29357

-- (1.00 secs, 2,548,868,384 bytes)

-- λ> length (show (pentanacci4 !! (10^5)))

-- 29357

-- (0.96 secs, 2,580,065,520 bytes)

El código se encuentra en GitHub.

Referencias

Tito III Piezas y Eric Weisstein, Pentanacci number.
N. J. A. Sloane, Sucesión A001591 de la OEIS.

Número de representaciones de n como suma de dos cuadrados

PorJosé A. Alonso 4-agosto-202214-diciembre-2022

Sea $n$ un número natural cuya factorización prima es
$$n = 2^{a} \times p(1)^{b(1)} \times \dots \times p(n)^{b(n)} \times q(1)^{c(1)} \times \dots \times q(m)^{c(m)}$$
donde los $p(i)$ son los divisores primos de $n$ congruentes con 3 módulo 4 y los $q(j)$ son los divisores primos de $n$ congruentes con 1 módulo 4. Entonces, el número de forma de descomponer $n$ como suma de dos
cuadrados es 0, si algún $b(i)$ es impar y es el techo (es decir, el número entero más próximo por exceso) de
$$\frac{(1+c(1)) \times \dots \times (1+c(m))}{2}$$
en caso contrario. Por ejemplo, el número
$$2^{3} \times (3^{9} \times 7^{8}) \times (5^{3} \times 13^{6})$$
no se puede descomponer como sumas de dos cuadrados (porque el exponente de 3 es impar) y el número
$$2^{3} \times (3^{2} \times 7^{8}) \times (5^{3} \times 13^{6})$$
tiene 14 descomposiciones como suma de dos cuadrados (porque los exponentes de 3 y 7 son pares y el techo de
$$\frac{(1+3) \times (1+6)}{2}$$
es 14).

Definir la función

   nRepresentaciones :: Integer -> Integer

1	nRepresentaciones :: Integer -> Integer

tal que (nRepresentaciones n) es el número de formas de representar n como suma de dos cuadrados. Por ejemplo,

   nRepresentaciones (2^3*3^9*5^3*7^8*13^6)        ==  0
   nRepresentaciones (2^3*3^2*5^3*7^8*13^6)        ==  14
   head [n | n <- [1..], nRepresentaciones n > 8]  ==  71825

nRepresentaciones (2^3*3^9*5^3*7^8*13^6) == 0

nRepresentaciones (2^3*3^2*5^3*7^8*13^6) == 14

head [n | n <- [1..], nRepresentaciones n > 8] == 71825

Usando la función representaciones del ejercicio anterior, comprobar con QuickCheck la siguiente propiedad

   prop_representacion :: Positive Integer -> Bool
   prop_representacion (Positive n) =
     nRepresentaciones2 n == genericLength (representaciones n)

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

nRepresentaciones2 n == genericLength (representaciones n)

Soluciones

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

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

nRepresentaciones1 :: Integer -> Integer
nRepresentaciones1 n =
  ceiling (fromIntegral (product (map aux (factorizacion n))) / 2)
  where aux (p,e) | p == 2         = 1
                  | p `mod` 4 == 3 = if even e then 1 else 0
                  | otherwise      = e+1

-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

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

nRepresentaciones2 :: Integer -> Integer
nRepresentaciones2 n =
  (1 + product (map aux (factorizacion n))) `div`  2
  where aux (p,e) | p == 2         = 1
                  | p `mod` 4 == 3 = if even e then 1 else 0
                  | otherwise      = e+1

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

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

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

-- La comparación es
--    λ> head [n | n <- [1..], nRepresentaciones1 n > 8]
--    71825
--    (1.39 secs, 2,970,063,760 bytes)
--    λ> head [n | n <- [1..], nRepresentaciones2 n > 8]
--    71825
--    (1.71 secs, 2,943,788,424 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_representacion :: Positive Integer -> Bool
prop_representacion (Positive n) =
  nRepresentaciones2 n == genericLength (representaciones n)

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

esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))
    
-- La comprobación es
--    λ> quickCheck prop_representacion
--    +++ OK, passed 100 tests.

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

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

-- 1ª solución

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

nRepresentaciones1 :: Integer -> Integer

nRepresentaciones1 n =

ceiling (fromIntegral (product (map aux (factorizacion n))) / 2)

where aux (p,e) | p == 2 = 1

| p `mod` 4 == 3 = if even e then 1 else 0

| otherwise = e+1

-- (factorizacion n) es la factorización prima de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

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

factorizacion n =

map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

-- 2ª solución

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

nRepresentaciones2 :: Integer -> Integer

nRepresentaciones2 n =

(1 + product (map aux (factorizacion n))) `div` 2

where aux (p,e) | p == 2 = 1

| p `mod` 4 == 3 = if even e then 1 else 0

| otherwise = e+1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_nRepresentaciones :: Positive Integer -> Bool

prop_nRepresentaciones (Positive n) =

nRepresentaciones1 n == nRepresentaciones2 n

-- La comprobación es

-- λ> quickCheck prop_nRepresentaciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> head [n | n <- [1..], nRepresentaciones1 n > 8]

-- 71825

-- (1.39 secs, 2,970,063,760 bytes)

-- λ> head [n | n <- [1..], nRepresentaciones2 n > 8]

-- 71825

-- (1.71 secs, 2,943,788,424 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

nRepresentaciones2 n == genericLength (representaciones n)

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

representaciones n =

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

let z = n - x*x,

esCuadrado z]

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- La comprobación es

-- λ> quickCheck prop_representacion

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Referencias

W. Stein, Which numbers are the sum of two squares?.
E.W. Weisstein, Sum of squares function en MathWorld.
N.J.A. Sloane,Sucesión A004018 de OEIS.
Expressing a number as a sum of two squares.
Sum of squares.

Representaciones de un número como suma de dos cuadrados

PorJosé A. Alonso 3-agosto-202214-diciembre-2022

Definir la función

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

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

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

   representaciones  20              ==  [(2,4)]
   representaciones  25              ==  [(0,5),(3,4)]
   representaciones 325              ==  [(1,18),(6,17),(10,15)]
   length (representaciones (10^14)) == 8

representaciones 20 == [(2,4)]

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

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

length (representaciones (10^14)) == 8

Comprobar con QuickCheck que un número natural n se puede como suma de dos cuadrados si, y sólo si, en la factorización prima de n todos los exponentes de sus factores primos congruentes con 3 módulo 4 son pares.

Soluciones

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

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

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

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

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

-- (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 = floor . sqrt . fromIntegral
    
-- (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

-- 3ª solución
-- ===========
 
representaciones3 :: Integer -> [(Integer, Integer)]
representaciones3 n = aux 0 (raiz 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)

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

-- La propiedad es
prop_representaciones :: Positive Integer -> Bool
prop_representaciones (Positive n) =
  all (== representaciones1 n)
      [representaciones2 n,
       representaciones3 n]

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

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

-- La comparación es
--    λ> representaciones1 4000
--    [(20,60),(36,52)]
--    (4.95 secs, 2,434,929,624 bytes)
--    λ> representaciones2 4000
--    [(20,60),(36,52)]
--    (0.00 secs, 599,800 bytes)
--    λ> representaciones3 4000
--    [(20,60),(36,52)]
--    (0.01 secs, 591,184 bytes)
--    
--    λ> length (representaciones2 (10^14))
--    8
--    (6.64 secs, 5,600,837,088 bytes)
--    λ> length (representaciones3 (10^14))
--    8
--    (9.37 secs, 4,720,548,264 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_representacion :: Positive Integer -> Bool
prop_representacion (Positive n) =
  not (null (representaciones2 n)) == 
  all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)

-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

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

100

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

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

-- 1ª solución

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

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

representaciones1 n =

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

-- 2ª solución

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

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

representaciones2 n =

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

let z = n - x*x,

esCuadrado z]

-- (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 = floor . sqrt . fromIntegral

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

-- 3ª solución

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

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

representaciones3 n = aux 0 (raiz 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)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_representaciones :: Positive Integer -> Bool

prop_representaciones (Positive n) =

all (== representaciones1 n)

[representaciones2 n,

representaciones3 n]

-- La comprobación es

-- λ> quickCheck prop_representaciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> representaciones1 4000

-- [(20,60),(36,52)]

-- (4.95 secs, 2,434,929,624 bytes)

-- λ> representaciones2 4000

-- [(20,60),(36,52)]

-- (0.00 secs, 599,800 bytes)

-- λ> representaciones3 4000

-- [(20,60),(36,52)]

-- (0.01 secs, 591,184 bytes)

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

-- 8

-- (6.64 secs, 5,600,837,088 bytes)

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

-- 8

-- (9.37 secs, 4,720,548,264 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

not (null (representaciones2 n)) ==

all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)

-- (factorizacion n) es la factorización prima de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

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

factorizacion n =

map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

-- La comprobación es

-- λ> quickCheck prop_representacion

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Factorizaciones de números de Hilbert

PorJosé A. Alonso 2-agosto-202224-julio-2022

Un número de Hilbert es un entero positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, …

Un primo de Hilbert es un número de Hilbert n que no es por ningún número de Hilbert menor que n (salvo el 1). Los primeros primos de Hilbert son 5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73, 77, 89, 93, 97, 101, 109, 113, 121, 129, 133, 137, …

Definir la función

   factorizacionesH :: Integer -> [[Integer]]

1	factorizacionesH :: Integer -> [[Integer]]

tal que (factorizacionesH n) es la listas de primos de Hilbert cuyo producto es el número de Hilbert n. Por ejemplo,

  factorizacionesH  25    ==  [[5,5]]
  factorizacionesH  45    ==  [[5,9]]
  factorizacionesH 441    ==  [[9,49],[21,21]]
  factorizacionesH 80109  ==  [[9,9,989],[9,69,129]]

factorizacionesH 25 == [[5,5]]

factorizacionesH 45 == [[5,9]]

factorizacionesH 441 == [[9,49],[21,21]]

factorizacionesH 80109 == [[9,9,989],[9,69,129]]

Comprobar con QuickCheck que todos los números de Hilbert son factorizables como producto de primos de Hilbert (aunque la factorización, como para el 441, puede no ser única).

Soluciones

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

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

factorizacionesH1 :: Integer -> [[Integer]]
factorizacionesH1 = aux primosH1
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]

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

factorizacionesH2 :: Integer -> [[Integer]]
factorizacionesH2 = aux primosH2
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
  where esPrimoH n = all noDivideAn [5,9..m]
          where noDivideAn x = n `mod` x /= 0
                m            = ceiling (sqrt (fromIntegral n))

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo 
-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)
-- (ver en https://bit.ly/3zq7h4e ).

factorizacionesH3 :: Integer -> [[Integer]]
factorizacionesH3 = aux primosH3
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH3 :: [Integer]
primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]
  where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs
          where xs = primeFactors n

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

-- La propiedad es
prop_factorizacionesH :: Positive Integer -> Bool
prop_factorizacionesH (Positive n) =
  all (== factorizacionesH1 m)
      [factorizacionesH2 m,
       factorizacionesH3 m]
  where m = 1 + 4 * n
  
-- La comprobación es
--    λ> quickCheck prop_factorizacionesH
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> factorizacionesH1 80109
--    [[9,9,989],[9,69,129]]
--    (42.77 secs, 14,899,787,640 bytes)
--    λ> factorizacionesH2 80109
--    [[9,9,989],[9,69,129]]
--    (0.26 secs, 156,051,104 bytes)
--    λ> factorizacionesH3 80109
--    [[9,9,989],[9,69,129]]
--    (0.35 secs, 1,118,236,536 bytes)

-- Propiedad de factorización
-- ==========================

-- La propiedad es
prop_factorizable :: Positive Integer -> Bool
prop_factorizable (Positive n) =
  not (null (factorizacionesH1 (1 + 4 * n)))

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

100

101

102

103

104

105

106

107

108

109

110

111

import Data.Numbers.Primes (isPrime, primeFactors)

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

-- 1ª solución

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

factorizacionesH1 :: Integer -> [[Integer]]

factorizacionesH1 = aux primosH1

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH1 :: [Integer]

primosH1 = [n | n <- tail numerosH,

divisoresH n == [1,n]]

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

-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]

numerosH :: [Integer]

numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a

-- n. Por ejemplo,

-- divisoresH 117 == [1,9,13,117]

-- divisoresH 21 == [1,21]

divisoresH :: Integer -> [Integer]

divisoresH n = [x | x <- takeWhile (<=n) numerosH,

n `mod` x == 0]

-- 2ª solución

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

factorizacionesH2 :: Integer -> [[Integer]]

factorizacionesH2 = aux primosH2

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH2 :: [Integer]

primosH2 = filter esPrimoH (tail numerosH)

where esPrimoH n = all noDivideAn [5,9..m]

where noDivideAn x = n `mod` x /= 0

m = ceiling (sqrt (fromIntegral n))

-- 3ª solución

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo

-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)

-- (ver en https://bit.ly/3zq7h4e ).

factorizacionesH3 :: Integer -> [[Integer]]

factorizacionesH3 = aux primosH3

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH3 :: [Integer]

primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]

where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs

where xs = primeFactors n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_factorizacionesH :: Positive Integer -> Bool

prop_factorizacionesH (Positive n) =

all (== factorizacionesH1 m)

[factorizacionesH2 m,

factorizacionesH3 m]

where m = 1 + 4 * n

-- La comprobación es

-- λ> quickCheck prop_factorizacionesH

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> factorizacionesH1 80109

-- [[9,9,989],[9,69,129]]

-- (42.77 secs, 14,899,787,640 bytes)

-- λ> factorizacionesH2 80109

-- [[9,9,989],[9,69,129]]

-- (0.26 secs, 156,051,104 bytes)

-- λ> factorizacionesH3 80109

-- [[9,9,989],[9,69,129]]

-- (0.35 secs, 1,118,236,536 bytes)

-- Propiedad de factorización

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

-- La propiedad es

prop_factorizable :: Positive Integer -> Bool

prop_factorizable (Positive n) =

not (null (factorizacionesH1 (1 + 4 * n)))

-- La comprobación es

-- λ> quickCheck prop_factorizable

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Referencias

Basado en el artículo Failure of unique factorization (A example of the failure of the fundamental theorem of arithmetic) de R.J. Lipton en el blog Gödel’s Lost Letter and P=NP.

Otras referencias

Wikipedia, Hilbert number.
E.W. Weisstein, Hilbert number en MathWorld.
N.J.A. Sloane, Sucesión A057948 en la OEIS.
N.J.A. Sloane, Sucesión A057949 en la OEIS.

Números primos de Hilbert

PorJosé A. Alonso 1-agosto-202223-julio-2022

Un número de Hilbert es un positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, 73, 77, 81, 85, 89, 93, 97, …

Definir la sucesión

   primosH :: [Integer]

1	primosH :: [Integer]

tal que sus elementos son los primos de Hilbert. Por ejemplo,

   take 15 primosH     == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
   primosH !! (3*10^4) == 313661

1 2	take 15 primosH == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73] primosH !! (3*10^4) == 313661

Soluciones

import Data.Numbers.Primes (isPrime, primeFactors) 
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

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

primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]

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

primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
  where esPrimoH n = all noDivideAn [5,9..m]
          where noDivideAn x = n `mod` x /= 0
                m            = ceiling (sqrt (fromIntegral n))

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo 
-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)
-- (ver en https://bit.ly/3zq7h4e ).

primosH3 :: [Integer]
primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]
  where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs
          where xs = primeFactors n

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

-- La propiedad es
prop_primosH :: NonNegative Int -> Bool
prop_primosH (NonNegative n) =
  all (== primosH1 !! n)
      [primosH2 !! n,
       primosH3 !! n]

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

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

-- La comparación es
--    λ> primosH1 !! 2000
--    16957
--    (2.16 secs, 752,085,752 bytes)
--    λ> primosH2 !! 2000
--    16957
--    (0.03 secs, 19,771,008 bytes)
--    λ> primosH3 !! 2000
--    16957
--    (0.07 secs, 152,029,168 bytes)
--    
--    λ> primosH2 !! (3*10^4)
--    313661
--    (1.44 secs, 989,761,888 bytes)
--    λ> primosH3 !! (3*10^4)
--    313661
--    (2.06 secs, 6,554,068,992 bytes)

import Data.Numbers.Primes (isPrime, primeFactors)

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

-- 1ª solución

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

primosH1 :: [Integer]

primosH1 = [n | n <- tail numerosH,

divisoresH n == [1,n]]

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

-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]

numerosH :: [Integer]

numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a

-- n. Por ejemplo,

-- divisoresH 117 == [1,9,13,117]

-- divisoresH 21 == [1,21]

divisoresH :: Integer -> [Integer]

divisoresH n = [x | x <- takeWhile (<=n) numerosH,

n `mod` x == 0]

-- 2ª solución

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

primosH2 :: [Integer]

primosH2 = filter esPrimoH (tail numerosH)

where esPrimoH n = all noDivideAn [5,9..m]

where noDivideAn x = n `mod` x /= 0

m = ceiling (sqrt (fromIntegral n))

-- 3ª solución

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo

-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)

-- (ver en https://bit.ly/3zq7h4e ).

primosH3 :: [Integer]

primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]

where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs

where xs = primeFactors n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosH :: NonNegative Int -> Bool

prop_primosH (NonNegative n) =

all (== primosH1 !! n)

[primosH2 !! n,

primosH3 !! n]

-- La comprobación es

-- λ> quickCheck prop_primosH

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosH1 !! 2000

-- 16957

-- (2.16 secs, 752,085,752 bytes)

-- λ> primosH2 !! 2000

-- 16957

-- (0.03 secs, 19,771,008 bytes)

-- λ> primosH3 !! 2000

-- 16957

-- (0.07 secs, 152,029,168 bytes)

-- λ> primosH2 !! (3*10^4)

-- 313661

-- (1.44 secs, 989,761,888 bytes)

-- λ> primosH3 !! (3*10^4)

-- 313661

-- (2.06 secs, 6,554,068,992 bytes)

El código se encuentra en GitHub.

Sumas de dos primos

PorJosé A. Alonso 29-julio-202222-julio-2022

Definir la sucesión

   sumasDeDosPrimos :: [Integer]

1	sumasDeDosPrimos :: [Integer]

cuyos elementos son los números que se pueden escribir como suma de dos números primos. Por ejemplo,

   λ> take 23 sumasDeDosPrimos
   [4,5,6,7,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,28,30,31]
   λ> sumasDeDosPrimos !! (5*10^5)
   862878

λ> take 23 sumasDeDosPrimos

[4,5,6,7,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,28,30,31]

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

862878

Soluciones

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

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

sumasDeDosPrimos1 :: [Integer]
sumasDeDosPrimos1 =
  [n | n <- [1..], not (null (sumaDeDosPrimos1 n))]

-- (sumaDeDosPrimos1 n) es la lista de pares de primos cuya suma es
-- n. Por ejemplo,
--    sumaDeDosPrimos  9  ==  [(2,7),(7,2)]
--    sumaDeDosPrimos 16  ==  [(3,13),(5,11),(11,5),(13,3)]
--    sumaDeDosPrimos 17  ==  []
sumaDeDosPrimos1 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos1 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (< n) primes

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

sumasDeDosPrimos2 :: [Integer]
sumasDeDosPrimos2 =
  [n | n <- [1..], not (null (sumaDeDosPrimos2 n))]

-- (sumasDeDosPrimos2 n) es la lista de pares (x,y) de primos cuya suma
-- es n y tales que x <= y. Por ejemplo,
--    sumaDeDosPrimos2  9  ==  [(2,7)]
--    sumaDeDosPrimos2 16  ==  [(3,13),(5,11)]
--    sumaDeDosPrimos2 17  ==  []
sumaDeDosPrimos2 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos2 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (<= (n `div` 2)) primes

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

sumasDeDosPrimos3 :: [Integer]
sumasDeDosPrimos3 = filter esSumaDeDosPrimos3 [4..]

-- (esSumaDeDosPrimos3 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos3  9  ==  True
--    esSumaDeDosPrimos3 16  ==  True
--    esSumaDeDosPrimos3 17  ==  False
esSumaDeDosPrimos3 :: Integer -> Bool
esSumaDeDosPrimos3 n
  | odd n     = isPrime (n-2)
  | otherwise = any isPrime [n-x | x <- takeWhile (<= (n `div` 2)) primes]

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

-- Usando la conjetura de Goldbach que dice que "Todo número par mayor
-- que 2 puede escribirse como suma de dos números primos" .

sumasDeDosPrimos4 :: [Integer]
sumasDeDosPrimos4 = filter esSumaDeDosPrimos4 [4..]

-- (esSumaDeDosPrimos4 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos4  9  ==  True
--    esSumaDeDosPrimos4 16  ==  True
--    esSumaDeDosPrimos4 17  ==  False
esSumaDeDosPrimos4 :: Integer -> Bool
esSumaDeDosPrimos4 n = even n || isPrime (n-2)

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

-- La propiedad es
prop_sumasDeDosPrimos :: NonNegative Int -> Bool
prop_sumasDeDosPrimos (NonNegative n) =
  all (== sumasDeDosPrimos1 !! n)
      [sumasDeDosPrimos2 !! n,
       sumasDeDosPrimos3 !! n,
       sumasDeDosPrimos4 !! n]
  
-- La comprobación es
--    λ> quickCheck prop_sumasDeDosPrimos
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> sumasDeDosPrimos1 !! 5000
--    7994
--    (2.61 secs, 9,299,106,792 bytes)
--    λ> sumasDeDosPrimos2 !! 5000
--    7994
--    (1.48 secs, 5,190,651,760 bytes)
--    λ> sumasDeDosPrimos3 !! 5000
--    7994
--    (0.12 secs, 351,667,104 bytes)
--    λ> sumasDeDosPrimos4 !! 5000
--    7994
--    (0.04 secs, 63,464,320 bytes)
--
--    λ> sumasDeDosPrimos3 !! (5*10^4)
--    83674
--    (2.23 secs, 7,776,049,264 bytes)
--    λ> sumasDeDosPrimos4 !! (5*10^4)
--    83674
--    (0.34 secs, 1,183,604,984 bytes)

100

101

102

103

104

105

106

107

108

import Data.Numbers.Primes (isPrime, primes)

import Test.QuickCheck

-- 1ª solución

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

sumasDeDosPrimos1 :: [Integer]

sumasDeDosPrimos1 =

[n | n <- [1..], not (null (sumaDeDosPrimos1 n))]

-- (sumaDeDosPrimos1 n) es la lista de pares de primos cuya suma es

-- n. Por ejemplo,

-- sumaDeDosPrimos 9 == [(2,7),(7,2)]

-- sumaDeDosPrimos 16 == [(3,13),(5,11),(11,5),(13,3)]

-- sumaDeDosPrimos 17 == []

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

sumaDeDosPrimos1 n =

[(x,n-x) | x <- primosN, isPrime (n-x)]

where primosN = takeWhile (< n) primes

-- 2ª solución

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

sumasDeDosPrimos2 :: [Integer]

sumasDeDosPrimos2 =

[n | n <- [1..], not (null (sumaDeDosPrimos2 n))]

-- (sumasDeDosPrimos2 n) es la lista de pares (x,y) de primos cuya suma

-- es n y tales que x <= y. Por ejemplo,

-- sumaDeDosPrimos2 9 == [(2,7)]

-- sumaDeDosPrimos2 16 == [(3,13),(5,11)]

-- sumaDeDosPrimos2 17 == []

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

sumaDeDosPrimos2 n =

[(x,n-x) | x <- primosN, isPrime (n-x)]

where primosN = takeWhile (<= (n `div` 2)) primes

-- 3ª solución

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

sumasDeDosPrimos3 :: [Integer]

sumasDeDosPrimos3 = filter esSumaDeDosPrimos3 [4..]

-- (esSumaDeDosPrimos3 n) se verifica si n es suma de dos primos. Por

-- ejemplo,

-- esSumaDeDosPrimos3 9 == True

-- esSumaDeDosPrimos3 16 == True

-- esSumaDeDosPrimos3 17 == False

esSumaDeDosPrimos3 :: Integer -> Bool

esSumaDeDosPrimos3 n

| odd n = isPrime (n-2)

| otherwise = any isPrime [n-x | x <- takeWhile (<= (n `div` 2)) primes]

-- 4ª solución

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

-- Usando la conjetura de Goldbach que dice que "Todo número par mayor

-- que 2 puede escribirse como suma de dos números primos" .

sumasDeDosPrimos4 :: [Integer]

sumasDeDosPrimos4 = filter esSumaDeDosPrimos4 [4..]

-- (esSumaDeDosPrimos4 n) se verifica si n es suma de dos primos. Por

-- ejemplo,

-- esSumaDeDosPrimos4 9 == True

-- esSumaDeDosPrimos4 16 == True

-- esSumaDeDosPrimos4 17 == False

esSumaDeDosPrimos4 :: Integer -> Bool

esSumaDeDosPrimos4 n = even n || isPrime (n-2)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sumasDeDosPrimos :: NonNegative Int -> Bool

prop_sumasDeDosPrimos (NonNegative n) =

all (== sumasDeDosPrimos1 !! n)

[sumasDeDosPrimos2 !! n,

sumasDeDosPrimos3 !! n,

sumasDeDosPrimos4 !! n]

-- La comprobación es

-- λ> quickCheck prop_sumasDeDosPrimos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sumasDeDosPrimos1 !! 5000

-- 7994

-- (2.61 secs, 9,299,106,792 bytes)

-- λ> sumasDeDosPrimos2 !! 5000

-- 7994

-- (1.48 secs, 5,190,651,760 bytes)

-- λ> sumasDeDosPrimos3 !! 5000

-- 7994

-- (0.12 secs, 351,667,104 bytes)

-- λ> sumasDeDosPrimos4 !! 5000

-- 7994

-- (0.04 secs, 63,464,320 bytes)

-- λ> sumasDeDosPrimos3 !! (5*10^4)

-- 83674

-- (2.23 secs, 7,776,049,264 bytes)

-- λ> sumasDeDosPrimos4 !! (5*10^4)

-- 83674

-- (0.34 secs, 1,183,604,984 bytes)

El código se encuentra en GitHub.

Referencia

N.J.A. Sloane, Sucesión A014091 en OEIS.

La sucesión de Thue-Morse

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

La serie de Thue-Morse comienza con el término [0] y sus siguientes términos se construyen añadiéndole al anterior su complementario. Los primeros términos de la serie son

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

[0]

[0,1]

[0,1,1,0]

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

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

De esta forma se va formando una sucesión

   0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,...

1	0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,...

que se conoce como la sucesión de Thue-Morse.

Definir la sucesión

   sucThueMorse :: [Int]

1	sucThueMorse :: [Int]

cuyos elementos son los de la sucesión de Thue-Morse. Por ejemplo,

   λ> take 30 sucThueMorse
   [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0]
   λ> map (sucThueMorse4 !!) [1234567..1234596] 
   [1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0]
   λ> map (sucThueMorse4 !!) [4000000..4000030] 
   [1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1]

λ> take 30 sucThueMorse

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

λ> map (sucThueMorse4 !!) [1234567..1234596]

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

λ> map (sucThueMorse4 !!) [4000000..4000030]

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

Comprobar con QuickCheck que si s(n) representa el término n-ésimo de la sucesión de Thue-Morse, entonces

   s(2n)   = s(n)
   s(2n+1) = 1 - s(n)

1 2	s(2n) = s(n) s(2n+1) = 1 - s(n)

Soluciones

import Test.QuickCheck

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

sucThueMorse1 :: [Int]
sucThueMorse1 = map termSucThueMorse1 [0..]

-- (termSucThueMorse1 n) es el n-ésimo término de la sucesión de
-- Thue-Morse. Por ejemplo, 
--    termSucThueMorse1 0  ==  0
--    termSucThueMorse1 1  ==  1
--    termSucThueMorse1 2  ==  1
--    termSucThueMorse1 3  ==  0
--    termSucThueMorse1 4  ==  1
termSucThueMorse1 :: Int -> Int
termSucThueMorse1 0 = 0
termSucThueMorse1 n = 
  (serieThueMorse !! k) !! n
  where k = 1 + floor (logBase 2 (fromIntegral n))

-- serieThueMorse es la lista cuyos elementos son los términos de la
-- serie de Thue-Morse. Por ejemplo, 
--    λ> take 4 serieThueMorse3
--    [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]
serieThueMorse :: [[Int]]
serieThueMorse = iterate paso [0]
  where paso xs = xs ++ map (1-) xs

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

sucThueMorse2 :: [Int]
sucThueMorse2 = 
  0 : intercala (map (1-) sucThueMorse2) (tail sucThueMorse2)

-- (intercala xs ys) es la lista obtenida intercalando los elementos de
-- las listas infinitas xs e ys. Por ejemplo, 
--    take 10 (intercala [1,5..] [2,4..])  ==  [1,2,5,4,9,6,13,8,17,10]
intercala :: [a] -> [a] -> [a]
intercala (x:xs) ys = x : intercala ys xs 

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

sucThueMorse3 :: [Int]
sucThueMorse3 = 0 : 1 : aux (tail sucThueMorse3) 
  where aux (x : xs) = x : (1 - x) : aux xs



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

sucThueMorse4 :: [Int]
sucThueMorse4 = 0 : aux [1]
  where aux xs = xs ++ aux (xs ++ map (1-) xs) 

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_termSucThueMorse :: NonNegative Int -> Bool
prop_termSucThueMorse (NonNegative n) =
  sucThueMorse1 !! (2*n)   == sn &&
  sucThueMorse1 !! (2*n+1) == 1 - sn 
  where sn = sucThueMorse1 !! n

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

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

sucThueMorse5 :: [Int]
sucThueMorse5 = map termSucThueMorse5 [0..]

-- (termSucThueMorse5 n) es el n-ésimo término de la sucesión de
-- Thue-Morse. Por ejemplo, 
--    termSucThueMorse5 0  ==  0
--    termSucThueMorse5 1  ==  1
--    termSucThueMorse5 2  ==  1
--    termSucThueMorse5 3  ==  0
--    termSucThueMorse5 4  ==  1
termSucThueMorse5 :: Int -> Int
termSucThueMorse5 0 = 0
termSucThueMorse5 n 
  | even n    = termSucThueMorse5 (n `div` 2)
  | otherwise = 1 - termSucThueMorse5 (n `div` 2)

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

-- La propiedad es
prop_sucThueMorse :: NonNegative Int -> Bool
prop_sucThueMorse (NonNegative n) =
  all (== sucThueMorse1 !! n)
      [sucThueMorse2 !! n,
       sucThueMorse3 !! n,
       sucThueMorse4 !! n,
       sucThueMorse5 !! n]

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

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

-- La comparación es
--    λ> sucThueMorse1 !! (10^7)
--    0
--    (3.28 secs, 3,420,080,168 bytes)
--    λ> sucThueMorse2 !! (10^7)
--    0
--    (3.01 secs, 1,720,549,640 bytes)
--    λ> sucThueMorse3 !! (10^7)
--    0
--    (1.80 secs, 1,360,550,040 bytes)
--    λ> sucThueMorse4 !! (10^7)
--    0
--    (0.88 secs, 1,254,772,768 bytes)
--    λ> sucThueMorse5 !! (10^7)
--    0
--    (0.62 secs, 1,600,557,072 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

import Test.QuickCheck

-- 1ª solución

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

sucThueMorse1 :: [Int]

sucThueMorse1 = map termSucThueMorse1 [0..]

-- (termSucThueMorse1 n) es el n-ésimo término de la sucesión de

-- Thue-Morse. Por ejemplo,

-- termSucThueMorse1 0 == 0

-- termSucThueMorse1 1 == 1

-- termSucThueMorse1 2 == 1

-- termSucThueMorse1 3 == 0

-- termSucThueMorse1 4 == 1

termSucThueMorse1 :: Int -> Int

termSucThueMorse1 0 = 0

termSucThueMorse1 n =

(serieThueMorse !! k) !! n

where k = 1 + floor (logBase 2 (fromIntegral n))

-- serieThueMorse es la lista cuyos elementos son los términos de la

-- serie de Thue-Morse. Por ejemplo,

-- λ> take 4 serieThueMorse3

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

serieThueMorse :: [[Int]]

serieThueMorse = iterate paso [0]

where paso xs = xs ++ map (1-) xs

-- 2ª solución

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

sucThueMorse2 :: [Int]

sucThueMorse2 =

0 : intercala (map (1-) sucThueMorse2) (tail sucThueMorse2)

-- (intercala xs ys) es la lista obtenida intercalando los elementos de

-- las listas infinitas xs e ys. Por ejemplo,

-- take 10 (intercala [1,5..] [2,4..]) == [1,2,5,4,9,6,13,8,17,10]

intercala :: [a] -> [a] -> [a]

intercala (x:xs) ys = x : intercala ys xs

-- 3ª solución

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

sucThueMorse3 :: [Int]

sucThueMorse3 = 0 : 1 : aux (tail sucThueMorse3)

where aux (x : xs) = x : (1 - x) : aux xs

-- 4ª solución

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

sucThueMorse4 :: [Int]

sucThueMorse4 = 0 : aux [1]

where aux xs = xs ++ aux (xs ++ map (1-) xs)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_termSucThueMorse :: NonNegative Int -> Bool

prop_termSucThueMorse (NonNegative n) =

sucThueMorse1 !! (2*n) == sn &&

sucThueMorse1 !! (2*n+1) == 1 - sn

where sn = sucThueMorse1 !! n

-- La comprobación es

-- λ> quickCheck prop_termSucThueMorse

-- +++ OK, passed 100 tests.

-- 5ª solución

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

sucThueMorse5 :: [Int]

sucThueMorse5 = map termSucThueMorse5 [0..]

-- (termSucThueMorse5 n) es el n-ésimo término de la sucesión de

-- Thue-Morse. Por ejemplo,

-- termSucThueMorse5 0 == 0

-- termSucThueMorse5 1 == 1

-- termSucThueMorse5 2 == 1

-- termSucThueMorse5 3 == 0

-- termSucThueMorse5 4 == 1

termSucThueMorse5 :: Int -> Int

termSucThueMorse5 0 = 0

termSucThueMorse5 n

| even n = termSucThueMorse5 (n `div` 2)

| otherwise = 1 - termSucThueMorse5 (n `div` 2)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sucThueMorse :: NonNegative Int -> Bool

prop_sucThueMorse (NonNegative n) =

all (== sucThueMorse1 !! n)

[sucThueMorse2 !! n,

sucThueMorse3 !! n,

sucThueMorse4 !! n,

sucThueMorse5 !! n]

-- La comprobación es

-- λ> quickCheck prop_sucThueMorse

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sucThueMorse1 !! (10^7)

-- 0

-- (3.28 secs, 3,420,080,168 bytes)

-- λ> sucThueMorse2 !! (10^7)

-- 0

-- (3.01 secs, 1,720,549,640 bytes)

-- λ> sucThueMorse3 !! (10^7)

-- 0

-- (1.80 secs, 1,360,550,040 bytes)

-- λ> sucThueMorse4 !! (10^7)

-- 0

-- (0.88 secs, 1,254,772,768 bytes)

-- λ> sucThueMorse5 !! (10^7)

-- 0

-- (0.62 secs, 1,600,557,072 bytes)

El código se encuentra en GitHub.

Referencias

N.J.A. Sloane Sucesión A010060 en OEIS.
Programming Praxis Thue-Morse sequence.
Wikipedia Thue–Morse sequence

La serie de Thue-Morse

PorJosé A. Alonso 27-julio-202220-julio-2022

La serie de Thue-Morse comienza con el término [0] y sus siguientes términos se construyen añadiéndole al anterior su complementario (es decir, la lista obtenida cambiando el 0 por 1 y el 1 por 0). Los primeros términos de la serie son

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

[0]

[0,1]

[0,1,1,0]

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

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

Definir la lista

   serieThueMorse :: [[Int]]

1	serieThueMorse :: [[Int]]

tal que sus elementos son los términos de la serie de Thue-Morse. Por ejemplo,

   λ> take 4 serieThueMorse
   [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]

1 2	λ> take 4 serieThueMorse [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]

Soluciones

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

serieThueMorse1 :: [[Int]]
serieThueMorse1 = map termSerieThueMorse [0..]

-- (termSerieThueMorse n) es el término n-ésimo de la serie de
-- Thue-Morse. Por ejemplo, 
--    termSerieThueMorse 1  ==  [0,1]
--    termSerieThueMorse 2  ==  [0,1,1,0]
--    termSerieThueMorse 3  ==  [0,1,1,0,1,0,0,1]
--    termSerieThueMorse 4  ==  [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0]
termSerieThueMorse :: Int -> [Int]
termSerieThueMorse 0 = [0]
termSerieThueMorse n = xs ++ complementaria xs
  where xs = termSerieThueMorse (n-1)

-- (complementaria xs) es la complementaria de la lista xs (formada por
-- ceros y unos); es decir, la lista obtenida cambiando el 0 por 1 y el
-- 1 por 0. Por ejemplo, 
--    complementaria [1,0,0,1,1,0,1]  ==  [0,1,1,0,0,1,0]
complementaria :: [Int] -> [Int]
complementaria = map (1-)

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

serieThueMorse2 :: [[Int]]
serieThueMorse2 = [0] : map paso serieThueMorse2
  where paso xs = xs ++ complementaria xs

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

serieThueMorse3 :: [[Int]]
serieThueMorse3 = iterate paso [0]
  where paso xs = xs ++ complementaria xs

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

-- Observando que cada término de la serie de Thue-Morse se obtiene del
-- anterior sustituyendo los 1 por 1, 0 y los 0 por  0, 1. 

serieThueMorse4 :: [[Int]]
serieThueMorse4 = [0] : map (concatMap paso4) serieThueMorse4
  where paso4 x = [x,1-x]

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

-- La propiedad es
prop_serieThueMorse :: NonNegative Int -> Bool 
prop_serieThueMorse  (NonNegative n) =
  all (== serieThueMorse1 !! n)
      [serieThueMorse2 !! n,
       serieThueMorse3 !! n,
       serieThueMorse4 !! n]

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

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

-- La comparación es
--    λ> (serieThueMorse1 !! 23) !! (2^22)
--    1
--    (1.08 secs, 839,419,224 bytes)
--    λ> (serieThueMorse2 !! 23) !! (2^22)
--    1
--    (0.61 secs, 839,413,592 bytes)
--    λ> (serieThueMorse3 !! 23) !! (2^22)
--    1
--    (1.43 secs, 839,413,592 bytes)
--    λ> (serieThueMorse4 !! 23) !! (2^22)
--    1
--    (1.57 secs, 1,007,190,024 bytes)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

serieThueMorse1 :: [[Int]]

serieThueMorse1 = map termSerieThueMorse [0..]

-- (termSerieThueMorse n) es el término n-ésimo de la serie de

-- Thue-Morse. Por ejemplo,

-- termSerieThueMorse 1 == [0,1]

-- termSerieThueMorse 2 == [0,1,1,0]

-- termSerieThueMorse 3 == [0,1,1,0,1,0,0,1]

-- termSerieThueMorse 4 == [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0]

termSerieThueMorse :: Int -> [Int]

termSerieThueMorse 0 = [0]

termSerieThueMorse n = xs ++ complementaria xs

where xs = termSerieThueMorse (n-1)

-- (complementaria xs) es la complementaria de la lista xs (formada por

-- ceros y unos); es decir, la lista obtenida cambiando el 0 por 1 y el

-- 1 por 0. Por ejemplo,

-- complementaria [1,0,0,1,1,0,1] == [0,1,1,0,0,1,0]

complementaria :: [Int] -> [Int]

complementaria = map (1-)

-- 2ª solución

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

serieThueMorse2 :: [[Int]]

serieThueMorse2 = [0] : map paso serieThueMorse2

where paso xs = xs ++ complementaria xs

-- 3ª solución

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

serieThueMorse3 :: [[Int]]

serieThueMorse3 = iterate paso [0]

where paso xs = xs ++ complementaria xs

-- 4ª solución

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

-- Observando que cada término de la serie de Thue-Morse se obtiene del

-- anterior sustituyendo los 1 por 1, 0 y los 0 por 0, 1.

serieThueMorse4 :: [[Int]]

serieThueMorse4 = [0] : map (concatMap paso4) serieThueMorse4

where paso4 x = [x,1-x]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_serieThueMorse :: NonNegative Int -> Bool

prop_serieThueMorse (NonNegative n) =

all (== serieThueMorse1 !! n)

[serieThueMorse2 !! n,

serieThueMorse3 !! n,

serieThueMorse4 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_serieThueMorse

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (serieThueMorse1 !! 23) !! (2^22)

-- 1

-- (1.08 secs, 839,419,224 bytes)

-- λ> (serieThueMorse2 !! 23) !! (2^22)

-- 1

-- (0.61 secs, 839,413,592 bytes)

-- λ> (serieThueMorse3 !! 23) !! (2^22)

-- 1

-- (1.43 secs, 839,413,592 bytes)

-- λ> (serieThueMorse4 !! 23) !! (2^22)

-- 1

-- (1.57 secs, 1,007,190,024 bytes)

El código se encuentra en [GitHub](https://github.com/jaalonso/Exercitium/blob/main/src/

Referencias

N.J.A. Sloane Sucesión A010060 en OEIS.
Programming Praxis Thue-Morse sequence.
Wikipedia Thue–Morse sequence

Números belgas

PorJosé A. Alonso 26-julio-202219-julio-2022

Un número n es k-belga si la sucesión cuyo primer elemento es k y cuyos elementos se obtienen sumando reiteradamente las cifras de n contiene a n.

El 18 es 0-belga, porque a partir del 0 vamos a ir sumando sucesivamente 1, 8, 1, 8, … hasta llegar o sobrepasar el 18: 0, 1, 9, 10, 18, … Como se alcanza el 18, resulta que el 18 es 0-belga.

El 19 no es 1-belga, porque a partir del 1 vamos a ir sucesivamente 1, 9, 1, 9, … hasta llegar o sobrepasar el 18: 0, 1, 10, 11, 20, 21, … Como no se alcanza el 19, resulta que el 19 no es 1-belga.

Definir la función

   esBelga :: Int -> Int -> Bool

1	esBelga :: Int -> Int -> Bool

tal que (esBelga k n) se verifica si n es k-belga. Por ejemplo,

   esBelga 0 18                              ==  True
   esBelga 1 19                              ==  False
   esBelga 0 2016                            ==  True
   [x | x <- [0..30], esBelga 7 x]           ==  [7,10,11,21,27,29]
   [x | x <- [0..30], esBelga 10 x]          ==  [10,11,20,21,22,24,26]
   length [n | n <- [1..10^6], esBelga 0 n]  ==  272049

esBelga 0 18 == True

esBelga 1 19 == False

esBelga 0 2016 == True

[x | x <- [0..30], esBelga 7 x] == [7,10,11,21,27,29]

[x | x <- [0..30], esBelga 10 x] == [10,11,20,21,22,24,26]

length [n | n <- [1..10^6], esBelga 0 n] == 272049

Comprobar con QuickCheck que para todo número entero positivo n, si k es el resto de n entre la suma de los dígitos de n, entonces n es k-belga.

Soluciones

import Data.Char (digitToInt)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

esBelga1 :: Int -> Int -> Bool
esBelga1 k n =
  n == head (dropWhile (<n) (scanl (+) k (cycle (digitos n))))

digitos :: Int -> [Int]
digitos n = map digitToInt (show n)

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

esBelga2 :: Int -> Int -> Bool
esBelga2 k n =
  k <= n && n == head (dropWhile (<n) (scanl (+) (k + q * s) ds))
  where ds = digitos n
        s  = sum ds
        q  = (n - k) `div` s

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

-- La propiedad es
prop_esBelga :: Positive Int -> Positive Int -> Bool
prop_esBelga (Positive k) (Positive n) = 
  esBelga1 k n == esBelga2 k n

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

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

-- La comparación es
--    λ> length [n | n <- [1..2*10^4], esBelga1 0 n]
--    6521
--    (6.27 secs, 6,508,102,192 bytes)
--    λ> length [n | n <- [1..2*10^4], esBelga2 0 n]
--    6521
--    (0.07 secs, 46,741,144 bytes)

-- Verificación de la propiedad
-- ============================

-- La propiedad es
prop_Belga :: Positive Int -> Bool
prop_Belga (Positive n) = 
  esBelga2 k n
  where k = n `mod` sum (digitos n)

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

import Data.Char (digitToInt)

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

-- 1ª solución

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

esBelga1 :: Int -> Int -> Bool

esBelga1 k n =

n == head (dropWhile (<n) (scanl (+) k (cycle (digitos n))))

digitos :: Int -> [Int]

digitos n = map digitToInt (show n)

-- 2ª solución

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

esBelga2 :: Int -> Int -> Bool

esBelga2 k n =

k <= n && n == head (dropWhile (<n) (scanl (+) (k + q * s) ds))

where ds = digitos n

s = sum ds

q = (n - k) `div` s

-- Equivalencia

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

-- La propiedad es

prop_esBelga :: Positive Int -> Positive Int -> Bool

prop_esBelga (Positive k) (Positive n) =

esBelga1 k n == esBelga2 k n

-- La comprobación es

-- λ> quickCheck prop_esBelga

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length [n | n <- [1..2*10^4], esBelga1 0 n]

-- 6521

-- (6.27 secs, 6,508,102,192 bytes)

-- λ> length [n | n <- [1..2*10^4], esBelga2 0 n]

-- 6521

-- (0.07 secs, 46,741,144 bytes)

-- Verificación de la propiedad

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

-- La propiedad es

prop_Belga :: Positive Int -> Bool

prop_Belga (Positive n) =

esBelga2 k n

where k = n `mod` sum (digitos n)

-- La comprobación es

-- λ> quickCheck prop_Belga

-- +++ OK, passed 100 tests.

Referencias

Basado en el artículo Números belgas del blog Números y hoja de cálculo de Antonio Roldán Martínez.

El código se encuentra en GitHub.

Huecos maximales entre primos

PorJosé A. Alonso 25-julio-202218-julio-2022

El hueco de un número primo p es la distancia entre p y primo siguiente de p. Por ejemplo, el hueco de 7 es 4 porque el primo siguiente de 7 es 11 y 4 = 11-7. Los huecos de los primeros números son

   Primo Hueco
    2    1
    3    2
    7    4
   11    2

Primo Hueco

2 1

3 2

7 4

11 2

El hueco de un número primo p es maximal si es mayor que huecos de todos los números menores que p. Por ejemplo, 4 es un hueco maximal de 7 ya que los huecos de los primos menores que 7 son 1 y 2 y ambos son menores que 4. La tabla de los primeros huecos maximales es

   Primo Hueco
     2    1
     3    2
     7    4
    23    6
    89    8
   113   14
   523   18
   887   20

Primo Hueco

2 1

3 2

7 4

23 6

89 8

113 14

523 18

887 20

Definir la sucesión

   primosYhuecosMaximales :: [(Integer,Integer)]

1	primosYhuecosMaximales :: [(Integer,Integer)]

cuyos elementos son los números primos con huecos maximales junto son sus huecos. Por ejemplo,

   λ> take 8 primosYhuecosMaximales
   [(2,1),(3,2),(7,4),(23,6),(89,8),(113,14),(523,18),(887,20)]
   λ> primosYhuecosMaximales !! 20
   (2010733,148)

λ> take 8 primosYhuecosMaximales

[(2,1),(3,2),(7,4),(23,6),(89,8),(113,14),(523,18),(887,20)]

λ> primosYhuecosMaximales !! 20

(2010733,148)

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

primosYhuecosMaximales1 :: [(Integer,Integer)]
primosYhuecosMaximales1 = 
  [(p,huecoPrimo p) | p <- primes, esMaximalHuecoPrimo p]

-- (siguientePrimo x) es el menor primo mayor que x. Por ejemplo,
--    siguientePrimo 7  ==  11
--    siguientePrimo 8  ==  11
siguientePrimo :: Integer -> Integer
siguientePrimo p =
  head (dropWhile (<= p) primes)

-- (huecoPrimo p) es la distancia del primo p hasta el siguiente
-- primo. Por ejemplo,
--    huecoPrimo 7  ==  4
huecoPrimo :: Integer -> Integer
huecoPrimo p = siguientePrimo p - p

-- (esMaximalHuecoPrimo p) se verifica si el hueco primo de p es
-- maximal. Por ejemplo,
--    esMaximalHuecoPrimo  7  ==  True
--    esMaximalHuecoPrimo 11  ==  False
esMaximalHuecoPrimo :: Integer -> Bool
esMaximalHuecoPrimo p =
  and [huecoPrimo n < h | n <- takeWhile (< p) primes]
  where h = huecoPrimo p

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

primosYhuecosMaximales2 :: [(Integer,Integer)]
primosYhuecosMaximales2 = aux primosYhuecos
  where aux ((x,y):ps) = (x,y) : aux (dropWhile (\(_,b) -> b <= y) ps)

-- primosYhuecos es la lista de los números primos junto son sus
-- huecos. Por ejemplo, 
--    λ> take 10 primosYhuecos
--    [(2,1),(3,2),(5,2),(7,4),(11,2),(13,4),(17,2),(19,4),(23,6),(29,2)]
primosYhuecos :: [(Integer,Integer)]
primosYhuecos =
  [(x,y-x) | (x,y) <- zip primes (tail primes)]

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

primosYhuecosMaximales3 :: [(Integer,Integer)]
primosYhuecosMaximales3 = aux 0 primes
  where aux n (x:y:zs) | y-x > n   = (x,y-x) : aux (y-x) (y:zs)
                       | otherwise = aux n (y:zs)

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

-- La propiedad es
prop_primosYhuecosMaximales :: NonNegative Int -> Bool
prop_primosYhuecosMaximales (NonNegative n) =
  all (== primosYhuecosMaximales1 !! n)
      [primosYhuecosMaximales2 !! n,
       primosYhuecosMaximales3 !! n]

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

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

-- La comparación es
--    λ> primosYhuecosMaximales1 !! 10
--    (9551,36)
--    (2.63 secs, 7,400,316,112 bytes)
--    λ> primosYhuecosMaximales2 !! 10
--    (9551,36)
--    (0.01 secs, 7,060,744 bytes)
--    λ> primosYhuecosMaximales3 !! 10
--    (9551,36)
--    (0.01 secs, 4,000,368 bytes)
--    
--    λ> primosYhuecosMaximales2 !! 22
--    (17051707,180)
--    (7.90 secs, 17,275,407,712 bytes)
--    λ> primosYhuecosMaximales3 !! 22
--    (17051707,180)
--    (3.78 secs, 8,808,779,096 bytes)

import Data.Numbers.Primes (primes)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

primosYhuecosMaximales1 :: [(Integer,Integer)]

primosYhuecosMaximales1 =

[(p,huecoPrimo p) | p <- primes, esMaximalHuecoPrimo p]

-- (siguientePrimo x) es el menor primo mayor que x. Por ejemplo,

-- siguientePrimo 7 == 11

-- siguientePrimo 8 == 11

siguientePrimo :: Integer -> Integer

siguientePrimo p =

head (dropWhile (<= p) primes)

-- (huecoPrimo p) es la distancia del primo p hasta el siguiente

-- primo. Por ejemplo,

-- huecoPrimo 7 == 4

huecoPrimo :: Integer -> Integer

huecoPrimo p = siguientePrimo p - p

-- (esMaximalHuecoPrimo p) se verifica si el hueco primo de p es

-- maximal. Por ejemplo,

-- esMaximalHuecoPrimo 7 == True

-- esMaximalHuecoPrimo 11 == False

esMaximalHuecoPrimo :: Integer -> Bool

esMaximalHuecoPrimo p =

and [huecoPrimo n < h | n <- takeWhile (< p) primes]

where h = huecoPrimo p

-- 2ª solución

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

primosYhuecosMaximales2 :: [(Integer,Integer)]

primosYhuecosMaximales2 = aux primosYhuecos

where aux ((x,y):ps) = (x,y) : aux (dropWhile (\(_,b) -> b <= y) ps)

-- primosYhuecos es la lista de los números primos junto son sus

-- huecos. Por ejemplo,

-- λ> take 10 primosYhuecos

-- [(2,1),(3,2),(5,2),(7,4),(11,2),(13,4),(17,2),(19,4),(23,6),(29,2)]

primosYhuecos :: [(Integer,Integer)]

primosYhuecos =

[(x,y-x) | (x,y) <- zip primes (tail primes)]

-- 3ª solución

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

primosYhuecosMaximales3 :: [(Integer,Integer)]

primosYhuecosMaximales3 = aux 0 primes

where aux n (x:y:zs) | y-x > n = (x,y-x) : aux (y-x) (y:zs)

| otherwise = aux n (y:zs)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosYhuecosMaximales :: NonNegative Int -> Bool

prop_primosYhuecosMaximales (NonNegative n) =

all (== primosYhuecosMaximales1 !! n)

[primosYhuecosMaximales2 !! n,

primosYhuecosMaximales3 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=12}) prop_primosYhuecosMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosYhuecosMaximales1 !! 10

-- (9551,36)

-- (2.63 secs, 7,400,316,112 bytes)

-- λ> primosYhuecosMaximales2 !! 10

-- (9551,36)

-- (0.01 secs, 7,060,744 bytes)

-- λ> primosYhuecosMaximales3 !! 10

-- (9551,36)

-- (0.01 secs, 4,000,368 bytes)

-- λ> primosYhuecosMaximales2 !! 22

-- (17051707,180)

-- (7.90 secs, 17,275,407,712 bytes)

-- λ> primosYhuecosMaximales3 !! 22

-- (17051707,180)

-- (3.78 secs, 8,808,779,096 bytes)

Referencias

Basado en el ejercicio Maximal prime gaps de
Programming Praxis.

Otras referencias

C. Caldwell The gaps between primes.
J.K. Andersen Maximal prime gaps.
N.J.A. Sloane Sequence A002386 en OEIS.
N.J.A. Sloane Sequence A005250 en OEIS.
E.W. Weisstein Prime gaps en MathWorld.

El código se encuentra en GitHub.

La función indicatriz de Euler

PorJosé A. Alonso 22-julio-202217-julio-2022

La indicatriz de Euler (también función φ de Euler) es una función importante en teoría de números. Si n es un entero positivo, entonces φ(n) se define como el número de enteros positivos menores o iguales a n y coprimos con n. Por ejemplo, φ(36) = 12 ya que los números menores o iguales a 36 y coprimos con 36 son doce: 1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, y 35.

Definir la función

   phi :: Integer -> Integer

1	phi :: Integer -> Integer

tal que (phi n) es igual a φ(n). Por ejemplo,

   phi 36                          ==  12
   map phi [10..20]                ==  [4,10,4,12,6,8,8,16,6,18,8]
   phi (3^10^5) `mod` (10^9)       ==  681333334
   length (show (phi (10^(10^5)))) == 100000

phi 36 == 12

map phi [10..20] == [4,10,4,12,6,8,8,16,6,18,8]

phi (3^10^5) `mod` (10^9) == 681333334

length (show (phi (10^(10^5)))) == 100000

Comprobar con QuickCheck que, para todo n > 0, φ(10^n) tiene n dígitos.

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Math.NumberTheory.ArithmeticFunctions (totient)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

phi1 :: Integer -> Integer
phi1 n = genericLength [x | x <- [1..n], gcd x n == 1] 

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

phi2 :: Integer -> Integer
phi2 n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n] 
    
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

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

phi3 :: Integer -> Integer
phi3 n = 
  product [(x-1) * product xs | (x:xs) <- group (primeFactors n)]

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

phi4 :: Integer -> Integer
phi4 = totient 

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

-- La propiedad es
prop_phi :: Positive Integer -> Bool 
prop_phi (Positive n) =
  all (== phi1 n)
      [phi2 n,
       phi3 n,
       phi4 n]

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

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

-- La comparación es
--    λ> phi1 (2*10^6)
--    800000
--    (2.49 secs, 2,117,853,856 bytes)
--    λ> phi2 (2*10^6)
--    800000
--    (0.02 secs, 565,664 bytes)
--    
--    λ> length (show (phi2 (10^100000)))
--    100000
--    (2.80 secs, 5,110,043,208 bytes)
--    λ> length (show (phi3 (10^100000)))
--    100000
--    (4.81 secs, 7,249,353,896 bytes)
--    λ> length (show (phi4 (10^100000)))
--    100000
--    (0.78 secs, 1,467,573,768 bytes)

-- Verificación de la propiedad
-- ============================

-- La propiedad es
prop_phi2 :: Positive Integer -> Bool
prop_phi2 (Positive n) =
  genericLength (show (phi4 (10^n))) == n

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

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Math.NumberTheory.ArithmeticFunctions (totient)

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

-- 1ª solución

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

phi1 :: Integer -> Integer

phi1 n = genericLength [x | x <- [1..n], gcd x n == 1]

-- 2ª solución

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

phi2 :: Integer -> Integer

phi2 n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n]

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

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 3ª solución

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

phi3 :: Integer -> Integer

phi3 n =

product [(x-1) * product xs | (x:xs) <- group (primeFactors n)]

-- 4ª solución

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

phi4 :: Integer -> Integer

phi4 = totient

-- Comprobación de equivalencia

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

-- La propiedad es

prop_phi :: Positive Integer -> Bool

prop_phi (Positive n) =

all (== phi1 n)

[phi2 n,

phi3 n,

phi4 n]

-- La comprobación es

-- λ> quickCheck prop_phi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> phi1 (2*10^6)

-- 800000

-- (2.49 secs, 2,117,853,856 bytes)

-- λ> phi2 (2*10^6)

-- 800000

-- (0.02 secs, 565,664 bytes)

-- λ> length (show (phi2 (10^100000)))

-- 100000

-- (2.80 secs, 5,110,043,208 bytes)

-- λ> length (show (phi3 (10^100000)))

-- 100000

-- (4.81 secs, 7,249,353,896 bytes)

-- λ> length (show (phi4 (10^100000)))

-- 100000

-- (0.78 secs, 1,467,573,768 bytes)

-- Verificación de la propiedad

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

-- La propiedad es

prop_phi2 :: Positive Integer -> Bool

prop_phi2 (Positive n) =

genericLength (show (phi4 (10^n))) == n

-- La comprobación es

-- λ> quickCheck prop_phi2

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Sucesión de suma de cuadrados de los dígitos

PorJosé A. Alonso 21-julio-202216-julio-2022

Definir la función

   sucSumaCuadradosDigitos :: Integer -> [Integer]

1	sucSumaCuadradosDigitos :: Integer -> [Integer]

tal que (sucSumaCuadradosDigitos n) es la sucesión cuyo primer término es n y los restantes se obtienen sumando los cuadrados de los dígitos de su término anterior. Por ejemplo,

   λ> take 20 (sucSumaCuadradosDigitos1 2000)
   [2000,4,16,37,58,89,145,42,20,4,16,37,58,89,145,42,20,4,16,37]
   λ> take 20 (sucSumaCuadradosDigitos 1976)
   [1976,167,86,100,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
   λ> sucSumaCuadradosDigitos 2000 !! (10^9)
   20

λ> take 20 (sucSumaCuadradosDigitos1 2000)

[2000,4,16,37,58,89,145,42,20,4,16,37,58,89,145,42,20,4,16,37]

λ> take 20 (sucSumaCuadradosDigitos 1976)

[1976,167,86,100,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

λ> sucSumaCuadradosDigitos 2000 !! (10^9)

Soluciones

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

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

sucSumaCuadradosDigitos1 :: Integer -> [Integer]
sucSumaCuadradosDigitos1 n =
  n : [sumaCuadradosDigitos x | x <- sucSumaCuadradosDigitos1 n]

-- (sumaCuadradosDigitos n) es la suma de los cuadrados de los dígitos
-- de n. Por ejemplo, 
--    sumaCuadradosDigitos 2016  ==  41
sumaCuadradosDigitos :: Integer -> Integer
sumaCuadradosDigitos n = sum (map (^2) (digitos n))

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

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

sucSumaCuadradosDigitos2 :: Integer -> [Integer]
sucSumaCuadradosDigitos2 = iterate sumaCuadradosDigitos

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

-- A partir de los cálculos con las definiciones anteriores, se observa
-- que para todo n (sucSumaCuadradosDigitos n) tiene una parte pura y
-- otra periódica. Por ejemplo, para n = 2016, 
--    λ> take 20 (sucSumaCuadradosDigitos 2016)
--    [2016,41,17,50,25,29,85,89,145,42,20,4,16,37,58,89,145,42,20,4]
-- la parte pura es
--    [2016,41,17,50,25,29,85]
-- y la parte periódica es
--    [89,145,42,20,4,16,37,58])

sucSumaCuadradosDigitos3 :: Integer -> [Integer]
sucSumaCuadradosDigitos3 n = xs ++ cycle ys
  where (xs,ys) = sucCompactaSumaCuadradosDigitos n

-- (sucCompactaSumaCuadradosDigitos n) es el par formado por la parte
-- pura y la periódica de (sucSumaCuadradosDigitos n). Por ejemplo, 
--    λ> sucCompactaSumaCuadradosDigitos 2016
--    ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])
--    λ> sucCompactaSumaCuadradosDigitos 1976
--    ([1976,167,86,100],[1])
sucCompactaSumaCuadradosDigitos :: Integer -> ([Integer],[Integer])
sucCompactaSumaCuadradosDigitos = 
  partePuraPeriodica . sucSumaCuadradosDigitos1

-- (partePuraPeriodica xs) es el par formado por la parte pura y la
-- periódica de xs. Por ejemplo,
--    λ> partePuraPeriodica (sucSumaCuadradosDigitos 2016)
--    ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])
--    λ> partePuraPeriodica (sucSumaCuadradosDigitos 1976)
--    ([1976,167,86,100],[1])
partePuraPeriodica :: [Integer] -> ([Integer],[Integer])
partePuraPeriodica = aux [] 
  where aux as (b:bs) | b `elem` as = span (/=b) (reverse as)
                      | otherwise = aux (b:as) bs

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

sucSumaCuadradosDigitos4 :: Integer -> [Integer]
sucSumaCuadradosDigitos4 1  = repeat 1
sucSumaCuadradosDigitos4 89 = cycle [89,145,42,20,4,16,37,58]
sucSumaCuadradosDigitos4 n  =
  n : sucSumaCuadradosDigitos4 (sumaCuadradosDigitos n)

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

-- La propiedad es
prop_sucSumaCuadradosDigitos :: Positive Integer -> NonNegative Int -> Bool
prop_sucSumaCuadradosDigitos (Positive n) (NonNegative k) =
  all (== sucSumaCuadradosDigitos1 n !! k)
      [sucSumaCuadradosDigitos2 n !! k,
       sucSumaCuadradosDigitos3 n !! k,
       sucSumaCuadradosDigitos4 n !! k]

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

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

-- La comparación es
--    λ> sucSumaCuadradosDigitos1 2000 !! (10^4)
--    20
--    (6.96 secs, 8,049,886,312 bytes)
--    λ> sucSumaCuadradosDigitos2 2000 !! (10^4)
--    20
--    (0.08 secs, 91,024,688 bytes)
--    λ> sucSumaCuadradosDigitos3 2000 !! (10^4)
--    20
--    (0.01 secs, 995,560 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (10^4)
--    20
--    (0.02 secs, 587,040 bytes)
--    
--    λ> sucSumaCuadradosDigitos2 2000 !! (3*10^5)
--    20
--    (1.96 secs, 2,715,501,416 bytes)
--    λ> sucSumaCuadradosDigitos3 2000 !! (3*10^5)
--    20
--    (0.02 secs, 995,872 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (3*10^5)
--    20
--    (0.02 secs, 587,352 bytes)
--    
--    λ> sucSumaCuadradosDigitos3 2000 !! (10^9)
--    20
--    (2.85 secs, 996,016 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (10^9)
--    20
--    (2.54 secs, 587,496 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 (Positive (Positive), NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

sucSumaCuadradosDigitos1 :: Integer -> [Integer]

sucSumaCuadradosDigitos1 n =

n : [sumaCuadradosDigitos x | x <- sucSumaCuadradosDigitos1 n]

-- (sumaCuadradosDigitos n) es la suma de los cuadrados de los dígitos

-- de n. Por ejemplo,

-- sumaCuadradosDigitos 2016 == 41

sumaCuadradosDigitos :: Integer -> Integer

sumaCuadradosDigitos n = sum (map (^2) (digitos n))

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

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

digitos :: Integer -> [Integer]

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

-- 2ª solución

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

sucSumaCuadradosDigitos2 :: Integer -> [Integer]

sucSumaCuadradosDigitos2 = iterate sumaCuadradosDigitos

-- 3ª solución

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

-- A partir de los cálculos con las definiciones anteriores, se observa

-- que para todo n (sucSumaCuadradosDigitos n) tiene una parte pura y

-- otra periódica. Por ejemplo, para n = 2016,

-- λ> take 20 (sucSumaCuadradosDigitos 2016)

-- [2016,41,17,50,25,29,85,89,145,42,20,4,16,37,58,89,145,42,20,4]

-- la parte pura es

-- [2016,41,17,50,25,29,85]

-- y la parte periódica es

-- [89,145,42,20,4,16,37,58])

sucSumaCuadradosDigitos3 :: Integer -> [Integer]

sucSumaCuadradosDigitos3 n = xs ++ cycle ys

where (xs,ys) = sucCompactaSumaCuadradosDigitos n

-- (sucCompactaSumaCuadradosDigitos n) es el par formado por la parte

-- pura y la periódica de (sucSumaCuadradosDigitos n). Por ejemplo,

-- λ> sucCompactaSumaCuadradosDigitos 2016

-- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])

-- λ> sucCompactaSumaCuadradosDigitos 1976

-- ([1976,167,86,100],[1])

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

sucCompactaSumaCuadradosDigitos =

partePuraPeriodica . sucSumaCuadradosDigitos1

-- (partePuraPeriodica xs) es el par formado por la parte pura y la

-- periódica de xs. Por ejemplo,

-- λ> partePuraPeriodica (sucSumaCuadradosDigitos 2016)

-- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])

-- λ> partePuraPeriodica (sucSumaCuadradosDigitos 1976)

-- ([1976,167,86,100],[1])

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

partePuraPeriodica = aux []

where aux as (b:bs) | b `elem` as = span (/=b) (reverse as)

| otherwise = aux (b:as) bs

-- 4ª solución

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

sucSumaCuadradosDigitos4 :: Integer -> [Integer]

sucSumaCuadradosDigitos4 1 = repeat 1

sucSumaCuadradosDigitos4 89 = cycle [89,145,42,20,4,16,37,58]

sucSumaCuadradosDigitos4 n =

n : sucSumaCuadradosDigitos4 (sumaCuadradosDigitos n)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sucSumaCuadradosDigitos :: Positive Integer -> NonNegative Int -> Bool

prop_sucSumaCuadradosDigitos (Positive n) (NonNegative k) =

all (== sucSumaCuadradosDigitos1 n !! k)

[sucSumaCuadradosDigitos2 n !! k,

sucSumaCuadradosDigitos3 n !! k,

sucSumaCuadradosDigitos4 n !! k]

-- La comprobación es

-- λ> quickCheck prop_sucSumaCuadradosDigitos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sucSumaCuadradosDigitos1 2000 !! (10^4)

-- 20

-- (6.96 secs, 8,049,886,312 bytes)

-- λ> sucSumaCuadradosDigitos2 2000 !! (10^4)

-- 20

-- (0.08 secs, 91,024,688 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (10^4)

-- 20

-- (0.01 secs, 995,560 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (10^4)

-- 20

-- (0.02 secs, 587,040 bytes)

-- λ> sucSumaCuadradosDigitos2 2000 !! (3*10^5)

-- 20

-- (1.96 secs, 2,715,501,416 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (3*10^5)

-- 20

-- (0.02 secs, 995,872 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (3*10^5)

-- 20

-- (0.02 secs, 587,352 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (10^9)

-- 20

-- (2.85 secs, 996,016 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (10^9)

-- 20

-- (2.54 secs, 587,496 bytes)

El código se encuentra en GitHub.

Potencias perfectas

PorJosé A. Alonso 20-julio-202213-julio-2022

Un número natural n es una potencia perfecta si existen dos números naturales m > 1 y k > 1 tales que n = m^k. Las primeras potencias perfectas son

   4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 
   36 = 6², 49 = 7², 64 = 2⁶, ...

1 2	4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 36 = 6², 49 = 7², 64 = 2⁶, ...

Definir la sucesión

   potenciasPerfectas :: [Integer]

1	potenciasPerfectas :: [Integer]

cuyos términos son las potencias perfectas. Por ejemplo,

   take 10 potenciasPerfectas  ==  [4,8,9,16,25,27,32,36,49,64]
   potenciasPerfectas !! 3000  ==  7778521

1 2	take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 3000 == 7778521

Definir el procedimiento

   grafica :: Int -> IO ()

1	grafica :: Int -> IO ()

tal que (grafica n) es la representación gráfica de las n primeras potencias perfectas. Por ejemplo, para (grafica 30) dibuja

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple (Attribute (Key, PNG), plotList)
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

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

potenciasPerfectas1 :: [Integer]
potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta = not . null. potenciasPerfectasDe 

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que 
-- x = a^b. Por ejemplo,
--    potenciasPerfectasDe 64  ==  [(2,6),(4,3),(8,2)]
--    potenciasPerfectasDe 72  ==  []
potenciasPerfectasDe :: Integer -> [(Integer,Integer)]
potenciasPerfectasDe n = 
  [(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]
         , k <- takeWhile (\x -> m^x <= n) [2..]
         , m^k == n]

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

potenciasPerfectas2 :: [Integer]
potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta2 36  ==  True
--    esPotenciaPerfecta2 72  ==  False
esPotenciaPerfecta2 :: Integer -> Bool
esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Int]
exponentes x = [length ys | ys <- group (primeFactors x)] 

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,
--    mcd [4,6,10]  ==  2
--    mcd [4,5,10]  ==  1
mcd :: [Int] -> Int
mcd = foldl1 gcd

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

potenciasPerfectas3 :: [Integer]
potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números
-- mayores que 1 con exponentes mayores que 1. Por ejemplo,
--    λ> map (take 3) (take 4 potencias)
--    [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]
potencias:: [[Integer]]
potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xss. Por ejemplo,
--    take 7 (mezclaTodas potencias)  ==  [4,8,9,16,25,27,32]
mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
  where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xs e ys. Por ejemplo,
--    take 7 (mezcla [2,5..] [4,6..])  ==  [2,4,5,6,8,10,11]
mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y  = x : mezcla xs (y:ys)
                     | x == y = x : mezcla xs ys
                     | x > y  = y : mezcla (x:xs) ys

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

-- La propiedad es
prop_potenciasPerfectas :: NonNegative Int -> Bool
prop_potenciasPerfectas (NonNegative n) =
  all (== potenciasPerfectas1 !! n)
      [potenciasPerfectas2 !! n,
       potenciasPerfectas3 !! n]

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

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

-- La comparación es
--    λ> potenciasPerfectas1 !! 200
--    28224
--    (10.56 secs, 8,434,647,368 bytes)
--    λ> potenciasPerfectas2 !! 200
--    28224
--    (0.36 secs, 825,040,416 bytes)
--    λ> potenciasPerfectas3 !! 200
--    28224
--    (0.05 secs, 7,474,280 bytes)
--    
--    λ> potenciasPerfectas2 !! 500
--    191844
--    (4.16 secs, 9,899,367,112 bytes)
--    λ> potenciasPerfectas3 !! 500
--    191844
--    (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución
potenciasPerfectas :: [Integer]
potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica
-- ======================

grafica :: Int -> IO ()
grafica n = 
  plotList [ Key Nothing
           , PNG "Potencias_perfectas.png"
           ]
           (take n potenciasPerfectas)

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

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

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

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

-- 1ª solución

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

potenciasPerfectas1 :: [Integer]

potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta = not . null. potenciasPerfectasDe

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que

-- x = a^b. Por ejemplo,

-- potenciasPerfectasDe 64 == [(2,6),(4,3),(8,2)]

-- potenciasPerfectasDe 72 == []

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

potenciasPerfectasDe n =

[(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]

, k <- takeWhile (\x -> m^x <= n) [2..]

, m^k == n]

-- 2ª solución

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

potenciasPerfectas2 :: [Integer]

potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta2 36 == True

-- esPotenciaPerfecta2 72 == False

esPotenciaPerfecta2 :: Integer -> Bool

esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Int]

exponentes x = [length ys | ys <- group (primeFactors x)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,

-- mcd [4,6,10] == 2

-- mcd [4,5,10] == 1

mcd :: [Int] -> Int

mcd = foldl1 gcd

-- 3ª solución

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

potenciasPerfectas3 :: [Integer]

potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números

-- mayores que 1 con exponentes mayores que 1. Por ejemplo,

-- λ> map (take 3) (take 4 potencias)

-- [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]

potencias:: [[Integer]]

potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xss. Por ejemplo,

-- take 7 (mezclaTodas potencias) == [4,8,9,16,25,27,32]

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

mezclaTodas = foldr1 xmezcla

where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xs e ys. Por ejemplo,

-- take 7 (mezcla [2,5..] [4,6..]) == [2,4,5,6,8,10,11]

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

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

| x == y = x : mezcla xs ys

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_potenciasPerfectas :: NonNegative Int -> Bool

prop_potenciasPerfectas (NonNegative n) =

all (== potenciasPerfectas1 !! n)

[potenciasPerfectas2 !! n,

potenciasPerfectas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_potenciasPerfectas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> potenciasPerfectas1 !! 200

-- 28224

-- (10.56 secs, 8,434,647,368 bytes)

-- λ> potenciasPerfectas2 !! 200

-- 28224

-- (0.36 secs, 825,040,416 bytes)

-- λ> potenciasPerfectas3 !! 200

-- 28224

-- (0.05 secs, 7,474,280 bytes)

-- λ> potenciasPerfectas2 !! 500

-- 191844

-- (4.16 secs, 9,899,367,112 bytes)

-- λ> potenciasPerfectas3 !! 500

-- 191844

-- (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución

potenciasPerfectas :: [Integer]

potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica

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

grafica :: Int -> IO ()

grafica n =

plotList [ Key Nothing

, PNG "Potencias_perfectas.png"

]

(take n potenciasPerfectas)

El código se encuentra en GitHub.

Sumas alternas de factoriales

PorJosé A. Alonso 19-julio-202212-julio-2022

Las primeras sumas alternas de los factoriales son números primos; en efecto,

   3! - 2! + 1! = 5
   4! - 3! + 2! - 1! = 19
   5! - 4! + 3! - 2! + 1! = 101
   6! - 5! + 4! - 3! + 2! - 1! = 619
   7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421
   8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

3! - 2! + 1! = 5

4! - 3! + 2! - 1! = 19

5! - 4! + 3! - 2! + 1! = 101

6! - 5! + 4! - 3! + 2! - 1! = 619

7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421

8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

son primos, pero

   9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

1	9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

no es primo.

Definir las funciones

   sumaAlterna         :: Integer -> Integer
   sumasAlternas       :: [Integer]
   conSumaAlternaPrima :: [Integer]

sumaAlterna :: Integer -> Integer

sumasAlternas :: [Integer]

conSumaAlternaPrima :: [Integer]

tales que

(sumaAlterna n) es la suma alterna de los factoriales desde n hasta 1. Por ejemplo,

     sumaAlterna 3  ==  5
     sumaAlterna 4  ==  19
     sumaAlterna 5  ==  101
     sumaAlterna 6  ==  619
     sumaAlterna 7  ==  4421
     sumaAlterna 8  ==  35899
     sumaAlterna 9  ==  326981
     sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumaAlterna 3 == 5

sumaAlterna 4 == 19

sumaAlterna 5 == 101

sumaAlterna 6 == 619

sumaAlterna 7 == 4421

sumaAlterna 8 == 35899

sumaAlterna 9 == 326981

sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumasAlternas es la sucesión de las sumas alternas de factoriales. Por ejemplo,

     λ> take 10 sumasAlternas1
     [0,1,1,5,19,101,619,4421,35899,326981]

1 2	λ> take 10 sumasAlternas1 [0,1,1,5,19,101,619,4421,35899,326981]

conSumaAlternaPrima es la sucesión de los números cuya suma alterna de factoriales es prima. Por ejemplo,

     λ> take 8 conSumaAlternaPrima
     [3,4,5,6,7,8,10,15]

1 2	λ> take 8 conSumaAlternaPrima [3,4,5,6,7,8,10,15]

Soluciones

import Data.List (genericTake)
import Data.Numbers.Primes (isPrime)
import Test.QuickCheck

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

sumaAlterna1 :: Integer -> Integer
sumaAlterna1 1 = 1
sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

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

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

sumaAlterna2 :: Integer -> Integer
sumaAlterna2 n =
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where
    signos | odd n     = concat (repeat [1,-1])
           | otherwise = concat (repeat [-1,1])

-- factoriales es la lista de los factoriales. Por ejemplo,
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

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

sumaAlterna3 :: Integer -> Integer
sumaAlterna3 n = 
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where signos | odd n     = cycle [1,-1]
               | otherwise = cycle [-1,1]

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

sumaAlterna4 :: Integer -> Integer
sumaAlterna4 n =
  foldl (flip (-)) 0 (scanl1 (*) [1..n])

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

-- La propiedad es
prop_sumaAlterna :: Positive Integer -> Bool 
prop_sumaAlterna (Positive n) =
  all (== sumaAlterna1 n)
      [sumaAlterna2 n,
       sumaAlterna3 n,
       sumaAlterna4 n]

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

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

-- La comparación es
--    λ> sumaAlterna1 4000 `mod` (10^6)
--    577019
--    (6.21 secs, 16,154,113,192 bytes)
--    λ> sumaAlterna2 4000 `mod` (10^6)
--    577019
--    (0.01 secs, 24,844,664 bytes)
--    
--    λ> sumaAlterna2 (5*10^4) `mod` (10^6)
--    577019
--    (1.81 secs, 4,729,583,864 bytes)
--    λ> sumaAlterna3 (5*10^4) `mod` (10^6)
--    577019
--    (0.89 secs, 4,725,983,928 bytes)
--    λ> sumaAlterna4 (5*10^4) `mod` (10^6)
--    577019
--    (0.70 secs, 4,710,770,592 bytes)

-- En lo que sigue se usa la 3ª definición
sumaAlterna :: Integer -> Integer
sumaAlterna = sumaAlterna3

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

sumasAlternas1 :: [Integer]
sumasAlternas1 =
  map sumaAlterna [0..]

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

sumasAlternas2 :: [Integer]
sumasAlternas2 =
  0 : zipWith (-) (tail factoriales) sumasAlternas2

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

sumasAlternas3 :: [Integer]
sumasAlternas3 =
  scanl (flip (-)) 0 $ scanl1 (*) [1..]

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

-- La propiedad es
prop_sumasAlternas :: NonNegative Int -> Bool
prop_sumasAlternas (NonNegative n) =
  all (== sumasAlternas1 !! n)
      [sumasAlternas2 !! n,
       sumasAlternas3 !! n]
  
-- La comprobación es
--    λ> quickCheck prop_sumasAlternas
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> length (show (sumasAlternas1 !! (5*10^4)))
--    213237
--    (4.90 secs, 4,731,620,600 bytes)
--    λ> length (show (sumasAlternas2 !! (5*10^4)))
--    213237
--    (2.39 secs, 4,726,820,456 bytes)
--    λ> length (show (sumasAlternas3 !! (5*10^4)))
--    213237
--    (1.78 secs, 4,726,820,280 bytes)

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

conSumaAlternaPrima1 :: [Integer]
conSumaAlternaPrima1 =
  [n | n <- [0..], isPrime (sumaAlterna n)]

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

conSumaAlternaPrima2 :: [Integer]
conSumaAlternaPrima2 =
  [x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

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

conSumaAlternaPrima3 :: [Integer]
conSumaAlternaPrima3 =
  filter (isPrime . sumaAlterna) [0..]

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

-- La propiedad es
prop_conSumaAlternaPrima :: NonNegative Int -> Bool
prop_conSumaAlternaPrima (NonNegative n) =
  all (== conSumaAlternaPrima1 !! n)
      [conSumaAlternaPrima2 !! n,
       conSumaAlternaPrima3 !! n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=5}) prop_conSumaAlternaPrima
--    +++ 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.List (genericTake)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª definición de sumaAlterna

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

sumaAlterna1 :: Integer -> Integer

sumaAlterna1 1 = 1

sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª definición de sumaAlterna

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

sumaAlterna2 :: Integer -> Integer

sumaAlterna2 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where

signos | odd n = concat (repeat [1,-1])

| otherwise = concat (repeat [-1,1])

-- factoriales es la lista de los factoriales. Por ejemplo,

-- take 7 factoriales == [1,1,2,6,24,120,720]

factoriales :: [Integer]

factoriales = 1 : scanl1 (*) [1..]

-- 3ª definición de sumaAlterna

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

sumaAlterna3 :: Integer -> Integer

sumaAlterna3 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where signos | odd n = cycle [1,-1]

| otherwise = cycle [-1,1]

-- 3ª definición de sumaAlterna

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

sumaAlterna4 :: Integer -> Integer

sumaAlterna4 n =

foldl (flip (-)) 0 (scanl1 (*) [1..n])

-- Comprobación de equivalencia de sumaAlterna

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

-- La propiedad es

prop_sumaAlterna :: Positive Integer -> Bool

prop_sumaAlterna (Positive n) =

all (== sumaAlterna1 n)

[sumaAlterna2 n,

sumaAlterna3 n,

sumaAlterna4 n]

-- La comprobación es

-- λ> quickCheck prop_sumaAlterna

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumaAlterna

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

-- La comparación es

-- λ> sumaAlterna1 4000 `mod` (10^6)

-- 577019

-- (6.21 secs, 16,154,113,192 bytes)

-- λ> sumaAlterna2 4000 `mod` (10^6)

-- 577019

-- (0.01 secs, 24,844,664 bytes)

-- λ> sumaAlterna2 (5*10^4) `mod` (10^6)

-- 577019

-- (1.81 secs, 4,729,583,864 bytes)

-- λ> sumaAlterna3 (5*10^4) `mod` (10^6)

-- 577019

-- (0.89 secs, 4,725,983,928 bytes)

-- λ> sumaAlterna4 (5*10^4) `mod` (10^6)

-- 577019

-- (0.70 secs, 4,710,770,592 bytes)

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

sumaAlterna :: Integer -> Integer

sumaAlterna = sumaAlterna3

-- 1ª definición de sumasAlternas

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

sumasAlternas1 :: [Integer]

sumasAlternas1 =

map sumaAlterna [0..]

-- 2ª definición de sumasAlternas

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

sumasAlternas2 :: [Integer]

sumasAlternas2 =

0 : zipWith (-) (tail factoriales) sumasAlternas2

-- 3ª definición de sumasAlternas

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

sumasAlternas3 :: [Integer]

sumasAlternas3 =

scanl (flip (-)) 0 $ scanl1 (*) [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sumasAlternas :: NonNegative Int -> Bool

prop_sumasAlternas (NonNegative n) =

all (== sumasAlternas1 !! n)

[sumasAlternas2 !! n,

sumasAlternas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_sumasAlternas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (show (sumasAlternas1 !! (5*10^4)))

-- 213237

-- (4.90 secs, 4,731,620,600 bytes)

-- λ> length (show (sumasAlternas2 !! (5*10^4)))

-- 213237

-- (2.39 secs, 4,726,820,456 bytes)

-- λ> length (show (sumasAlternas3 !! (5*10^4)))

-- 213237

-- (1.78 secs, 4,726,820,280 bytes)

-- 1ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima1 :: [Integer]

conSumaAlternaPrima1 =

[n | n <- [0..], isPrime (sumaAlterna n)]

-- 2ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima2 :: [Integer]

conSumaAlternaPrima2 =

[x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

-- 3ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima3 :: [Integer]

conSumaAlternaPrima3 =

filter (isPrime . sumaAlterna) [0..]

-- Comprobación de equivalencia de conSumaAlternaPrima

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

-- La propiedad es

prop_conSumaAlternaPrima :: NonNegative Int -> Bool

prop_conSumaAlternaPrima (NonNegative n) =

all (== conSumaAlternaPrima1 !! n)

[conSumaAlternaPrima2 !! n,

conSumaAlternaPrima3 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos con cubos

PorJosé A. Alonso 18-julio-202211-julio-2022

Un primo con cubo es un número primo p para el que existe algún entero positivo n tal que la expresión n^3 + n^2p es un cubo perfecto. Por ejemplo, 19 es un primo con cubo ya que 8^3 + 8^2×19 = 12^3.

Definir la sucesión

   primosConCubos :: [Integer]

1	primosConCubos :: [Integer]

tal que sus elementos son los primos con cubo. Por ejemplo,

   λ> take 6 primosConCubos
   [7,19,37,61,127,271]
   λ> length (takeWhile (< 1000000) primosConCubos)
   173

λ> take 6 primosConCubos

[7,19,37,61,127,271]

λ> length (takeWhile (< 1000000) primosConCubos)

173

Soluciones

import Data.Numbers.Primes (isPrime)
import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs) 

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

primosConCubos1 :: [Integer]
primosConCubos1 =
  [p | x <- [1..],
       n <- [1..x],
       (x^3 - n^3) `mod` (n^2) == 0,
       let p = (x^3 - n^3) `div` (n^2),
       isPrime p]

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

-- Para analizar la respuesta, en esta solución se calculan los pares
-- (p,n) tales que p es un primo con cubo y n es un número positivo tal
-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]
primosConCubos2' =
  [(p,n) | x <- [1..],
           n <- [1..x],
           (x^3 - n^3) `mod` (n^2) == 0,
           let p = (x^3 - n^3) `div` (n^2),
           isPrime p]

-- El cálculo es
--    λ> take 7 primosConCubos2
--    [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]
-- es la de los cubos y que los primeros elementos se obtienen restando
-- los segundos elementos consecutivos; es decir,
--     7 =  8 -  1 = 2^3 - 1^3  
--    19 = 27 -  8 = 3^3 - 2^3
--    37 = 64 - 27 = 4^3 - 3^3
-- Continuando el patrón,
--     61 =  5^3 - 4^3 =  125 -   64
--     91 =  6^3 - 5^3 =  216 -  125
--    127 =  7^3 - 6^3 =  343 -  216
--    271 = 10^3 - 9^3 = 1000 -  729
--    331 = 11^3 -10^3 = 1331 - 1000
-- Por tanto, los primos con cubos son diferencias de dos cubos
-- consecutivos; es decir, coinciden con los números cubanos del
-- ejercicio anterior. A partir de la conjetura se obtienen las
-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente
-- definición 
primosConCubos2 :: [Integer]
primosConCubos2 = 
  filter isPrime [(x+1)^3 - x^3 | x <- [1..]] 

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

primosConCubos3 :: [Integer]
primosConCubos3 =
  filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]
diferenciasCubosConsecutivos =
  zipWith (-) (tail cubos) cubos 

cubos :: [Integer]
cubos = map (^3) [0..]

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,
-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]
primosConCubos4 =
  [p | x <- [1..],
       let p = 3*x^2+ 3*x + 1,
       isPrime p]

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

-- La propiedad es
prop_primosConCubos :: NonNegative Int -> Bool
prop_primosConCubos (NonNegative n) =
  all (== primosConCubos1 !! n)
      [primosConCubos2 !! n,
       primosConCubos3 !! n,
       primosConCubos4 !! n]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_primosConCubos
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> primosConCubos1 !! 6
--    331
--    (1.15 secs, 1,105,612,336 bytes)
--    λ> primosConCubos1 !! 7
--    397
--    (1.96 secs, 1,909,369,592 bytes)
--    λ> primosConCubos2 !! 7
--    397
--    
--    (0.01 secs, 648,840 bytes)
--    λ> primosConCubos2 !! (10^3)
--    65580901
--    (0.53 secs, 1,726,837,688 bytes)
--    λ> primosConCubos3 !! (10^3)
--    65580901
--    (0.49 secs, 1,724,258,632 bytes)
--    λ> primosConCubos4 !! (10^3)
--    65580901
--    (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura
-- ============================

-- Vamos a demostrar que los primos con cubos son diferencias de dos
-- cubos consecutivos.
--
-- Sea p un primo con cubo. Por la definición, existe un entero
-- positivo n tal que n³ + n²p es un cubo. 
--
-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).
-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que 
-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,
-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de
-- a³ es a³+3a²+3a+1.
--
-- Lema 2: Los números n² y n+p son coprimos.
-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;
-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se
-- tiene que k = 1.
--
-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene
-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es
-- decir, existen enteros positivos x e y tales que n = x³ y 
-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya
-- que 
--    p = y³-x³
--      = (n+k)³-n³
--      = 3k+3k²+k³
-- no es primo para k > 1.
--
-- Por consiguiente, p = (x+1)³-x³.

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

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs)

-- 1ª solución

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

primosConCubos1 :: [Integer]

primosConCubos1 =

[p | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

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

isPrime p]

-- 2ª solución

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

-- Para analizar la respuesta, en esta solución se calculan los pares

-- (p,n) tales que p es un primo con cubo y n es un número positivo tal

-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]

primosConCubos2' =

[(p,n) | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

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

isPrime p]

-- El cálculo es

-- λ> take 7 primosConCubos2

-- [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]

-- es la de los cubos y que los primeros elementos se obtienen restando

-- los segundos elementos consecutivos; es decir,

-- 7 = 8 - 1 = 2^3 - 1^3

-- 19 = 27 - 8 = 3^3 - 2^3

-- 37 = 64 - 27 = 4^3 - 3^3

-- Continuando el patrón,

-- 61 = 5^3 - 4^3 = 125 - 64

-- 91 = 6^3 - 5^3 = 216 - 125

-- 127 = 7^3 - 6^3 = 343 - 216

-- 271 = 10^3 - 9^3 = 1000 - 729

-- 331 = 11^3 -10^3 = 1331 - 1000

-- Por tanto, los primos con cubos son diferencias de dos cubos

-- consecutivos; es decir, coinciden con los números cubanos del

-- ejercicio anterior. A partir de la conjetura se obtienen las

-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente

-- definición

primosConCubos2 :: [Integer]

primosConCubos2 =

filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª definición

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

primosConCubos3 :: [Integer]

primosConCubos3 =

filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]

diferenciasCubosConsecutivos =

zipWith (-) (tail cubos) cubos

cubos :: [Integer]

cubos = map (^3) [0..]

-- 4ª solución

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,

-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]

primosConCubos4 =

[p | x <- [1..],

let p = 3*x^2+ 3*x + 1,

isPrime p]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosConCubos :: NonNegative Int -> Bool

prop_primosConCubos (NonNegative n) =

all (== primosConCubos1 !! n)

[primosConCubos2 !! n,

primosConCubos3 !! n,

primosConCubos4 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosConCubos1 !! 6

-- 331

-- (1.15 secs, 1,105,612,336 bytes)

-- λ> primosConCubos1 !! 7

-- 397

-- (1.96 secs, 1,909,369,592 bytes)

-- λ> primosConCubos2 !! 7

-- 397

-- (0.01 secs, 648,840 bytes)

-- λ> primosConCubos2 !! (10^3)

-- 65580901

-- (0.53 secs, 1,726,837,688 bytes)

-- λ> primosConCubos3 !! (10^3)

-- 65580901

-- (0.49 secs, 1,724,258,632 bytes)

-- λ> primosConCubos4 !! (10^3)

-- 65580901

-- (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura

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

-- Vamos a demostrar que los primos con cubos son diferencias de dos

-- cubos consecutivos.

-- Sea p un primo con cubo. Por la definición, existe un entero

-- positivo n tal que n³ + n²p es un cubo.

-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).

-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que

-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,

-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de

-- a³ es a³+3a²+3a+1.

-- Lema 2: Los números n² y n+p son coprimos.

-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;

-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se

-- tiene que k = 1.

-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene

-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es

-- decir, existen enteros positivos x e y tales que n = x³ y

-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya

-- que

-- p = y³-x³

-- = (n+k)³-n³

-- = 3k+3k²+k³

-- no es primo para k > 1.

-- Por consiguiente, p = (x+1)³-x³.

El código se encuentra en GitHub.

Primos cubanos

PorJosé A. Alonso 15-julio-20229-julio-2022

Un primo cubano es un número primo que se puede escribir como diferencia de dos cubos consecutivos. Por ejemplo, el 61 es un primo cubano porque es primo y 61 = 5³-4³.

Definir la sucesión

   cubanos :: [Integer]

1	cubanos :: [Integer]

tal que sus elementos son los números cubanos. Por ejemplo,

   λ> take 15 cubanos
   [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

1 2	λ> take 15 cubanos [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

Soluciones

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

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

cubanos1 :: [Integer]
cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos) 

-- cubos es la lista de los cubos. Por ejemplo,
--    λ> take 10 cubos
--    [1,8,27,64,125,216,343,512,729,1000]
cubos :: [Integer]
cubos = map (^3) [1..]

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

cubanos2 :: [Integer]
cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

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

cubanos3 :: [Integer]
cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

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

-- La propiedad es
prop_cubanos :: NonNegative Int -> Bool
prop_cubanos (NonNegative n) =
  all (== cubanos1 !! n)
      [cubanos2 !! n,
       cubanos3 !! n]

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

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

-- La comparación es
--    λ> cubanos1 !! 3000
--    795066361
--    (4.21 secs, 16,953,612,192 bytes)
--    λ> cubanos2 !! 3000
--    795066361
--    (4.27 secs, 16,962,597,288 bytes)
--    λ> cubanos3 !! 3000
--    795066361
--    (4.29 secs, 16,956,085,672 bytes)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª solución

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

cubanos1 :: [Integer]

cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos)

-- cubos es la lista de los cubos. Por ejemplo,

-- λ> take 10 cubos

-- [1,8,27,64,125,216,343,512,729,1000]

cubos :: [Integer]

cubos = map (^3) [1..]

-- 2ª solución

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

cubanos2 :: [Integer]

cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª solución

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

cubanos3 :: [Integer]

cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_cubanos :: NonNegative Int -> Bool

prop_cubanos (NonNegative n) =

all (== cubanos1 !! n)

[cubanos2 !! n,

cubanos3 !! n]

-- La comprobación es

-- λ> quickCheck prop_cubanos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> cubanos1 !! 3000

-- 795066361

-- (4.21 secs, 16,953,612,192 bytes)

-- λ> cubanos2 !! 3000

-- 795066361

-- (4.27 secs, 16,962,597,288 bytes)

-- λ> cubanos3 !! 3000

-- 795066361

-- (4.29 secs, 16,956,085,672 bytes)

El código se encuentra en GitHub.

Mayor semiprimo menor que n

PorJosé A. Alonso 7-julio-20226-julio-2022

Un número semiprimo es un número natural es producto de dos números primos no necesariamente distintos. Por ejemplo, 26 es semiprimo (porque 26 = 2·13) y 49 también lo es (porque 49 = 7·7).

Definir la función

   mayorSemiprimoMenor :: Integer -> Integer

1	mayorSemiprimoMenor :: Integer -> Integer

tal que (mayorSemiprimoMenor n) es el mayor semiprimo menor que n (suponiendo que n > 4). Por ejemplo,

   mayorSemiprimoMenor 27      ==  26
   mayorSemiprimoMenor 50      ==  49
   mayorSemiprimoMenor 49      ==  46
   mayorSemiprimoMenor (10^15) == 999999999999998

mayorSemiprimoMenor 27 == 26

mayorSemiprimoMenor 50 == 49

mayorSemiprimoMenor 49 == 46

mayorSemiprimoMenor (10^15) == 999999999999998

Soluciones

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

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

mayorSemiprimoMenor1 :: Integer -> Integer
mayorSemiprimoMenor1 n =
  head [x | x <- [n-1,n-2..2], semiPrimo x]

semiPrimo :: Integer -> Bool
semiPrimo n =
  not (null [x | x <- [n,n-1..2], 
                 primo x,
                 n `mod` x == 0,
                 primo (n `div` x)])

primo :: Integer -> Bool
primo n = [x | x <- [1..n], n `mod` x == 0] == [1,n] 

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

mayorSemiprimoMenor2 :: Integer -> Integer
mayorSemiprimoMenor2 n =
  head [x | x <- [n-1,n-2..2], semiPrimo2 x]

semiPrimo2 :: Integer -> Bool
semiPrimo2 n =
  not (null [x | x <- [n-1,n-2..2], 
                 isPrime x,
                 n `mod` x == 0,
                 isPrime (n `div` x)])

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

mayorSemiprimoMenor3 :: Integer -> Integer
mayorSemiprimoMenor3 n =
  head [x | x <- [n-1,n-2..2], semiPrimo3 x]

semiPrimo3 :: Integer -> Bool
semiPrimo3 n =
  not (null [x | x <- reverse (takeWhile (<n) primes),
                 n `mod` x == 0,
                 isPrime (n `div` x)])

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

mayorSemiprimoMenor4 :: Integer -> Integer
mayorSemiprimoMenor4 n =
  head [ p | p <- [n-1,n-2..2]
           , (length . primeFactors) p == 2]

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

mayorSemiprimoMenor5 :: Integer -> Integer
mayorSemiprimoMenor5 n
  | semiPrimo5 (n-1) = n-1
  | otherwise        = mayorSemiprimoMenor5 (n-1)

semiPrimo5 :: Integer -> Bool
semiPrimo5 x = length (primeFactors x) == 2

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

-- La propiedad es
prop_mayorSemiprimoMenor :: Integer -> Property
prop_mayorSemiprimoMenor n =
  n > 4 ==>
  all (== mayorSemiprimoMenor1 n)
      [mayorSemiprimoMenor2 n,
       mayorSemiprimoMenor3 n,
       mayorSemiprimoMenor4 n,
       mayorSemiprimoMenor5 n]

-- La comprobación es
--    λ> quickCheck prop_mayorSemiprimoMenor
--    +++ OK, passed 100 tests; 353 discarded.

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

-- La comparación es
--    λ> mayorSemiprimoMenor1 5000
--    4997
--    (1.92 secs, 945,507,880 bytes)
--    λ> mayorSemiprimoMenor2 5000
--    4997
--    (0.05 secs, 123,031,264 bytes)
--    λ> mayorSemiprimoMenor3 5000
--    4997
--    (0.01 secs, 5,865,120 bytes)
--    λ> mayorSemiprimoMenor4 5000
--    4997
--    (0.00 secs, 593,528 bytes)
--    λ> mayorSemiprimoMenor5 5000
--    4997
--    (0.00 secs, 593,200 bytes)
--
--    λ> mayorSemiprimoMenor3 (3*10^6)
--    2999995
--    (2.34 secs, 6,713,620,000 bytes)
--    λ> mayorSemiprimoMenor4 (2*10^6)
--    1999997
--    (0.01 secs, 728,936 bytes)
--    λ> mayorSemiprimoMenor5 (2*10^6)
--    1999997
--    (0.01 secs, 728,608 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

import Data.Numbers.Primes (primeFactors, isPrime, primes)

import Test.QuickCheck

-- 1ª solución

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

mayorSemiprimoMenor1 :: Integer -> Integer

mayorSemiprimoMenor1 n =

head [x | x <- [n-1,n-2..2], semiPrimo x]

semiPrimo :: Integer -> Bool

semiPrimo n =

not (null [x | x <- [n,n-1..2],

primo x,

n `mod` x == 0,

primo (n `div` x)])

primo :: Integer -> Bool

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

-- 2ª solución

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

mayorSemiprimoMenor2 :: Integer -> Integer

mayorSemiprimoMenor2 n =

head [x | x <- [n-1,n-2..2], semiPrimo2 x]

semiPrimo2 :: Integer -> Bool

semiPrimo2 n =

not (null [x | x <- [n-1,n-2..2],

isPrime x,

n `mod` x == 0,

isPrime (n `div` x)])

-- 3ª solución

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

mayorSemiprimoMenor3 :: Integer -> Integer

mayorSemiprimoMenor3 n =

head [x | x <- [n-1,n-2..2], semiPrimo3 x]

semiPrimo3 :: Integer -> Bool

semiPrimo3 n =

not (null [x | x <- reverse (takeWhile (<n) primes),

n `mod` x == 0,

isPrime (n `div` x)])

-- 4ª solución

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

mayorSemiprimoMenor4 :: Integer -> Integer

mayorSemiprimoMenor4 n =

head [ p | p <- [n-1,n-2..2]

, (length . primeFactors) p == 2]

-- 5ª solución

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

mayorSemiprimoMenor5 :: Integer -> Integer

mayorSemiprimoMenor5 n

| semiPrimo5 (n-1) = n-1

| otherwise = mayorSemiprimoMenor5 (n-1)

semiPrimo5 :: Integer -> Bool

semiPrimo5 x = length (primeFactors x) == 2

-- Comprobación de equivalencia

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

-- La propiedad es

prop_mayorSemiprimoMenor :: Integer -> Property

prop_mayorSemiprimoMenor n =

n > 4 ==>

all (== mayorSemiprimoMenor1 n)

[mayorSemiprimoMenor2 n,

mayorSemiprimoMenor3 n,

mayorSemiprimoMenor4 n,

mayorSemiprimoMenor5 n]

-- La comprobación es

-- λ> quickCheck prop_mayorSemiprimoMenor

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

-- Comparación de eficiencia

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

-- La comparación es

-- λ> mayorSemiprimoMenor1 5000

-- 4997

-- (1.92 secs, 945,507,880 bytes)

-- λ> mayorSemiprimoMenor2 5000

-- 4997

-- (0.05 secs, 123,031,264 bytes)

-- λ> mayorSemiprimoMenor3 5000

-- 4997

-- (0.01 secs, 5,865,120 bytes)

-- λ> mayorSemiprimoMenor4 5000

-- 4997

-- (0.00 secs, 593,528 bytes)

-- λ> mayorSemiprimoMenor5 5000

-- 4997

-- (0.00 secs, 593,200 bytes)

-- λ> mayorSemiprimoMenor3 (3*10^6)

-- 2999995

-- (2.34 secs, 6,713,620,000 bytes)

-- λ> mayorSemiprimoMenor4 (2*10^6)

-- 1999997

-- (0.01 secs, 728,936 bytes)

-- λ> mayorSemiprimoMenor5 (2*10^6)

-- 1999997

-- (0.01 secs, 728,608 bytes)

El código se encuentra en GitHub.

Ceros finales del factorial

PorJosé A. Alonso 4-julio-20226-julio-2022

Definir la función

   cerosDelFactorial :: Integer -> Integer

1	cerosDelFactorial :: Integer -> Integer

tal que (cerosDelFactorial n) es el número de ceros en que termina el factorial de n. Por ejemplo,

   cerosDelFactorial 24                         == 4
   cerosDelFactorial 25                         == 6
   length (show (cerosDelFactorial (10^70000))) == 70000

cerosDelFactorial 24 == 4

cerosDelFactorial 25 == 6

length (show (cerosDelFactorial (10^70000))) == 70000

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

cerosDelFactorial1 :: Integer -> Integer
cerosDelFactorial1 n = ceros (factorial n)

-- (factorial n) es el factorial n. Por ejemplo,
--    factorial 3  ==  6
factorial :: Integer -> Integer
factorial n = product [1..n]

-- (ceros n) es el número de ceros en los que termina el número n. Por
-- ejemplo, 
--    ceros 320000  ==  4
ceros :: Integer -> Integer
ceros n | rem n 10 /= 0 = 0
        | otherwise     = 1 + ceros (div n 10)

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

cerosDelFactorial2 :: Integer -> Integer
cerosDelFactorial2 = ceros2 . factorial 

ceros2 :: Integer -> Integer
ceros2 n = genericLength (takeWhile (=='0') (reverse (show n)))

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

cerosDelFactorial3 :: Integer -> Integer
cerosDelFactorial3 n
  | n < 5     = 0
  | otherwise = m + cerosDelFactorial3 m
  where m = n `div` 5

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

-- La propiedad es
prop_cerosDelFactorial :: Positive Integer -> Bool
prop_cerosDelFactorial (Positive n) =
  all (== cerosDelFactorial1 n)
      [cerosDelFactorial2 n,
       cerosDelFactorial3 n]

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

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

-- La comparación es
--    λ> cerosDelFactorial1 (4*10^4)
--    9998
--    (1.93 secs, 2,296,317,904 bytes)
--    λ> cerosDelFactorial2 (4*10^4)
--    9998
--    (1.57 secs, 1,636,242,040 bytes)
--    λ> cerosDelFactorial3 (4*10^4)
--    9998
--    (0.02 secs, 527,584 bytes)

import Data.List (genericLength)

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

-- 1ª solución

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

cerosDelFactorial1 :: Integer -> Integer

cerosDelFactorial1 n = ceros (factorial n)

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

-- factorial 3 == 6

factorial :: Integer -> Integer

factorial n = product [1..n]

-- (ceros n) es el número de ceros en los que termina el número n. Por

-- ejemplo,

-- ceros 320000 == 4

ceros :: Integer -> Integer

ceros n | rem n 10 /= 0 = 0

| otherwise = 1 + ceros (div n 10)

-- 2ª solución

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

cerosDelFactorial2 :: Integer -> Integer

cerosDelFactorial2 = ceros2 . factorial

ceros2 :: Integer -> Integer

ceros2 n = genericLength (takeWhile (=='0') (reverse (show n)))

-- 3ª solución

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

cerosDelFactorial3 :: Integer -> Integer

cerosDelFactorial3 n

| n < 5 = 0

| otherwise = m + cerosDelFactorial3 m

where m = n `div` 5

-- Comprobación de equivalencia

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

-- La propiedad es

prop_cerosDelFactorial :: Positive Integer -> Bool

prop_cerosDelFactorial (Positive n) =

all (== cerosDelFactorial1 n)

[cerosDelFactorial2 n,

cerosDelFactorial3 n]

-- La comprobación es

-- λ> quickCheck prop_cerosDelFactorial

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> cerosDelFactorial1 (4*10^4)

-- 9998

-- (1.93 secs, 2,296,317,904 bytes)

-- λ> cerosDelFactorial2 (4*10^4)

-- 9998

-- (1.57 secs, 1,636,242,040 bytes)

-- λ> cerosDelFactorial3 (4*10^4)

-- 9998

-- (0.02 secs, 527,584 bytes)

El código se encuentra en GitHub.

Números autodescriptivos

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

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

Definir la función

   autodescriptivo :: Integer -> Bool

1	autodescriptivo :: Integer -> Bool

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

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

λ> autodescriptivo 1210

True

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

[1210,2020,21200]

λ> autodescriptivo 9210000001000

True

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

import Data.Char (digitToInt)

import Data.List (genericLength)

import Test.QuickCheck

-- 1ª solución

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

autodescriptivo1 :: Integer -> Bool

autodescriptivo1 n = autodescriptiva (digitos n)

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

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

digitos :: Integer -> [Integer]

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

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

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

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

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

-- ejemplo,

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

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

autodescriptiva :: [Integer] -> Bool

autodescriptiva ns =

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

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

-- ejemplo,

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

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

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

-- 2ª solución

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

autodescriptivo2 :: Integer -> Bool

autodescriptivo2 n =

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

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

where xs = show n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_autodescriptivo :: Positive Integer -> Bool

prop_autodescriptivo (Positive n) =

autodescriptivo1 n == autodescriptivo2 n

-- La comprobación es

-- λ> quickCheck prop_autodescriptivo

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

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

-- [1210,2020,21200]

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

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

-- [1210,2020,21200]

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

El código se encuentra en GitHub.

Pandigitales primos

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

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

Definir la lista

   pandigitalesPrimos :: [Int]

1	pandigitalesPrimos :: [Int]

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

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

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

2143 `elem` pandigitalesPrimos == True

length pandigitalesPrimos == 538

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

import Data.List (permutations, sort)

import Data.Char (intToDigit)

import Data.Numbers.Primes (isPrime)

-- 1ª solución

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

pandigitalesPrimos1 :: [Int]

pandigitalesPrimos1 =

concatMap nPandigitalesPrimos1 [9,8..1]

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

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

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

-- nPandigitalesPrimos 5 == []

nPandigitalesPrimos1 :: Int -> [Int]

nPandigitalesPrimos1 n = filter isPrime (pandigitales n)

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

-- ordenada de mayor a menor. Por ejemplo,

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

pandigitales :: Int -> [Int]

pandigitales n =

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

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

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

digitosAentero :: [Int] -> Int

digitosAentero = read . map intToDigit

-- 2ª solución

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

pandigitalesPrimos2 :: [Int]

pandigitalesPrimos2 =

concatMap nPandigitalesPrimos2 [9,8..1]

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

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

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

-- son primos.

nPandigitalesPrimos2 :: Int -> [Int]

nPandigitalesPrimos2 n

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

| otherwise = filter isPrime (pandigitales n)

-- 2ª solución

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

pandigitalesPrimos3 :: [Int]

pandigitalesPrimos3 =

concatMap nPandigitalesPrimos3 [9,8..1]

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

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

-- [1,4,7]

nPandigitalesPrimos3 :: Int -> [Int]

nPandigitalesPrimos3 n

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

| otherwise = []

-- Comprobación de equivalencia

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

-- La propiedad es

prop_pandigitalesPrimos :: Bool

prop_pandigitalesPrimos =

all (== pandigitalesPrimos1)

[pandigitalesPrimos2,

pandigitalesPrimos3]

-- La comprobación es

-- λ> prop_pandigitalesPrimos

-- True

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (pandigitalesPrimos1)

-- 538

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

-- λ> length (pandigitalesPrimos2)

-- 538

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

-- λ> length (pandigitalesPrimos3)

-- 538

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

El código se encuentra en GitHub.

Codificación de Fibonacci

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

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

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

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

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

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

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

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

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

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

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

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

1 = 1 = F(2) ≡ 11

2 = 2 = F(3) ≡ 011

3 = 3 = F(4) ≡ 0011

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

5 = 5 = F(5) ≡ 00011

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

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

8 = 8 = F(6) ≡ 000011

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

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

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

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

13 = 13 = F(7) ≡ 0000011

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

Definir la función

   codigoFib :: Integer -> String

1	codigoFib :: Integer -> String

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

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

λ> codigoFib 65

"0100100011"

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

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

Comprobar con QuickCheck las siguientes propiedades:

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

import Data.List (isInfixOf)

import Data.Array (Array, accumArray, elems)

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

-- 1ª solución

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

codigoFib1 :: Integer -> String

codigoFib1 = concatMap show . codificaFibLista

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

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

-- λ> codificaFibLista 65

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

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

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

codificaFibLista :: Integer -> [Integer]

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

where xs = map fst (descomposicion n)

f i | i `elem` xs = 1

| otherwise = 0

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

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

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

-- ejemplo,

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

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

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

descomposicion 0 = []

descomposicion 1 = [(2,1)]

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

where (i,x) = fibAnterior n

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

-- n. Por ejemplo,

-- fibAnterior 33 == (8,21)

-- fibAnterior 34 == (9,34)

fibAnterior :: Integer -> (Integer, Integer)

fibAnterior n = last (takeWhile p fibsConIndice)

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

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

-- sus índices. Por ejemplo,

-- λ> take 10 fibsConIndice

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

fibsConIndice :: [(Integer, Integer)]

fibsConIndice = zip [0..] fibs

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

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

fibs :: [Integer]

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

--- 2ª solución

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

codigoFib2 :: Integer -> String

codigoFib2 = concatMap show . elems . codificaFibVec

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

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

-- λ> codificaFibVec 65

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

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

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

codificaFibVec :: Integer -> Array Integer Integer

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

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

a = fst (head is)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_codigoFib :: Positive Integer -> Bool

prop_codigoFib (Positive n) =

codigoFib1 n == codigoFib2 n

-- La comprobación es

-- λ> quickCheck prop_codigoFib

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

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

-- 121393

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

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

-- 121393

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

-- Propiedades

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

-- Usaremos la 2ª definición

codigoFib :: Integer -> String

codigoFib = codigoFib2

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

prop_descomposicion_correcta :: Positive Integer -> Bool

prop_descomposicion_correcta (Positive n) =

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

-- La comprobación es

-- λ> quickCheck prop_descomposicion_correcta

-- +++ OK, passed 100 tests.

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

-- la sucesión de Fibonacci.

prop_descomposicion :: Positive Integer -> Bool

prop_descomposicion (Positive n) =

not (null (descomposicion n))

-- La comprobación es

-- λ> quickCheck prop_descomposicion

-- +++ OK, passed 100 tests.

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

prop_length_codigoFib :: Positive Integer -> Bool

prop_length_codigoFib (Positive n) =

length (codigoFib n) >= 2

-- La comprobación es

-- λ> quickCheck prop_length_codigoFib

-- +++ OK, passed 100 tests.

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

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

prop3_cadena_11_en_codigoFib :: Positive Integer -> Bool

prop3_cadena_11_en_codigoFib (Positive n) =

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

where xs = reverse (codigoFib n)

-- La comprobación es

-- λ> quickCheck prop3_cadena_11_en_codigoFib

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

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

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

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

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

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

Definir la sucesión

   especiales :: [Integer]

1	especiales :: [Integer]

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

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

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

especiales !! 50 == 84

especiales !! 500 == 638

especiales !! 5000 == 5806

especiales !! 50000 == 55746

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- La comparación es

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

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

-- 1ª solución

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

especiales1 :: [Integer]

especiales1 = filter especial1 [1..]

especial1 :: Integer -> Bool

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

mcm1 :: [Integer] -> Integer

mcm1 [] = 1

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

-- 2ª solución

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

especiales2 :: [Integer]

especiales2 = filter especial2 [1..]

especial2 :: Integer -> Bool

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

mcm2 :: [Integer] -> Integer

mcm2 = foldr lcm 1

-- 3ª solución

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

especiales3 :: [Integer]

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

, x == y]

mcms :: [(Integer,Integer)]

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_especiales :: NonNegative Int -> Bool

prop_especiales (NonNegative n) =

all (== especiales1 !! n)

[especiales2 !! n,

especiales3 !! n]

-- La comprobación es

-- λ> quickCheck prop_especiales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- Comparación

-- λ> especiales1 !! 2000

-- 2390

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

-- λ> especiales2 !! 2000

-- 2390

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

-- λ> especiales3 !! 2000

-- 2390

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

El código se encuentra en GitHub.

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

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

Definir la función

   suma :: Integer -> Integer

1	suma :: Integer -> Integer

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

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

suma 1 == 1

suma 2 == 5

suma 3 == 23

suma 4 == 119

suma 5 == 719

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- 1ª solución

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

suma1 :: Integer -> Integer

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

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª solución

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

suma2 :: Integer -> Integer

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

factoriales :: [Integer]

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

-- 3ª solución

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

-- Basada en los siguientes cálculos

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

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

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

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

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

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

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

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

suma3 :: Integer -> Integer

suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_suma :: Positive Integer -> Bool

prop_suma (Positive n) =

all (== suma1 n)

[suma2 n,

suma3 n]

-- La comprobación es

-- λ> quickCheck prop_suma

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

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

-- "73170"

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

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

-- "73170"

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

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

-- "73170"

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

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

-- "83669"

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

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

-- "83669"

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

El código se encuentra en GitHub.

Menor número con una cantidad dada de divisores

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

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

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

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

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

Definir la función

   menor :: Integer -> Integer

1	menor :: Integer -> Integer

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

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

menor 1 == 2

menor 2 == 6

menor 3 == 24

menor 4 == 120

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

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

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

100

101

102

103

import Data.List (genericLength, genericTake, group)

import Data.Numbers.Primes (primeFactors, primes)

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

quickCheck, quickCheckWith)

-- 1ª solución

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

menor1 :: Integer -> Integer

menor1 n =

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

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

-- numeroDivisores 12 == 6

numeroDivisores :: Integer -> Integer

numeroDivisores =

genericLength . divisores

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

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

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

divisores :: Integer -> [Integer]

divisores n =

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

-- 2ª solución

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

menor2 :: Integer -> Integer

menor2 n =

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

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

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

-- 3ª solución

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

menor3 :: Integer -> Integer

menor3 n = product (genericTake n potencias)

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

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

-- menor a mayor. Por ejemplo,

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

potencias :: [Integer]

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

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

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

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

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

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

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

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_menor :: Positive Integer -> Bool

prop_menor (Positive n) =

all (== menor1 n)

[menor2 n,

menor3 n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> menor 8

-- 1081080

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

-- λ> menor2 8

-- 1081080

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

-- λ> menor3 8

-- 1081080

-- (0.00 secs, 116,960 bytes)

-- Definición de menor

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

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

menor :: Integer -> Integer

menor = menor3

-- Propiedad

-- =========

-- La propiedad es

prop_menor_divide :: Positive Integer -> Bool

prop_menor_divide (Positive n) =

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

-- La comprobación es

-- λ> quickCheck prop_menor_divide

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos circulares

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

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

Definir la lista

   circulares :: [Integer]

1	circulares :: [Integer]

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

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

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

100

101

102

103

104

105

106

107

import Data.Numbers.Primes (isPrime, primes)

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

-- 1ª solución

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

circulares1 :: [Integer]

circulares1 = filter esCircular1 primes

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

-- esCircular1 197 == True

-- esCircular1 157 == False

esCircular1 :: Integer -> Bool

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

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

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

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

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

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

where aux 1 yss = yss

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

-- 2ª solución

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

circulares2 :: [Integer]

circulares2 = filter esCircular2 primes

esCircular2 :: Integer -> Bool

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

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

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

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

-- final. Por ejemplo,

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

rota :: [a] -> [a]

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

-- 3ª solución

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

circulares3 :: [Integer]

circulares3 = filter (all isPrime . rotaciones3) primes

rotaciones3 :: Integer -> [Integer]

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

where s = show n

m = length s

-- 4ª definición

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

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

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

-- impares, salvo para n = 2.

circulares4 :: [Integer]

circulares4 = 2 : filter esCircular4 primes

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

-- esCircular4 197 == True

-- esCircular4 157 == False

esCircular4 :: Integer -> Bool

esCircular4 n = digitosImpares n &&

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

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

-- impares. Por ejemplo,

-- digitosImpares 7351 == True

-- digitosImpares 7341 == False

digitosImpares :: Integer -> Bool

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_circulares :: Positive Int -> Bool

prop_circulares (Positive n) =

all (== circulares1 !! n)

[circulares2 !! n,

circulares3 !! n,

circulares4 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> circulares1 !! 46

-- 331999

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

-- λ> circulares2 !! 46

-- 331999

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

-- λ> circulares3 !! 46

-- 331999

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

-- λ> circulares4 !! 46

-- 331999

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

El código se encuentra en GitHub.

Término ausente en una progresión aritmética

PorJosé A. Alonso 1-junio-202230-mayo-2022

Una progresión aritmética es una sucesión de números tales que la diferencia de dos términos sucesivos cualesquiera de la sucesión es constante.

Definir la función

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

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

tal que (ausente xs) es el único término ausente de la progresión aritmética xs. Por ejemplo,

   ausente [3,7,9,11]               ==  5
   ausente [3,5,9,11]               ==  7
   ausente [3,5,7,11]               ==  9
   ausente ([1..9]++[11..])         ==  10
   ausente ([1..10^6] ++ [2+10^6])  ==  1000001

ausente [3,7,9,11] == 5

ausente [3,5,9,11] == 7

ausente [3,5,7,11] == 9

ausente ([1..9]++[11..]) == 10

ausente ([1..10^6] ++ [2+10^6]) == 1000001

Nota. Se supone que la lista tiene al menos 3 elementos, que puede ser infinita y que sólo hay un término de la progresión aritmética que no está en la lista.

Soluciones

import Data.List (group, genericLength)
import Test.QuickCheck (Arbitrary, Gen,
                        arbitrary, frequency, suchThat, quickCheck) 

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

ausente1 :: Integral a => [a] -> a
ausente1 (x1:xs@(x2:x3:_))
  | d1 == d2     = ausente1 xs
  | d1 == 2 * d2 = x1 + d2
  | d2 == 2 * d1 = x2 + d1
  where d1 = x2 - x1
        d2 = x3 - x2          

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

ausente2 :: Integral a => [a] -> a
ausente2 s@(x1:x2:x3:_) 
  | x1 + x3 /= 2 * x2 = x1 + (x3 - x2)
  | otherwise         = head [a | (a,b) <- zip [x1,x2..] s
                                , a /= b]

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

ausente3 :: Integral a => [a] -> a
ausente3  xs@(x1:x2:_) 
  | null us   = x1 + v
  | otherwise = x2 + u * genericLength (u:us) 
  where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

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

-- Tipo de progresiones aritméticas con un término ausente.
newtype PAconAusente = PA [Integer]
  deriving Show

-- Generación de progresiones aritméticas con un término ausente. 
progresionConAusenteArbitraria :: Gen PAconAusente
progresionConAusenteArbitraria = do
  x <- arbitrary
  d <- arbitrary `suchThat` (/= 0)
  n <- arbitrary `suchThat` (> 2)
  k <- arbitrary `suchThat` (> 0)
  let (as,_:bs) = splitAt k [x,x+d..]
  frequency [(1,return (PA (as ++ bs))),
             (1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.
instance Arbitrary PAconAusente where
  arbitrary = progresionConAusenteArbitraria

-- La propiedad es
prop_ausente :: PAconAusente -> Bool 
prop_ausente (PA xs) =
  all (== ausente1 xs)
      [ausente2 xs,
       ausente3 xs]

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

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

-- La comparación es
--    λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])
--    1000001
--    (1.15 secs, 560,529,520 bytes)
--    λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])
--    1000001
--    (0.33 secs, 336,530,680 bytes)
--    λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])
--    1000001
--    (0.50 secs, 498,047,584 bytes)

import Data.List (group, genericLength)

import Test.QuickCheck (Arbitrary, Gen,

arbitrary, frequency, suchThat, quickCheck)

-- 1ª solución

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

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

ausente1 (x1:xs@(x2:x3:_))

| d1 == d2 = ausente1 xs

| d1 == 2 * d2 = x1 + d2

| d2 == 2 * d1 = x2 + d1

where d1 = x2 - x1

d2 = x3 - x2

-- 2ª solución

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

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

ausente2 s@(x1:x2:x3:_)

| x1 + x3 /= 2 * x2 = x1 + (x3 - x2)

| otherwise = head [a | (a,b) <- zip [x1,x2..] s

, a /= b]

-- 3ª solución

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

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

ausente3 xs@(x1:x2:_)

| null us = x1 + v

| otherwise = x2 + u * genericLength (u:us)

where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

-- Comprobación de equivalencia

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

-- Tipo de progresiones aritméticas con un término ausente.

newtype PAconAusente = PA [Integer]

deriving Show

-- Generación de progresiones aritméticas con un término ausente.

progresionConAusenteArbitraria :: Gen PAconAusente

progresionConAusenteArbitraria = do

x <- arbitrary

d <- arbitrary `suchThat` (/= 0)

n <- arbitrary `suchThat` (> 2)

k <- arbitrary `suchThat` (> 0)

let (as,_:bs) = splitAt k [x,x+d..]

frequency [(1,return (PA (as ++ bs))),

(1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.

instance Arbitrary PAconAusente where

arbitrary = progresionConAusenteArbitraria

-- La propiedad es

prop_ausente :: PAconAusente -> Bool

prop_ausente (PA xs) =

all (== ausente1 xs)

[ausente2 xs,

ausente3 xs]

-- La comprobación es

-- λ> quickCheck prop_ausente

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])

-- 1000001

-- (1.15 secs, 560,529,520 bytes)

-- λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])

-- 1000001

-- (0.33 secs, 336,530,680 bytes)

-- λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])

-- 1000001

-- (0.50 secs, 498,047,584 bytes)

El código se encuentra en GitHub.

Ordenación de los racionales

PorJosé A. Alonso 30-mayo-202229-mayo-2022

En este ejercicio, representamos las fracciones mediante pares de números de enteros.

Definir la función

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

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

tal que (fraccionesOrd n) es la lista con las fracciones propias positivas ordenadas, con denominador menor o igual que n. Por ejemplo,

   λ> fraccionesOrd 4
   [(1,4),(1,3),(1,2),(2,3),(3,4)]
   λ> fraccionesOrd 5
   [(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

λ> fraccionesOrd 4

[(1,4),(1,3),(1,2),(2,3),(3,4)]

λ> fraccionesOrd 5

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

Soluciones

import Data.List (sort, sortBy)
import Data.Ratio ((%), numerator, denominator)
import Test.QuickCheck

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

fraccionesOrd1 :: Integer -> [(Integer,Integer)]
fraccionesOrd1 n = 
  [(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))
                              | y <- [2..n], 
                                x <- [1..y-1], 
                                gcd x y == 1]]

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

fraccionesOrd2 :: Integer -> [(Integer,Integer)]
fraccionesOrd2 n = 
  map snd (sort [(fromIntegral x/fromIntegral y,(x,y))
                 | y <- [2..n], 
                   x <- [1..y-1], 
                   gcd x y == 1])

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

fraccionesOrd3 :: Integer -> [(Integer,Integer)]
fraccionesOrd3 n = 
  sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]
  where comp (a,b) (c,d) = compare (a*d) (b*c)

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

fraccionesOrd4 :: Integer -> [(Integer,Integer)]
fraccionesOrd4 n = 
  [(numerator x, denominator x) | x <- racionalesOrd4 n]
  
-- (racionalesOrd4 n) es la lista con los racionales ordenados, con
-- denominador menor o igual que n. Por ejemplo,
--    λ> racionalesOrd4 4
--    [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]
racionalesOrd4 :: Integer -> [Rational]
racionalesOrd4 n =
  sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

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

-- La propiedad es
prop_fraccionesOrd :: Positive Integer -> Bool 
prop_fraccionesOrd (Positive n) =
  all (== fraccionesOrd1 n)
      [fraccionesOrd2 n,
       fraccionesOrd3 n,
       fraccionesOrd4 n]

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

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

-- La comparación es
--    λ> length (fraccionesOrd1 2000)
--    1216587
--    (3.65 secs, 2,879,842,368 bytes)
--    λ> length (fraccionesOrd2 2000)
--    1216587
--    (3.36 secs, 2,870,109,640 bytes)
--    λ> length (fraccionesOrd3 2000)
--    1216587
--    (8.83 secs, 5,700,519,584 bytes)
--    λ> length (fraccionesOrd4 2000)
--    1216587
--    (4.12 secs, 5,181,904,336 bytes)

import Data.List (sort, sortBy)

import Data.Ratio ((%), numerator, denominator)

import Test.QuickCheck

-- 1ª solución

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

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

fraccionesOrd1 n =

[(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

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

gcd x y == 1]]

-- 2ª solución

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

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

fraccionesOrd2 n =

map snd (sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

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

gcd x y == 1])

-- 3ª solución

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

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

fraccionesOrd3 n =

sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]

where comp (a,b) (c,d) = compare (a*d) (b*c)

-- 4ª solución

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

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

fraccionesOrd4 n =

[(numerator x, denominator x) | x <- racionalesOrd4 n]

-- (racionalesOrd4 n) es la lista con los racionales ordenados, con

-- denominador menor o igual que n. Por ejemplo,

-- λ> racionalesOrd4 4

-- [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]

racionalesOrd4 :: Integer -> [Rational]

racionalesOrd4 n =

sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_fraccionesOrd :: Positive Integer -> Bool

prop_fraccionesOrd (Positive n) =

all (== fraccionesOrd1 n)

[fraccionesOrd2 n,

fraccionesOrd3 n,

fraccionesOrd4 n]

-- La comprobación es

-- λ> quickCheck prop_fraccionesOrd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (fraccionesOrd1 2000)

-- 1216587

-- (3.65 secs, 2,879,842,368 bytes)

-- λ> length (fraccionesOrd2 2000)

-- 1216587

-- (3.36 secs, 2,870,109,640 bytes)

-- λ> length (fraccionesOrd3 2000)

-- 1216587

-- (8.83 secs, 5,700,519,584 bytes)

-- λ> length (fraccionesOrd4 2000)

-- 1216587

-- (4.12 secs, 5,181,904,336 bytes)

El código se encuentra en GitHub.

Polinomios cuadráticos generadores de primos

PorJosé A. Alonso 27-mayo-2022

En 1772, Euler publicó que el polinomio n² + n + 41 genera 40 números primos para todos los valores de n entre 0 y 39. Sin embargo, cuando n = 40, 40²+40+41 = 40(40+1)+41 es divisible por 41.

Usando ordenadores, se descubrió que el polinomio n² – 79n + 1601 genera 80 números primos para todos los valores de n entre 0 y 79.

Definir la función

   generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

1	generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

tal que (generadoresMaximales n) es el par (m,xs) donde

xs es la lista de pares (x,y) tales que n²+xn+y es uno de polinomios que genera un número máximo de números primos consecutivos a partir de cero entre todos los polinomios de la forma n²+an+b, con |a| ≤ n y |b| ≤ n y
m es dicho número máximo.

Por ejemplo,

   generadoresMaximales    4  ==  ( 3,[(-2,3),(-1,3),(3,3)])
   generadoresMaximales    6  ==  ( 5,[(-1,5),(5,5)])
   generadoresMaximales   41  ==  (41,[(-1,41)])
   generadoresMaximales   50  ==  (43,[(-5,47)])
   generadoresMaximales  100  ==  (48,[(-15,97)])
   generadoresMaximales  200  ==  (53,[(-25,197)])
   generadoresMaximales 1650  ==  (80,[(-79,1601)])

generadoresMaximales 4 == ( 3,[(-2,3),(-1,3),(3,3)])

generadoresMaximales 6 == ( 5,[(-1,5),(5,5)])

generadoresMaximales 41 == (41,[(-1,41)])

generadoresMaximales 50 == (43,[(-5,47)])

generadoresMaximales 100 == (48,[(-15,97)])

generadoresMaximales 200 == (53,[(-25,197)])

generadoresMaximales 1650 == (80,[(-79,1601)])

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (primes, isPrime)
import I1M.PolOperaciones (valor, consPol, polCero)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales1 n =
  (m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])
  where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos (-1) 41     ==  41
--    nPrimos (-79) 1601  ==  80
nPrimos :: Integer -> Integer -> Int
nPrimos a b =
  length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

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

-- Notas:
-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que
--    0²+a*0+b = b es primo.
-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para
--    n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales2 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    a <- [-n..n],
                                    isPrime(1+a+b)]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales3 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    p <- takeWhile (<=1+2*n) primes,
                                    let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)
-- ===========================================

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales4 n = (m,map snd zs)
  where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,
                                     p <- takeWhile (<=1+2*n) primes,
                                     let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos2 (-1) 41     ==  41
--    nPrimos2 (-79) 1601  ==  80
nPrimos2 :: Integer -> Integer -> Int
nPrimos2 a b =
  length (takeWhile isPrime [valor p n | n <- [0..]])
  where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

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

-- La propiedad es
prop_generadoresMaximales :: Positive Integer -> Bool
prop_generadoresMaximales (Positive n) =
  all (equivalentes (generadoresMaximales1 n'))
      [generadoresMaximales2 n',
       generadoresMaximales3 n',
       generadoresMaximales4 n']
  where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool
equivalentes (n,xs) (m,ys) =
  n == m && sort xs == sort ys

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

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

-- La comparación es
--    λ> generadoresMaximales1 300
--    (56,[(-31,281)])
--    (2.10 secs, 2,744,382,760 bytes)
--    λ> generadoresMaximales2 300
--    (56,[(-31,281)])
--    (0.17 secs, 382,103,656 bytes)
--    λ> generadoresMaximales3 300
--    (56,[(-31,281)])
--    (0.19 secs, 346,725,872 bytes)
--    λ> generadoresMaximales4 300
--    (56,[(-31,281)])
--    (0.20 secs, 388,509,808 bytes)

100

101

102

103

104

105

106

107

108

import Data.List (sort)

import Data.Numbers.Primes (primes, isPrime)

import I1M.PolOperaciones (valor, consPol, polCero)

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

-- 1ª solución

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales1 n =

(m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])

where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos (-1) 41 == 41

-- nPrimos (-79) 1601 == 80

nPrimos :: Integer -> Integer -> Int

nPrimos a b =

length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

-- 2ª solución

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

-- Notas:

-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que

-- 0²+a*0+b = b es primo.

-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para

-- n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales2 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

a <- [-n..n],

isPrime(1+a+b)]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 3ª solución

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales3 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)

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

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales4 n = (m,map snd zs)

where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos2 (-1) 41 == 41

-- nPrimos2 (-79) 1601 == 80

nPrimos2 :: Integer -> Integer -> Int

nPrimos2 a b =

length (takeWhile isPrime [valor p n | n <- [0..]])

where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_generadoresMaximales :: Positive Integer -> Bool

prop_generadoresMaximales (Positive n) =

all (equivalentes (generadoresMaximales1 n'))

[generadoresMaximales2 n',

generadoresMaximales3 n',

generadoresMaximales4 n']

where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool

equivalentes (n,xs) (m,ys) =

n == m && sort xs == sort ys

-- La comprobación es

-- λ> quickCheck prop_generadoresMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> generadoresMaximales1 300

-- (56,[(-31,281)])

-- (2.10 secs, 2,744,382,760 bytes)

-- λ> generadoresMaximales2 300

-- (56,[(-31,281)])

-- (0.17 secs, 382,103,656 bytes)

-- λ> generadoresMaximales3 300

-- (56,[(-31,281)])

-- (0.19 secs, 346,725,872 bytes)

-- λ> generadoresMaximales4 300

-- (56,[(-31,281)])

-- (0.20 secs, 388,509,808 bytes)

El código se encuentra en GitHub.

El triángulo de Floyd

PorJosé A. Alonso 26-mayo-2022

El triángulo de Floyd, llamado así en honor a Robert Floyd, es un triángulo rectángulo formado con números naturales. Para crear un triángulo de Floyd, se comienza con un 1 en la esquina superior izquierda, y se continúa escribiendo la secuencia de los números naturales de manera que cada línea contenga un número más que la anterior. Las 5 primeras líneas del triángulo de Floyd son

    1
    2   3
    4   5   6
    7   8   9  10
   11  12  13  14  15

2 3

4 5 6

7 8 9 10

11 12 13 14 15

Definir la función

   trianguloFloyd :: [[Integer]]

1	trianguloFloyd :: [[Integer]]

tal que trianguloFloyd es el triángulo de Floyd. Por ejemplo,

   λ> take 4 trianguloFloyd
   [[1],
    [2,3],
    [4,5,6],
    [7,8,9,10]]
  (trianguloFloyd !! (10^5)) !! 0  ==  5000050001
  (trianguloFloyd !! (10^6)) !! 0  ==  500000500001
  (trianguloFloyd !! (10^7)) !! 0  ==  50000005000001

λ> take 4 trianguloFloyd

[[1],

[2,3],

[4,5,6],

[7,8,9,10]]

(trianguloFloyd !! (10^5)) !! 0 == 5000050001

(trianguloFloyd !! (10^6)) !! 0 == 500000500001

(trianguloFloyd !! (10^7)) !! 0 == 50000005000001

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

trianguloFloyd1 :: [[Integer]]
trianguloFloyd1 = floyd 1 [1..]
  where floyd n xs = i : floyd (n+1) r
          where (i,r) = splitAt n xs

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

trianguloFloyd2 :: [[Integer]]
trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el
-- triángulo de Floyd. Por ejemplo,
--    siguienteF [2,3]    ==  [4,5,6]
--    siguienteF [4,5,6]  ==  [7,8,9,10]
siguienteF :: [Integer] -> [Integer]
siguienteF xs = [a..a+n]
    where a = 1 + last xs
          n = genericLength xs

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

trianguloFloyd3 :: [[Integer]]
trianguloFloyd3 =
  [[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

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

trianguloFloyd4 :: [[Integer]]
trianguloFloyd4 =
  scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

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

-- La propiedad es
prop_trianguloFloyd :: Positive Int -> Bool
prop_trianguloFloyd (Positive n) =
  all (== (trianguloFloyd1 !! n))
      [trianguloFloyd2 !! n,
       trianguloFloyd3 !! n,
       trianguloFloyd4 !! n]

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

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

-- La comparación es
--    λ> (trianguloFloyd1 !! 5000) !! 5000
--    12507501
--    (1.47 secs, 2,505,005,752 bytes)
--    λ> (trianguloFloyd2 !! 5000) !! 5000
--    12507501
--    (0.79 secs, 2,416,259,176 bytes)
--    λ> (trianguloFloyd3 !! 5000) !! 5000
--    12507501
--    (0.00 secs, 1,809,152 bytes)
--    λ> (trianguloFloyd4 !! 5000) !! 5000
--    12507501
--    (0.01 secs, 3,517,896 bytes)
--
--    λ> (trianguloFloyd3 !! (10^7)) !! 0
--    50000005000001
--    (2.45 secs, 1,656,534,080 bytes)
--    λ> (trianguloFloyd4 !! (10^7)) !! 0
--    50000005000001
--    (10.86 secs, 5,302,760,752 bytes)

import Data.List (genericLength)

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

-- 1ª solución

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

trianguloFloyd1 :: [[Integer]]

trianguloFloyd1 = floyd 1 [1..]

where floyd n xs = i : floyd (n+1) r

where (i,r) = splitAt n xs

-- 2ª solución

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

trianguloFloyd2 :: [[Integer]]

trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el

-- triángulo de Floyd. Por ejemplo,

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

-- siguienteF [4,5,6] == [7,8,9,10]

siguienteF :: [Integer] -> [Integer]

siguienteF xs = [a..a+n]

where a = 1 + last xs

n = genericLength xs

-- 3ª solución

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

trianguloFloyd3 :: [[Integer]]

trianguloFloyd3 =

[[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución

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

trianguloFloyd4 :: [[Integer]]

trianguloFloyd4 =

scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_trianguloFloyd :: Positive Int -> Bool

prop_trianguloFloyd (Positive n) =

all (== (trianguloFloyd1 !! n))

[trianguloFloyd2 !! n,

trianguloFloyd3 !! n,

trianguloFloyd4 !! n]

-- La comprobación es

-- λ> quickCheck prop_trianguloFloyd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (trianguloFloyd1 !! 5000) !! 5000

-- 12507501

-- (1.47 secs, 2,505,005,752 bytes)

-- λ> (trianguloFloyd2 !! 5000) !! 5000

-- 12507501

-- (0.79 secs, 2,416,259,176 bytes)

-- λ> (trianguloFloyd3 !! 5000) !! 5000

-- 12507501

-- (0.00 secs, 1,809,152 bytes)

-- λ> (trianguloFloyd4 !! 5000) !! 5000

-- 12507501

-- (0.01 secs, 3,517,896 bytes)

-- λ> (trianguloFloyd3 !! (10^7)) !! 0

-- 50000005000001

-- (2.45 secs, 1,656,534,080 bytes)

-- λ> (trianguloFloyd4 !! (10^7)) !! 0

-- 50000005000001

-- (10.86 secs, 5,302,760,752 bytes)

El código se encuentra en GitHub.

Numeración con base múltiple

PorJosé A. Alonso 25-mayo-202225-mayo-2022

Sea (b(i) | i ≥ 1) una sucesión infinita de números enteros mayores que 1. Entonces todo entero x mayor que cero se puede escribir de forma única como

   x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

1	x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

donde cada x(i) satisface la condición 0 ≤ x(i) < b(i+1). Se dice que [x(n),x(n-1),…,x(2),x(1),x(0)] es la representación de x en la base (b(i)). Por ejemplo, la representación de 377 en la base (2, 6, 8, …) es [7,5,0,1] ya que

   377 = 1 + 0*2 + 5*2*4 + 7*2*4*6

1	377 = 1 + 02 + 524 + 7246

y, además, 0 ≤ 1 < 2, 0 ≤ 0 < 4, 0 ≤ 5 < 6 y 0 ≤ 7 < 8.

Definir las funciones

   decimalAmultiple :: [Integer] -> Integer -> [Integer]
   multipleAdecimal :: [Integer] -> [Integer] -> Integer

1 2	decimalAmultiple :: [Integer] -> Integer -> [Integer] multipleAdecimal :: [Integer] -> [Integer] -> Integer

tales que

(decimalAmultiple bs x) es la representación del número x en la base bs. Por ejemplo,

     decimalAmultiple [2,4..] 377                      ==  [7,5,0,1]
     decimalAmultiple [2,5..] 377                      ==  [4,5,3,1]
     decimalAmultiple [2^n | n <- [1..]] 2015          ==  [1,15,3,3,1]
     decimalAmultiple (repeat 10) 2015                 ==  [2,0,1,5]

decimalAmultiple [2,4..] 377 == [7,5,0,1]

decimalAmultiple [2,5..] 377 == [4,5,3,1]

decimalAmultiple [2^n | n <- [1..]] 2015 == [1,15,3,3,1]

decimalAmultiple (repeat 10) 2015 == [2,0,1,5]

(multipleAdecimal bs cs) es el número decimal cuya representación en la base bs es cs. Por ejemplo,

     multipleAdecimal [2,4..] [7,5,0,1]                ==  377
     multipleAdecimal [2,5..] [4,5,3,1]                ==  377
     multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1]  ==  2015
     multipleAdecimal (repeat 10) [2,0,1,5]            ==  2015

multipleAdecimal [2,4..] [7,5,0,1] == 377

multipleAdecimal [2,5..] [4,5,3,1] == 377

multipleAdecimal [2^n | n <- [1..]] [1,15,3,3,1] == 2015

multipleAdecimal (repeat 10) [2,0,1,5] == 2015

Comprobar con QuickCheck que se verifican las siguientes propiedades

Para cualquier base bs y cualquier entero positivo n,

     multipleAdecimal bs (decimalAmultiple bs x) == x

1	multipleAdecimal bs (decimalAmultiple bs x) == x

Para cualquier base bs y cualquier entero positivo n, el coefiente i-ésimo de la representación múltiple de n en la base bs es un entero no negativo menos que el i-ésimo elemento de bs.

Soluciones

import Test.QuickCheck
import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple
-- ===============================

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]
decimalAmultiple1 bs n = reverse (aux bs n)
  where aux _ 0      = []
        aux (d:ds) m = r : aux ds q
          where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple
-- ===============================

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]
decimalAmultiple2 bs n = aux bs n []
  where aux _ 0  xs     = xs
        aux (d:ds) m xs = aux ds q (r:xs)
          where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple
-- ===============================

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]
decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))
  where f (_     ,0) = Nothing
        f ((y:ys),m) = Just (r,(ys,q))
                       where (q,r) = quotRem m y

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

-- La propiedad es
prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_decimalAmultiple (InfiniteList xs _) (Positive n) =
  all (== decimalAmultiple1 xs' n)
      [decimalAmultiple2 xs' n,
       decimalAmultiple3 xs' n]
  where xs' = map getPositive xs

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

-- La comparación es
--    λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))
--    21731
--    (0.45 secs, 486,085,256 bytes)
--    λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))
--    21731
--    (0.32 secs, 485,563,664 bytes)
--    λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))
--    21731
--    (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal
-- ===============================

multipleAdecimal1  :: [Integer] -> [Integer] -> Integer
multipleAdecimal1 xs ns = aux xs (reverse ns)
  where aux (y:ys) (m:ms) = m + y * (aux ys ms)
        aux _ _           = 0

-- 2ª solución de multipleAdecimal
-- ===============================

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer
multipleAdecimal2 bs xs =
  sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

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

-- La propiedad es
prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool
prop_multipleAdecimal (InfiniteList xs _) ys =
  multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'
  where xs' = map getPositive xs
        ys' = map getPositive ys

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

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

-- La comparación es
--    λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))
--    35660
--    (0.14 secs, 179,522,152 bytes)
--    λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))
--    35660
--    (0.22 secs, 243,368,664 bytes)

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

-- La primera propiedad es
prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_inversas (InfiniteList xs _) (Positive n) =
  multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n
  where xs' = map getPositive xs

-- Su comprobación es
--    λ> quickCheck prop_inversas
--    +++ OK, passed 100 tests.

-- la 2ª propiedad es
prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_coeficientes (InfiniteList xs _) (Positive n) =
  and [0 <= c && c < b | (c,b) <- zip cs xs']
  where xs' = map getPositive xs
        cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es
--    λ> quickCheck prop_coeficientes
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

import Test.QuickCheck

import Data.List (unfoldr)

-- 1ª solución de decimalAmultiple

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

decimalAmultiple1 :: [Integer] -> Integer -> [Integer]

decimalAmultiple1 bs n = reverse (aux bs n)

where aux _ 0 = []

aux (d:ds) m = r : aux ds q

where (q,r) = quotRem m d

-- 2ª solución de decimalAmultiple

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

decimalAmultiple2 :: [Integer] -> Integer -> [Integer]

decimalAmultiple2 bs n = aux bs n []

where aux _ 0 xs = xs

aux (d:ds) m xs = aux ds q (r:xs)

where (q,r) = quotRem m d

-- 3ª solución de decimalAmultiple

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

decimalAmultiple3 :: [Integer] -> Integer -> [Integer]

decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))

where f (_ ,0) = Nothing

f ((y:ys),m) = Just (r,(ys,q))

where (q,r) = quotRem m y

-- Comprobación de equivalencia

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

-- La propiedad es

prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_decimalAmultiple (InfiniteList xs _) (Positive n) =

all (== decimalAmultiple1 xs' n)

[decimalAmultiple2 xs' n,

decimalAmultiple3 xs' n]

where xs' = map getPositive xs

-- Comparación de eficiencia de decimalAmultiple

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

-- La comparación es

-- λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))

-- 21731

-- (0.45 secs, 486,085,256 bytes)

-- λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))

-- 21731

-- (0.32 secs, 485,563,664 bytes)

-- λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))

-- 21731

-- (0.44 secs, 487,649,768 bytes)

-- 1ª solución de multipleAdecimal

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

multipleAdecimal1 :: [Integer] -> [Integer] -> Integer

multipleAdecimal1 xs ns = aux xs (reverse ns)

where aux (y:ys) (m:ms) = m + y * (aux ys ms)

aux _ _ = 0

-- 2ª solución de multipleAdecimal

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

multipleAdecimal2 :: [Integer] -> [Integer] -> Integer

multipleAdecimal2 bs xs =

sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))

-- Comprobación de equivalencia de multipleAdecimal

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

-- La propiedad es

prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool

prop_multipleAdecimal (InfiniteList xs _) ys =

multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'

where xs' = map getPositive xs

ys' = map getPositive ys

-- La comprobación es

-- λ> quickCheck prop_multipleAdecimal

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de multipleAdecimal

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

-- La comparación es

-- λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))

-- 35660

-- (0.14 secs, 179,522,152 bytes)

-- λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))

-- 35660

-- (0.22 secs, 243,368,664 bytes)

-- Comprobación de las propiedades

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

-- La primera propiedad es

prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_inversas (InfiniteList xs _) (Positive n) =

multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n

where xs' = map getPositive xs

-- Su comprobación es

-- λ> quickCheck prop_inversas

-- +++ OK, passed 100 tests.

-- la 2ª propiedad es

prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool

prop_coeficientes (InfiniteList xs _) (Positive n) =

and [0 <= c && c < b | (c,b) <- zip cs xs']

where xs' = map getPositive xs

cs = reverse (decimalAmultiple1 xs' n)

-- Su comprobación es

-- λ> quickCheck prop_coeficientes

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Soluciones

Referencias

Soluciones

Soluciones

Referencia

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico