div

Máximos locales en los números de descomposiciones de Goldbach

PorJosé A. Alonso 20-enero-202027-enero-2020

La conjetura de Goldbach afirma que todo número entero mayor que 2 se puede expresar como suma de dos primos.

Las descomposiciones de Goldbach son las maneras de expresar un número como suma de dos primos. Por ejemplo, el número 10 tiene dos descomposiciones de Goldbach ya que se puede expresar como la suma de 3 y 7 y la suma de 5 y 5.

Definir las funciones

   descomposicionesGoldbach :: Integer -> [(Integer,Integer)]
   numeroGoldbach :: Integer -> Integer
   tieneMaximoLocalGoldbach :: Integer -> Bool

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

numeroGoldbach :: Integer -> Integer

tieneMaximoLocalGoldbach :: Integer -> Bool

tales que

(descomposicionesGoldbach n) es la lista de las descomposiciones de Goldbach del número n. Por ejemplo,

     descomposicionesGoldbach 5   ==  [(2,3)]
     descomposicionesGoldbach 10  ==  [(3,7),(5,5)]
     descomposicionesGoldbach 22  ==  [(3,19),(5,17),(11,11)]
     descomposicionesGoldbach 34  ==  [(3,31),(5,29),(11,23),(17,17)]
     descomposicionesGoldbach 35  ==  []
     descomposicionesGoldbach (9+10^9)  ==  [(2,1000000007)]

descomposicionesGoldbach 5 == [(2,3)]

descomposicionesGoldbach 10 == [(3,7),(5,5)]

descomposicionesGoldbach 22 == [(3,19),(5,17),(11,11)]

descomposicionesGoldbach 34 == [(3,31),(5,29),(11,23),(17,17)]

descomposicionesGoldbach 35 == []

descomposicionesGoldbach (9+10^9) == [(2,1000000007)]

(numeroGolbach n) es el número de descomposiciones de Goldbach del número n. Por ejemplo,

     numeroGoldbach 5         ==  1
     numeroGoldbach 10        ==  2
     numeroGoldbach 22        ==  3
     numeroGoldbach 34        ==  4
     numeroGoldbach 35        ==  0
     numeroGoldbach (9+10^9)  ==  1
     maximum [numeroGoldbach n | n <- [2..3000]]  ==  128

numeroGoldbach 5 == 1

numeroGoldbach 10 == 2

numeroGoldbach 22 == 3

numeroGoldbach 34 == 4

numeroGoldbach 35 == 0

numeroGoldbach (9+10^9) == 1

maximum [numeroGoldbach n | n <- [2..3000]] == 128

(tieneMaximoLocalGoldbach n) se verifica si en n se alcanza un máximo local en el número de descomposiciones de Goldbach; es decir, los números n tales que el número de descomposiciones de Goldbach de n es mayor o igual que las de n-1 y las de n+1. Por ejemplo,

     λ> filter tieneMaximoLocalGoldbach [1..45]
     [1,2,4,5,6,7,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44]

1 2	λ> filter tieneMaximoLocalGoldbach [1..45] [1,2,4,5,6,7,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44]

En el ejemplo anterior se comprueba que en los múltiplos de 6 (es decir, en 6, 12, 18, 24, 30, 36 y 42), el número de descomposiciones de Goldbach alcanza un máximo local. Comprobar con QuickCheck que esta propiedad se cumple en general; es decir, para todo entero positivo n, el número de descomposiciones de Goldbach en 6n es un máximo local.

Soluciones

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

-- Definiciones de descomposicionesGoldbach
-- ========================================

-- 1ª definición
descomposicionesGoldbach1 :: Integer -> [(Integer,Integer)]
descomposicionesGoldbach1 n =
  [(p,n-p) | p <- takeWhile (<= n `div` 2) primes
           , isPrime (n-p)]

-- 2ª definición
descomposicionesGoldbach2 :: Integer -> [(Integer,Integer)]
descomposicionesGoldbach2 n
  | odd n     = [(2,n-2) | isPrime (n-2)]
  | otherwise = [(p,n-p) | p <- takeWhile (<= n `div` 2) primes
                         , isPrime (n-p)]                               

