Ejercicio

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

Representaciones de un número como suma de dos cuadrados

Soluciones

Factorizaciones de números de Hilbert

Soluciones

Referencias

Números primos de Hilbert

Soluciones

Sumas de dos primos

Soluciones

Referencia

La sucesión de Thue-Morse

Soluciones

Referencias

Entradas recientes

Buscar

Correo electrónico