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,