-- Comparación de eficiencia 
--    λ> descomposicionesGoldbach1 (9+10^8)
--    [(2,100000007)]
--    (10.75 secs, 32,177,389,480 bytes)
--    λ> descomposicionesGoldbach2 (9+10^8)
--    [(2,100000007)]
--    (0.01 secs, 3,228,912 bytes)

-- En lo que sigue, usaremos la 2ª definición
descomposicionesGoldbach :: Integer -> [(Integer,Integer)]
descomposicionesGoldbach = descomposicionesGoldbach2

-- Definición de numeroGolbach
-- ===========================

numeroGoldbach :: Integer -> Integer
numeroGoldbach = genericLength . descomposicionesGoldbach

-- Definición de tieneMaximoLocalGoldbach
-- ======================================

tieneMaximoLocalGoldbach :: Integer -> Bool
tieneMaximoLocalGoldbach n =
  numeroGoldbach (n-1) <= x && x >= numeroGoldbach (n+1)
  where x = numeroGoldbach n

-- La propiedad es
prop_Goldbach :: Integer -> Property
prop_Goldbach n =
  n > 0 ==> tieneMaximoLocalGoldbach (6*n)

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

import Data.List (genericLength)

import Data.Numbers.Primes (primes, isPrime)

import Test.QuickCheck

-- Definiciones de descomposicionesGoldbach

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

-- 1ª definición

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

descomposicionesGoldbach1 n =

[(p,n-p) | p <- takeWhile (<= n `div` 2) primes

, isPrime (n-p)]

-- 2ª definición

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

descomposicionesGoldbach2 n

| odd n = [(2,n-2) | isPrime (n-2)]

| otherwise = [(p,n-p) | p <- takeWhile (<= n `div` 2) primes

, isPrime (n-p)]

-- Comparación de eficiencia

-- λ> descomposicionesGoldbach1 (9+10^8)

-- [(2,100000007)]

-- (10.75 secs, 32,177,389,480 bytes)

-- λ> descomposicionesGoldbach2 (9+10^8)

-- [(2,100000007)]

-- (0.01 secs, 3,228,912 bytes)

-- En lo que sigue, usaremos la 2ª definición

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

descomposicionesGoldbach = descomposicionesGoldbach2

-- Definición de numeroGolbach

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

numeroGoldbach :: Integer -> Integer

numeroGoldbach = genericLength . descomposicionesGoldbach

-- Definición de tieneMaximoLocalGoldbach

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

tieneMaximoLocalGoldbach :: Integer -> Bool

tieneMaximoLocalGoldbach n =

numeroGoldbach (n-1) <= x && x >= numeroGoldbach (n+1)

where x = numeroGoldbach n

-- La propiedad es

prop_Goldbach :: Integer -> Property

prop_Goldbach n =

n > 0 ==> tieneMaximoLocalGoldbach (6*n)

-- La comprobación es

-- λ> quickCheck prop_Goldbach

-- +++ OK, passed 100 tests.

Otras soluciones

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

Referencia

Local maxima of number of Goldbach decompositions en ProofWiki.
Sucesión A061358 de OEIS.

Pensamiento

Te abanicaras
con un madrigal que diga:
en amor el olvido pone la sal.

Antonio Machado

Enteros como sumas de tres coprimos.

PorJosé A. Alonso 8-enero-202015-enero-2020

Dos números enteros son coprimos (o primos entre sí) si no tienen ningún factor primo en común. Por ejemplo, 4 y 15 son coprimos.

Una terna coprima es una terna (a,b,c) tal que

a y b son coprimos,
a y c son coprimos y
b y c son coprimos.

Por ejemplo, (3,4,5) es una terna coprima.

Definir la función

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

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

tal que (sumas3coprimos n) es la lista de las ternas coprimas cuya suma es n. Por ejemplo,

   sumas3coprimos 10  ==  [(2,3,5)]
   sumas3coprimos 11  ==  []
   sumas3coprimos 12  ==  [(2,3,7),(3,4,5)]
   length (sumas3coprimos 4000)  ==  546146

sumas3coprimos 10 == [(2,3,5)]

sumas3coprimos 11 == []

sumas3coprimos 12 == [(2,3,7),(3,4,5)]

length (sumas3coprimos 4000) == 546146

