Teoría de números – Exercitium

Números de Pentanacci

José A. Alonso, 5-agosto-2022, Ejercicio

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

Los primeros números de Fibonacci son

   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

Los primeros números de Pentanacci son

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

Definir la sucesión

   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

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)

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.

Leer más →

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

José A. Alonso, 4-agosto-2022, Ejercicio

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

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

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)

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.

Leer más →

Representaciones de un número como suma de dos cuadrados

José A. Alonso, 3-agosto-2022, Ejercicio

Definir la función

   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

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.

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.

Leer más →

Factorizaciones de números de Hilbert

José A. Alonso, 2-agosto-2022, Ejercicio

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

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

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.

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.

Leer más →

Números primos de Hilbert

José A. Alonso, 1-agosto-2022, Ejercicio

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]

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

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.

Leer más →

Etiqueta: Teoría de números

Números de Pentanacci

Soluciones

Referencias

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

Soluciones

Referencias

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

Subscribe to Blog via Email

Entradas recientes

Etiquetas