Otra conjetura de Goldbach en Haskell

En 1752 Goldbach le escribió un carta a Euler en la que conjeturaba que todo número impar compuesto se puede escribir como la suma de un primo y el doble de un cuadrado. Por ejemplo,

   9 =  7 + 2×1^2
   15 =  7 + 2×2^2
   21 =  3 + 2×3^2
   25 =  7 + 2×3^2
   27 = 19 + 2×2^2
   33 = 31 + 2×1^2

9 = 7 + 2×1^2

15 = 7 + 2×2^2

21 = 3 + 2×3^2

25 = 7 + 2×3^2

27 = 19 + 2×2^2

33 = 31 + 2×1^2

En la siguiente relación de ejercicios (elaborada para la asignatura de Informática de 1º del Grado en Matemáticas se comprueba con Haskell que la conjetura es falsa.

-- ---------------------------------------------------------------------
-- Librerías auxiliares
-- ---------------------------------------------------------------------

import Graphics.Gnuplot.Simple
import Data.Numbers.Primes (isPrime, primes)
import qualified Data.Map as M

-- ---------------------------------------------------------------------
-- Ejercicio 1. Definir la función
--    descomposiciones :: Integer -> [(Integer,Integer)]
-- tal que (descomposiciones x) es la lista de los pares (p,y) tales que 
-- x = p+2*y^2 con p primo e y entero. Por ejemplo,
--    descomposiciones 15  ==  [(13,1),(7,2)]
-- ---------------------------------------------------------------------

descomposiciones :: Integer -> [(Integer,Integer)]
descomposiciones x =
    [(p,y) | y <- [0..n]
           , let p = x-2*y^2
           , isPrime p]
    where n = ceiling (sqrt ((fromIntegral x)/2))

-- ---------------------------------------------------------------------
-- Ejercicio 2. Definir, usando descomposiciones, la lista
--    contraejemplos1 :: [Integer]
-- cuyos elementos sean los contraejemplos de la anterior conjetura de
-- Goldbach. Por ejemplo, 
--    take 2 contraejemplos1  ==  [5777,5993]
-- ---------------------------------------------------------------------

contraejemplos1 :: [Integer]
contraejemplos1 = 
    [x | x <- [3,5..]
       , not (isPrime x)
       , null (descomposiciones x)]

-- ---------------------------------------------------------------------
-- Ejercicio 3. Definir la función
--    imparesCompuestos :: [Integer]
-- tal que imparesCompuestos es la lista de los números impares que son
-- compuestos. Por ejemplo,
--    take 10 imparesCompuestos  ==  [9,15,21,25,27,33,35,39,45,49]
-- ---------------------------------------------------------------------

imparesCompuestos :: [Integer]
imparesCompuestos = aux [3,5..] (tail primes)
    where aux (x:xs) (y:ys) | x == y    = aux xs ys
                            | otherwise = x : aux xs (y:ys)

-- ---------------------------------------------------------------------
-- Ejercicio 4. Definir, usando descomposiciones e imparesCompuestos, la
-- lista 
--    contraejemplos2 :: [Integer]
-- cuyos elementos sean los contraejemplos de la anterior conjetura de
-- Goldbach. Por ejemplo, 
--    take 2 contraejemplos2  ==  [5777,5993]
-- ---------------------------------------------------------------------

contraejemplos2 :: [Integer]
contraejemplos2 = 
    [x | x <- imparesCompuestos
       , null (descomposiciones x)]

-- ---------------------------------------------------------------------
-- Ejercicio 5. Definir la función
--    conjetura :: Integer -> Bool
-- tal que (conjetura n) se verifica si n es cumple la conjetura. Por
-- ejemplo, 
--    conjetura 2015  ==  True
--    conjetura 5777  ==  False
-- ---------------------------------------------------------------------

conjetura :: Integer -> Bool
conjetura n = 
    any isPrime (takeWhile (>0) (map (\i -> n - 2*i*i) [1..]))

-- ---------------------------------------------------------------------
-- Ejercicio 6. Definir, usando conjetura, la lista 
--    contraejemplos3 :: [Integer]
-- cuyos elementos sean los contraejemplos de la anterior conjetura de
-- Goldbach. Por ejemplo, 
--    take 2 contraejemplos3  ==  [5777,5993]
-- ---------------------------------------------------------------------

contraejemplos3 :: [Integer]
contraejemplos3 = 
    filter (not . conjetura) (filter (not . isPrime) [3,5..])

-- ---------------------------------------------------------------------
-- Ejercicio 7. Comparar la eficiencia de las 3 definiciones de
-- contraejemplos para calcular el primer contraejemplo.
-- ---------------------------------------------------------------------

-- La comparación es
--    ghci> head contraejemplos1
--    5777
--    (0.31 secs, 184659216 bytes)
--    ghci> head contraejemplos2
--    5777
--    (0.27 secs, 146877536 bytes)
--    ghci> head contraejemplos3
--    5777
--    (0.23 secs, 159419472 bytes)

-- ---------------------------------------------------------------------
-- Ejercicio 8. Definir la función
--    distribucion :: [Integer] -> [(Int,[Integer])]
-- tal que (distribucion xs) es la lista de pares (n,ys) tales que ys es
-- la lista de elementos de xs con n descomposiciones como suma de un
-- primo y es doble del cuadrado de un entero. Por ejemplo,
--    ghci> distribucion [3,5..100]
--    [(1,[3,9,17,27,33,57,65,95]),
--     (2,[5,7,11,15,23,29,35,41,47,51,53,59,71,77,83,87,93,99]),
--     (3,[13,19,21,25,39,43,45,63,67,81,89]),
--     (4,[31,37,49,69,73,75,85,97]),
--     (5,[55]),
--     (6,[61,79,91])]
--     ghci> distribucion [2,4..100]
--     [(0,[6,8,12,14,16,18,22,24,26,28,30,32,36,38,40,42,44,46,48,50,54,
--          56,58,60,62,64,66,68,70,72,76,78,80,82,84,86,88,90,92,94,96,98]),
--      (1,[2,4,10,20,34,52,74,100])]
------------------------------------------------------------------------

distribucion :: [Integer] -> [(Int,[Integer])]
distribucion xs = M.toList (M.map reverse (aux xs M.empty))
    where aux [] d = d
          aux (y:ys) d = aux ys (M.insertWith (++) k [y] d)
                         where k = length (descomposiciones y)

-- ---------------------------------------------------------------------
-- Ejercicio 8. Definir la función
--    frecuencias :: [Integer] -> [(Int,Int)]
-- tal que (frecuencias xs) es la lista de pares (n,y) tales que y es
-- el número de elementos de xs con n descomposiciones como suma de un
-- primo y es doble del cuadrado de un entero. Por ejemplo,
--    frecuencias [3,5..100]  ==  [(1,8),(2,18),(3,11),(4,8),(5,1),(6,3)]
--    frecuencias [2,4..100]  ==  [(0,42),(1,8)]
-- ---------------------------------------------------------------------

frecuencias :: [Integer] -> [(Int,Int)]
frecuencias xs =
    [(n,length ys) | (n,ys) <- distribucion xs] 

-- ---------------------------------------------------------------------
-- Ejercicio 9. Definir el procedimiento
--    dibujoFrecuencias :: Integer -> IO ()
-- tal que (dibujoFrecuencias n) dibuja las frecuencias de [3,5..n] para
-- n igual a 1000, 3000 y 6000.
-- ---------------------------------------------------------------------

dibujoFrecuencias :: Integer -> IO ()
dibujoFrecuencias n = 
    plotLists []
              [frecuencias' [3,5..k] | k <- [n,3*n,6*n]]
    where frecuencias' :: [Integer] -> [(Double,Double)]
          frecuencias' xs = [(fromIntegral x, fromIntegral y)
                             | (x,y) <- frecuencias xs]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

-- ---------------------------------------------------------------------

-- Librerías auxiliares

-- ---------------------------------------------------------------------

import Graphics.Gnuplot.Simple

import Data.Numbers.Primes (isPrime, primes)

import qualified Data.Map as M

-- ---------------------------------------------------------------------

-- Ejercicio 1. Definir la función

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

-- tal que (descomposiciones x) es la lista de los pares (p,y) tales que

-- x = p+2*y^2 con p primo e y entero. Por ejemplo,

-- descomposiciones 15 == [(13,1),(7,2)]

-- ---------------------------------------------------------------------

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

descomposiciones x =

[(p,y) | y <- [0..n]

, let p = x-2*y^2

, isPrime p]

where n = ceiling (sqrt ((fromIntegral x)/2))

-- ---------------------------------------------------------------------

-- Ejercicio 2. Definir, usando descomposiciones, la lista

-- contraejemplos1 :: [Integer]

-- cuyos elementos sean los contraejemplos de la anterior conjetura de

-- Goldbach. Por ejemplo,

-- take 2 contraejemplos1 == [5777,5993]

-- ---------------------------------------------------------------------

contraejemplos1 :: [Integer]

contraejemplos1 =

[x | x <- [3,5..]

, not (isPrime x)

, null (descomposiciones x)]

-- ---------------------------------------------------------------------

-- Ejercicio 3. Definir la función

-- imparesCompuestos :: [Integer]

-- tal que imparesCompuestos es la lista de los números impares que son

-- compuestos. Por ejemplo,

-- take 10 imparesCompuestos == [9,15,21,25,27,33,35,39,45,49]

-- ---------------------------------------------------------------------

imparesCompuestos :: [Integer]

imparesCompuestos = aux [3,5..] (tail primes)

where aux (x:xs) (y:ys) | x == y = aux xs ys

| otherwise = x : aux xs (y:ys)

-- ---------------------------------------------------------------------

-- Ejercicio 4. Definir, usando descomposiciones e imparesCompuestos, la

-- lista

-- contraejemplos2 :: [Integer]

-- cuyos elementos sean los contraejemplos de la anterior conjetura de

-- Goldbach. Por ejemplo,

-- take 2 contraejemplos2 == [5777,5993]

-- ---------------------------------------------------------------------

contraejemplos2 :: [Integer]

contraejemplos2 =

[x | x <- imparesCompuestos

, null (descomposiciones x)]

-- ---------------------------------------------------------------------

-- Ejercicio 5. Definir la función

-- conjetura :: Integer -> Bool

-- tal que (conjetura n) se verifica si n es cumple la conjetura. Por

-- ejemplo,

-- conjetura 2015 == True

-- conjetura 5777 == False

-- ---------------------------------------------------------------------

conjetura :: Integer -> Bool

conjetura n =

any isPrime (takeWhile (>0) (map (\i -> n - 2*i*i) [1..]))

-- ---------------------------------------------------------------------

-- Ejercicio 6. Definir, usando conjetura, la lista

-- contraejemplos3 :: [Integer]

-- cuyos elementos sean los contraejemplos de la anterior conjetura de

-- Goldbach. Por ejemplo,

-- take 2 contraejemplos3 == [5777,5993]

-- ---------------------------------------------------------------------

contraejemplos3 :: [Integer]

contraejemplos3 =

filter (not . conjetura) (filter (not . isPrime) [3,5..])

-- ---------------------------------------------------------------------

-- Ejercicio 7. Comparar la eficiencia de las 3 definiciones de

-- contraejemplos para calcular el primer contraejemplo.

-- ---------------------------------------------------------------------

-- La comparación es

-- ghci> head contraejemplos1

-- 5777

-- (0.31 secs, 184659216 bytes)

-- ghci> head contraejemplos2

-- 5777

-- (0.27 secs, 146877536 bytes)

-- ghci> head contraejemplos3

-- 5777

-- (0.23 secs, 159419472 bytes)

-- ---------------------------------------------------------------------

-- Ejercicio 8. Definir la función

-- distribucion :: [Integer] -> [(Int,[Integer])]

-- tal que (distribucion xs) es la lista de pares (n,ys) tales que ys es

-- la lista de elementos de xs con n descomposiciones como suma de un

-- primo y es doble del cuadrado de un entero. Por ejemplo,

-- ghci> distribucion [3,5..100]

-- [(1,[3,9,17,27,33,57,65,95]),

-- (2,[5,7,11,15,23,29,35,41,47,51,53,59,71,77,83,87,93,99]),

-- (3,[13,19,21,25,39,43,45,63,67,81,89]),

-- (4,[31,37,49,69,73,75,85,97]),

-- (5,[55]),

-- (6,[61,79,91])]

-- ghci> distribucion [2,4..100]

-- [(0,[6,8,12,14,16,18,22,24,26,28,30,32,36,38,40,42,44,46,48,50,54,

-- 56,58,60,62,64,66,68,70,72,76,78,80,82,84,86,88,90,92,94,96,98]),

-- (1,[2,4,10,20,34,52,74,100])]

------------------------------------------------------------------------

distribucion :: [Integer] -> [(Int,[Integer])]

distribucion xs = M.toList (M.map reverse (aux xs M.empty))

where aux [] d = d

aux (y:ys) d = aux ys (M.insertWith (++) k [y] d)

where k = length (descomposiciones y)

-- ---------------------------------------------------------------------

-- Ejercicio 8. Definir la función

-- frecuencias :: [Integer] -> [(Int,Int)]

-- tal que (frecuencias xs) es la lista de pares (n,y) tales que y es

-- el número de elementos de xs con n descomposiciones como suma de un

-- primo y es doble del cuadrado de un entero. Por ejemplo,

-- frecuencias [3,5..100] == [(1,8),(2,18),(3,11),(4,8),(5,1),(6,3)]

-- frecuencias [2,4..100] == [(0,42),(1,8)]

-- ---------------------------------------------------------------------

frecuencias :: [Integer] -> [(Int,Int)]

frecuencias xs =

[(n,length ys) | (n,ys) <- distribucion xs]

-- ---------------------------------------------------------------------

-- Ejercicio 9. Definir el procedimiento

-- dibujoFrecuencias :: Integer -> IO ()

-- tal que (dibujoFrecuencias n) dibuja las frecuencias de [3,5..n] para

-- n igual a 1000, 3000 y 6000.

-- ---------------------------------------------------------------------

dibujoFrecuencias :: Integer -> IO ()

dibujoFrecuencias n =

plotLists []

[frecuencias' [3,5..k] | k <- [n,3*n,6*n]]

where frecuencias' :: [Integer] -> [(Double,Double)]

frecuencias' xs = [(fromIntegral x, fromIntegral y)

| (x,y) <- frecuencias xs]

El dibujo es

Referencias

Problema 46 del Proyecto Euler.
A lesser-known Goldbach conjecture por L. Hodges.
Sucesión A060003 de OEIS.

Referencias

Se puede imprimir o compartir con