Comprobar con QuickCheck que todo número entero mayor que 17 se puede escribir como suma de alguna terna coprima; es decir, para todo entero n, (sumas3coprimos2 (18 + abs n)) tiene algún elemento.

Soluciones

import Test.QuickCheck

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

sumas3coprimos :: Integer -> [(Integer,Integer,Integer)]
sumas3coprimos n =
  [(a,b,c) | a <- [2..n]
           , b <- [a+1..n]
           , c <- [b+1..n]
           , a + b + c == n
           , gcd a b == 1
           , gcd a c == 1
           , gcd b c == 1]

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

sumas3coprimos2 :: Integer -> [(Integer,Integer,Integer)]
sumas3coprimos2 n =
  [(a,b,c) | a <- [2..n `div` 3]
           , b <- [a+1..(n - a) `div` 2]
           , gcd a b == 1
           , let c = n - a - b  
           , gcd a c == 1
           , gcd b c == 1]

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

-- La propiedad de equivalencia es
prop_sumas3coprimos_equiv :: Integer -> Property
prop_sumas3coprimos_equiv n =
  n > 0 ==> sumas3coprimos n == sumas3coprimos2 n

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

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

--    λ> length (sumas3coprimos 400)
--    5345
--    (4.16 secs, 2,894,799,744 bytes)
--    λ> length (sumas3coprimos2 400)
--    5345
--    (0.06 secs, 16,565,136 bytes)

-- Propiedad
-- =========

-- La propiedad
prop_sumas3coprimos :: Integer -> Bool
prop_sumas3coprimos n =
  not (null (sumas3coprimos2 (18 + abs n)))

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

import Test.QuickCheck

-- 1ª solución

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

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

sumas3coprimos n =

[(a,b,c) | a <- [2..n]

, b <- [a+1..n]

, c <- [b+1..n]

, a + b + c == n

, gcd a b == 1

, gcd a c == 1

, gcd b c == 1]

-- 2ª solución

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

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

sumas3coprimos2 n =

[(a,b,c) | a <- [2..n `div` 3]

, b <- [a+1..(n - a) `div` 2]

, gcd a b == 1

, let c = n - a - b

, gcd a c == 1

, gcd b c == 1]

-- Equivalencia de las definiciones

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

-- La propiedad de equivalencia es

prop_sumas3coprimos_equiv :: Integer -> Property

prop_sumas3coprimos_equiv n =

n > 0 ==> sumas3coprimos n == sumas3coprimos2 n

-- La comprobación es

-- λ> quickCheck prop_sumas3coprimos_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> length (sumas3coprimos 400)

-- 5345

-- (4.16 secs, 2,894,799,744 bytes)

-- λ> length (sumas3coprimos2 400)

-- 5345

-- (0.06 secs, 16,565,136 bytes)

-- Propiedad

-- =========

-- La propiedad

prop_sumas3coprimos :: Integer -> Bool

prop_sumas3coprimos n =

not (null (sumas3coprimos2 (18 + abs n)))

-- La comprobación es

-- λ> quickCheck prop_sumas3coprimos

-- +++ OK, passed 100 tests.

Referencias

Basado en Integers as sum of three pairwise coprime integers de ProofWiki.

Pensamiento

Todo amor es fantasía;
él inventa el año, el día,
la hora y su melodía;
inventa el amante y, más
la amada. No prueba nada,
contra el amor, que la amada
no haya existido jamás.

Antonio Machado

Factorizaciones de 4n+1

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

Sea S el conjunto

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

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

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

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

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

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

Definir las funciones

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

esIrreducible :: Integer -> Bool

factorizaciones :: Integer -> [[Integer]]

conFactorizacionNoUnica :: [Integer]

tales que

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

     esIrreducible 9   ==  True
     esIrreducible 45  ==  False

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

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

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

factorizaciones 9 == [[9]]

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

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

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

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

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

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

Soluciones

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

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

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

esIrreducible :: Integer -> Bool

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

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

-- ejemplo,

-- divisores 9 == [9]

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

divisores :: Integer -> [Integer]

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

factorizaciones :: Integer -> [[Integer]]

factorizaciones 1 = [[1]]

factorizaciones n

| esIrreducible n = [[n]]

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

, esIrreducible d

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

, d <= e]

conFactorizacionNoUnica :: [Integer]

conFactorizacionNoUnica =

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

Pensamiento

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

Antonio Machado

Derivada aritmética

PorJosé A. Alonso 31-diciembre-20197-enero-2020

La derivada aritmética es una función definida sobre los números naturales por analogía con la regla del producto para el cálculo de las derivadas usada en análisis.

Para un número natural n su derivada D(n) se define por

   D(0)  = 0
   D(1)  = 0
   D(p)  = 1, si p es primo
   D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

D(0) = 0

D(1) = 0

D(p) = 1, si p es primo

D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

Por ejemplo,

   D(6)  = D(2*3) = D(2)*3 + 2*D(3) = 1*3 + 2*1 =  5
   D(12) = D(2*6) = D(2)*6 + 2*D(6) = 1*6 + 2*5 = 16

1 2	D(6) = D(23) = D(2)3 + 2D(3) = 13 + 21 = 5 D(12) = D(26) = D(2)6 + 2D(6) = 16 + 25 = 16

Definir la función

   derivada :: Integer -> Integer

1	derivada :: Integer -> Integer

tal que (derivada n) es la derivada aritmética de n. Por ejemplo,

   derivada  6  ==  5
   derivada 12  ==  16
   maximum [derivada n | n <- [1..60000]]  ==  380928

derivada 6 == 5

derivada 12 == 16

maximum [derivada n | n <- [1..60000]] == 380928

Comprobar con QuickCheck que si x es un número entero positivo y su descomposición en factores primos es

   x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

1	x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

entonces la derivada de x es

   x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

1	x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

Nota: No usar en la definición la propiedad que hay que comprobar.

Soluciones

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

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

derivada :: Integer -> Integer
derivada 0 = 0
derivada 1 = 0
derivada n | esPrimo n = 1
           | otherwise = (derivada a) * b + a * (derivada b)
  where a = menorFactor n
        b = n `div` a

-- (esPrimo n) se verifica si n es primo. Por ejemplo,
--    esPrimo 5  ==  True
--    esPrimo 6  ==  False
esPrimo :: Integer -> Bool
esPrimo 0 = False
esPrimo 1 = False
esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por
-- ejemplo, 
--    menorFactor 6   ==  2
--    menorFactor 7   ==  7
--    menorFactor 15  ==  3
menorFactor :: Integer -> Integer
menorFactor n
  | even n = 2
  | otherwise = head [x | x <- [3,5..]
                        , n `mod` x == 0]

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

derivada2 :: Integer -> Integer
derivada2 0 = 0
derivada2 1 = 0
derivada2 n | isPrime n = 1
            | otherwise = (derivada2 a) * b + a * (derivada2 b)
  where (a:_) = primeFactors n
        b     = n `div` a

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

--    λ> maximum [derivada n | n <- [1..10000]]
--    53248
--    (1.59 secs, 1,091,452,552 bytes)
--    λ> maximum [derivada2 n | n <- [1..10000]]
--    53248
--    (0.17 secs, 457,819,120 bytes)

-- Propiedad
-- =========

-- La propiedad es
prop_derivada :: Integer -> Property
prop_derivada x =
  x > 0 ==>
  derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de
-- la descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- Su comprobación es
--    λ> quickCheck prop_derivada
--    +++ OK, passed 100 tests.

import Data.List (genericLength, group)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

-- 1ª solución

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

derivada :: Integer -> Integer

derivada 0 = 0

derivada 1 = 0

derivada n | esPrimo n = 1

| otherwise = (derivada a) * b + a * (derivada b)

where a = menorFactor n

b = n `div` a

-- (esPrimo n) se verifica si n es primo. Por ejemplo,

-- esPrimo 5 == True

-- esPrimo 6 == False

esPrimo :: Integer -> Bool

esPrimo 0 = False

esPrimo 1 = False

esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por

-- ejemplo,

-- menorFactor 6 == 2

-- menorFactor 7 == 7

-- menorFactor 15 == 3

menorFactor :: Integer -> Integer

menorFactor n

| even n = 2

| otherwise = head [x | x <- [3,5..]

, n `mod` x == 0]

-- 2ª solución

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

derivada2 :: Integer -> Integer

derivada2 0 = 0

derivada2 1 = 0

derivada2 n | isPrime n = 1

| otherwise = (derivada2 a) * b + a * (derivada2 b)

where (a:_) = primeFactors n

b = n `div` a

-- Comparación de eficiencia

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

-- λ> maximum [derivada n | n <- [1..10000]]

-- 53248

-- (1.59 secs, 1,091,452,552 bytes)

-- λ> maximum [derivada2 n | n <- [1..10000]]

-- 53248

-- (0.17 secs, 457,819,120 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_derivada :: Integer -> Property

prop_derivada x =

x > 0 ==>

derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de

-- la descomposición prima de x. Por ejemplo,

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

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

factorizacion n =

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

-- Su comprobación es

-- λ> quickCheck prop_derivada

-- +++ OK, passed 100 tests.

Referencias

Arithmetic derivative en Wikipedia.
Sucesión A003415 de OEIS.

Pensamiento

En ese jardín, Guiomar,
el mutuo jardín que inventan
dos corazones al par,
se funden y complementan
nuestras horas.

Antonio Machado

El teorema de Navidad de Fermat

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

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

Definir las funciones

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

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

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

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

representaciones 20 == [(2,4)]

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

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

representaciones 100000147984 == [(0,316228)]

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

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

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

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

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

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

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

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

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

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

Soluciones

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

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

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

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

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

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

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer 
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
             
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª definición de representaciones

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

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

representaciones n =

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

-- 2ª definición de representaciones

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

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

representaciones2 n =

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

, let z = n - x*x

, esCuadrado z]

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

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

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

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz 0 = 0

raiz 1 = 1

raiz x = aux (0,x)

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

| c == a = a

| d < x = aux (c,b)

| otherwise = aux (a,c)

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

d = c^2

-- 3ª definición de representaciones

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

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

representaciones3 n =

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

, let z = n - x*x

, esCuadrado3 z]

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

-- ejemplo,

-- esCuadrado3 25 == True

-- esCuadrado3 26 == False

esCuadrado3 :: Integer -> Bool

esCuadrado3 x = x == y * y

where y = raiz3 x

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

-- raiz3 25 == 5

-- raiz3 24 == 4

-- raiz3 26 == 5

raiz3 :: Integer -> Integer

raiz3 x = floor (sqrt (fromIntegral x))

-- 4ª definición de representaciones

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

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

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

where aux x y

| x > y = []

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

LT -> aux (x + 1) y

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

GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones

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

-- La propiedad es

prop_representaciones_equiv :: (Positive Integer) -> Bool

prop_representaciones_equiv (Positive n) =

representaciones n == representaciones2 n &&

representaciones2 n == representaciones3 n &&

representaciones3 n == representaciones4 n

-- La comprobación es

-- λ> quickCheck prop_representaciones_equiv

-- +++ OK, passed 100 tests.

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

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

-- λ> representaciones 3025

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

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

-- λ> representaciones2 3025

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

-- (0.00 secs, 867,944 bytes)

-- λ> representaciones3 3025

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

-- (0.00 secs, 173,512 bytes)

-- λ> representaciones4 3025

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

-- (0.00 secs, 423,424 bytes)

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

-- 6

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

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

-- 6

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

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

-- 6

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

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

-- 7

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

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

-- 7

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

-- Definición de primosImparesConRepresentacionUnica

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

primosImparesConRepresentacionUnica :: [Integer]

primosImparesConRepresentacionUnica =

[x | x <- tail primes

, length (representaciones4 x) == 1]

-- Definición de primos4nM1

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

primos4nM1 :: [Integer]

primos4nM1 = [x | x <- primes

, x `mod` 4 == 1]

-- Teorema de Navidad de Fermat

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

-- La propiedad es

prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool

prop_teoremaDeNavidadDeFermat (Positive n) =

primosImparesConRepresentacionUnica !! n == primos4nM1 !! n

-- La comprobación es

-- λ> quickCheck prop_teoremaDeNavidadDeFermat

-- +++ OK, passed 100 tests.

Pensamiento

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

Antonio Machado

Máximos locales en los números de descomposiciones de Goldbach

Soluciones

Otras soluciones

Referencia

Pensamiento

Enteros como sumas de tres coprimos.

Soluciones

Referencias

Pensamiento

Factorizaciones de 4n+1

Soluciones

Pensamiento

El teorema de Navidad de Fermat

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico