take

Cálculo de pi mediante la fracción continua de Lange

PorJosé A. Alonso 10-abril-202026-marzo-2022

En 1999, L.J. Lange publicó el artículo An elegant new continued fraction for π.

En el primer teorema del artículo se demuestra la siguiente expresión de π mediante una fracción continua
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange$

La primeras aproximaciones son

   a(1) = 3+1                = 4.0
   a(2) = 3+(1/(6+9))        = 3.066666666666667
   a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

a(1) = 3+1 = 4.0

a(2) = 3+(1/(6+9)) = 3.066666666666667

a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

Definir las funciones

   aproximacionPi :: Int -> Double
   grafica        :: [Int] -> IO ()

1 2	aproximacionPi :: Int -> Double grafica :: [Int] -> IO ()

tales que

(aproximacionPi n) es la n-ésima aproximación de pi con la fracción continua de Lange. Por ejemplo,

     aproximacionPi 1     ==  4.0
     aproximacionPi 2     ==  3.066666666666667
     aproximacionPi 3     ==  3.158974358974359
     aproximacionPi 10    ==  3.141287132741557
     aproximacionPi 100   ==  3.141592398533554
     aproximacionPi 1000  ==  3.1415926533392926

aproximacionPi 1 == 4.0

aproximacionPi 2 == 3.066666666666667

aproximacionPi 3 == 3.158974358974359

aproximacionPi 10 == 3.141287132741557

aproximacionPi 100 == 3.141592398533554

aproximacionPi 1000 == 3.1415926533392926

(grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi donde k toma los valores de la lista xs. Por ejemplo, (grafica [1..10]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_2$
(grafica [10..100]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_3$
y (grafica [100..200]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_4$

Soluciones

import Graphics.Gnuplot.Simple
 
-- fraccionPi es la representación de la fracción continua de pi como un
-- par de listas infinitas.
fraccionPi :: [(Integer, Integer)]
fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción
-- continua fc (como un par de listas).  
aproximacionFC :: Int -> [(Integer, Integer)] -> Double
aproximacionFC n =
  foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double
aproximacionPi n =
  aproximacionFC n fraccionPi
 
grafica :: [Int] -> IO ()
grafica xs = 
  plotList [Key Nothing]
           [(k,aproximacionPi k) | k <- xs]

import Graphics.Gnuplot.Simple

-- fraccionPi es la representación de la fracción continua de pi como un

-- par de listas infinitas.

fraccionPi :: [(Integer, Integer)]

fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción

-- continua fc (como un par de listas).

aproximacionFC :: Int -> [(Integer, Integer)] -> Double

aproximacionFC n =

foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double

aproximacionPi n =

aproximacionFC n fraccionPi

grafica :: [Int] -> IO ()

grafica xs =

plotList [Key Nothing]

[(k,aproximacionPi k) | k <- xs]

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>

Cálculo de dígitos de pi y su distribución

PorJosé A. Alonso 8-abril-202021-abril-2020

Se pueden generar los dígitos de Pi, como se explica en el artículo Unbounded spigot algorithms for the digits of pi c0on la función digitosPi definida por

   digitosPi :: [Integer]
   digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) = 
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

Por ejemplo,

   λ> take 25 digitosPi
   [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

1 2	λ> take 25 digitosPi [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

La distribución de los primeros 25 dígitos de pi es [0,2,3,5,3,3,3,1,2,3] ya que el 0 no aparece, el 1 ocurre 2 veces, el 3 ocurre 3 veces, el 4 ocurre 5 veces, …

Usando digitosPi, definir las siguientes funciones

   distribucionDigitosPi :: Int -> [Int]
   frecuenciaDigitosPi   :: Int -> [Double]

1 2	distribucionDigitosPi :: Int -> [Int] frecuenciaDigitosPi :: Int -> [Double]

tales que

(distribucionDigitosPi n) es la distribución de los n primeros dígitos de pi. Por ejemplo,

     λ> distribucionDigitosPi 10
     [0,2,1,2,1,2,1,0,0,1]
     λ> distribucionDigitosPi 100
     [8,8,12,12,10,8,9,8,12,13]
     λ> distribucionDigitosPi 1000
     [93,116,103,103,93,97,94,95,101,105]
     λ> distribucionDigitosPi 5000
     [466,531,496,460,508,525,513,488,492,521]

λ> distribucionDigitosPi 10

[0,2,1,2,1,2,1,0,0,1]

λ> distribucionDigitosPi 100

[8,8,12,12,10,8,9,8,12,13]

λ> distribucionDigitosPi 1000

[93,116,103,103,93,97,94,95,101,105]

λ> distribucionDigitosPi 5000

[466,531,496,460,508,525,513,488,492,521]

(frecuenciaDigitosPi n) es la frecuencia de los n primeros dígitos de pi. Por ejemplo,

   λ> frecuenciaDigitosPi 10
   [0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]
   λ> frecuenciaDigitosPi 100
   [8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]
   λ> frecuenciaDigitosPi 1000
   [9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]
   λ> frecuenciaDigitosPi 5000
   [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

λ> frecuenciaDigitosPi 10

[0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]

λ> frecuenciaDigitosPi 100

[8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]

λ> frecuenciaDigitosPi 1000

[9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]

λ> frecuenciaDigitosPi 5000

[9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

Soluciones

import Data.Array
import Data.List (group, sort)

digitosPi :: [Integer]
digitosPi = g(1,0,1,1,3,3) where
  g (q,r,t,k,n,l) = 
    if 4*q+r-t < n*t
    then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
    else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

-- 1ª definición
-- =============

distribucionDigitosPi :: Int -> [Int]
distribucionDigitosPi n =
  elems (accumArray (+) 0 (0,9) [ (i,1)
                                | i <- take n digitosPi]) 

frecuenciaDigitosPi :: Int -> [Double]
frecuenciaDigitosPi n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]
  where m = fromIntegral n

-- 2ª definición
-- =============

distribucionDigitosPi2 :: Int -> [Int]
distribucionDigitosPi2 n =
  [length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]
frecuenciaDigitosPi2 n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]
  where m = fromIntegral n

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

--    λ> last (take 5000 digitosPi)
--    2
--    (4.47 secs, 3,927,848,448 bytes)
--    λ> frecuenciaDigitosPi 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.01 secs, 0 bytes)
--    λ> frecuenciaDigitosPi2 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.02 secs, 0 bytes)

import Data.Array

import Data.List (group, sort)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

-- 1ª definición

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

distribucionDigitosPi :: Int -> [Int]

distribucionDigitosPi n =

elems (accumArray (+) 0 (0,9) [ (i,1)

| i <- take n digitosPi])

frecuenciaDigitosPi :: Int -> [Double]

frecuenciaDigitosPi n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]

where m = fromIntegral n

-- 2ª definición

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

distribucionDigitosPi2 :: Int -> [Int]

distribucionDigitosPi2 n =

[length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]

frecuenciaDigitosPi2 n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]

where m = fromIntegral n

-- Comparación de eficiencia

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

-- λ> last (take 5000 digitosPi)

-- 2

-- (4.47 secs, 3,927,848,448 bytes)

-- λ> frecuenciaDigitosPi 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.01 secs, 0 bytes)

-- λ> frecuenciaDigitosPi2 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.02 secs, 0 bytes)

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>

Cálculo de pi mediante los métodos de Gregory-Leibniz y de Beeler

PorJosé A. Alonso 7-abril-202026-marzo-2022

La fórmula de Gregory-Leibniz para calcular pi es

y la de Beeler es

Definir las funciones

   aproximaPiGL     :: Int -> Double
   aproximaPiBeeler :: Int -> Double
   graficas         :: [Int] -> IO ()

aproximaPiGL :: Int -> Double

aproximaPiBeeler :: Int -> Double

graficas :: [Int] -> IO ()

tales que

(aproximaPiGL n) es la aproximación de pi con los primeros n términos de la fórmula de Gregory-Leibniz. Por ejemplo,

     aproximaPiGL 1       ==  4.0
     aproximaPiGL 2       ==  2.666666666666667
     aproximaPiGL 3       ==  3.466666666666667
     aproximaPiGL 10      ==  3.0418396189294032
     aproximaPiGL 100     ==  3.1315929035585537
     aproximaPiGL 1000    ==  3.140592653839794
     aproximaPiGL 10000   ==  3.1414926535900345
     aproximaPiGL 100000  ==  3.1415826535897198

aproximaPiGL 1 == 4.0

aproximaPiGL 2 == 2.666666666666667

aproximaPiGL 3 == 3.466666666666667

aproximaPiGL 10 == 3.0418396189294032

aproximaPiGL 100 == 3.1315929035585537

aproximaPiGL 1000 == 3.140592653839794

aproximaPiGL 10000 == 3.1414926535900345

aproximaPiGL 100000 == 3.1415826535897198

(aproximaPiBeeler n) es la aproximación de pi con los primeros n términos de la fórmula de Beeler. Por ejemplo,

     aproximaPiBeeler 1   ==  2.0
     aproximaPiBeeler 2   ==  2.6666666666666665
     aproximaPiBeeler 3   ==  2.933333333333333
     aproximaPiBeeler 10  ==  3.140578169680337
     aproximaPiBeeler 60  ==  3.141592653589793
     pi                   ==  3.141592653589793

aproximaPiBeeler 1 == 2.0

aproximaPiBeeler 2 == 2.6666666666666665

aproximaPiBeeler 3 == 2.933333333333333

aproximaPiBeeler 10 == 3.140578169680337

aproximaPiBeeler 60 == 3.141592653589793

pi == 3.141592653589793

(graficas xs) dibuja la gráfica de las k-ésimas aproximaciones de pi, donde k toma los valores de la lista xs, con las fórmulas de Gregory-Leibniz y de Beeler. Por ejemplo, (graficas [1..25]) dibuja

donde la línea morada corresponde a la aproximación de Gregory-Leibniz y la verde a la de Beeler.

Soluciones

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL
-- ============================

-- 1ª definición de aproximaPiGL
aproximaPiGL :: Int -> Double
aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]
  where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL
aproximaPiGL2 :: Int -> Double
aproximaPiGL2 n =
  4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL
aproximaPiGL3 :: Int -> Double
aproximaPiGL3 n =
  4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL
aproximaPiGL4 :: Int -> Double
aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]
serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)
  where numeradores   = cycle [4,-4]
        denominadores = [1,3..]

-- Definición de aproximaPiBeeler
aproximaPiBeeler :: Int -> Double
aproximaPiBeeler n = 2 * aux (fromIntegral n) 1
  where
    aux :: Double -> Double -> Double 
    aux n k | n == k    = 1
            | otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas
graficas :: [Int] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
             [[(k,aproximaPiGL k)     | k <- xs],
              [(k,aproximaPiBeeler k) | k <- xs]]

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL

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

-- 1ª definición de aproximaPiGL

aproximaPiGL :: Int -> Double

aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]

where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL

aproximaPiGL2 :: Int -> Double

aproximaPiGL2 n =

4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL

aproximaPiGL3 :: Int -> Double

aproximaPiGL3 n =

4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL

aproximaPiGL4 :: Int -> Double

aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]

serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)

where numeradores = cycle [4,-4]

denominadores = [1,3..]

-- Definición de aproximaPiBeeler

aproximaPiBeeler :: Int -> Double

aproximaPiBeeler n = 2 * aux (fromIntegral n) 1

where

aux :: Double -> Double -> Double

aux n k | n == k = 1

| otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas

graficas :: [Int] -> IO ()

graficas xs =

plotLists [Key Nothing]

[[(k,aproximaPiGL k) | k <- xs],

[(k,aproximaPiBeeler k) | k <- xs]]

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>

Conjetura de Goldbach

PorJosé A. Alonso 1-abril-202026-marzo-2022

Una forma de la conjetura de Golbach afirma que todo entero mayor que 1 se puede escribir como la suma de uno, dos o tres números primos.

Si se define el índice de Goldbach de n > 1 como la mínima cantidad de primos necesarios para que su suma sea n, entonces la conjetura de Goldbach afirma que todos los índices de Goldbach de los enteros mayores que 1 son menores que 4.

Definir las siguientes funciones

   indiceGoldbach  :: Int -> Int
   graficaGoldbach :: Int -> IO ()

1 2	indiceGoldbach :: Int -> Int graficaGoldbach :: Int -> IO ()

tales que

(indiceGoldbach n) es el índice de Goldbach de n. Por ejemplo,

     indiceGoldbach 2                        ==  1
     indiceGoldbach 4                        ==  2
     indiceGoldbach 27                       ==  3
     sum (map indiceGoldbach [2..5000])      ==  10619
     maximum (map indiceGoldbach [2..5000])  ==  3

indiceGoldbach 2 == 1

indiceGoldbach 4 == 2

indiceGoldbach 27 == 3

sum (map indiceGoldbach [2..5000]) == 10619

maximum (map indiceGoldbach [2..5000]) == 3

(graficaGoldbach n) dibuja la gráfica de los índices de Goldbach de los números entre 2 y n. Por ejemplo, (graficaGoldbach 150) dibuja

Comprobar con QuickCheck la conjetura de Goldbach anterior.

Soluciones

import Data.Array
import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
import Test.QuickCheck


-- 1ª definición
-- =============

indiceGoldbach :: Int -> Int
indiceGoldbach n =
  minimum (map length (particiones n))

particiones :: Int -> [[Int]]
particiones n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
    where f 0 = [[]]
          f m = [x:y | x <- xs, 
                       y <- v ! (m-x), 
                       [x] >= take 1 y]
            where xs = reverse (takeWhile (<= m) primes)

-- 2ª definición
-- =============

indiceGoldbach2 :: Int -> Int
indiceGoldbach2 x =
  head [n | n <- [1..], esSumaDe x n]

-- (esSumaDe x n) se verifica si x se puede escribir como la suma de n
-- primos. Por ejemplo,
--    esSumaDe 2  1  ==  True
--    esSumaDe 4  1  ==  False
--    esSumaDe 4  2  ==  True
--    esSumaDe 27 2  ==  False
--    esSumaDe 27 3  ==  True
esSumaDe :: Int -> Int -> Bool
esSumaDe x 1 = isPrime x
esSumaDe x n = or [esSumaDe (x-p) (n-1) | p <- takeWhile (<= x) primes]

-- 3ª definición
-- =============

indiceGoldbach3 :: Int -> Int
indiceGoldbach3 x =
  head [n | n <- [1..], esSumaDe3 x n]

esSumaDe3 :: Int -> Int -> Bool
esSumaDe3 x n = a ! (x,n) where
  a = array ((2,1),(x,9)) [((i,j),f i j) | i <- [2..x], j <- [1..9]]
  f i 1 = isPrime i
  f i j = or [a!(i-k,j-1) | k <- takeWhile (<= i) primes]

-- 4ª definición
-- =============

indiceGoldbach4 :: Int -> Int
indiceGoldbach4 n = v ! n where
  v = array (2,n) [(i,f i) | i <- [2..n]]
  f i | isPrime i = 1
      | otherwise = 1 + minimum [v!(i-p) | p <- takeWhile (< (i-1)) primes]

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

--    λ> sum (map indiceGoldbach [2..80])
--    142
--    (2.66 secs, 1,194,330,496 bytes)
--    λ> sum (map indiceGoldbach2 [2..80])
--    142
--    (0.01 secs, 1,689,944 bytes)
--    λ> sum (map indiceGoldbach3 [2..80])
--    142
--    (0.03 secs, 27,319,296 bytes)
--    λ> sum (map indiceGoldbach4 [2..80])
--    142
--    (0.03 secs, 47,823,656 bytes)
--    
--    λ> sum (map indiceGoldbach2 [2..1000])
--    2030
--    (0.10 secs, 200,140,264 bytes)
--    λ> sum (map indiceGoldbach3 [2..1000])
--    2030
--    (3.10 secs, 4,687,467,664 bytes)

-- Gráfica
-- =======

graficaGoldbach :: Int -> IO ()
graficaGoldbach n =
  plotList [ Key Nothing
           , XRange (2,fromIntegral n)
           , PNG ("Conjetura_de_Goldbach_" ++ show n ++ ".png")
           ]
           [indiceGoldbach2 k | k <- [2..n]]

-- Comprobación de la conjetura de Goldbach
-- ========================================

-- La propiedad es
prop_Goldbach :: Int -> Property
prop_Goldbach x =
  x >= 2 ==> indiceGoldbach2 x < 4

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

100

101

102

103

104

105

106

107

import Data.Array

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

import Test.QuickCheck

-- 1ª definición

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

indiceGoldbach :: Int -> Int

indiceGoldbach n =

minimum (map length (particiones n))

particiones :: Int -> [[Int]]

particiones n = v ! n where

v = array (0,n) [(i,f i) | i <- [0..n]]

where f 0 = [[]]

f m = [x:y | x <- xs,

y <- v ! (m-x),

[x] >= take 1 y]

where xs = reverse (takeWhile (<= m) primes)

-- 2ª definición

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

indiceGoldbach2 :: Int -> Int

indiceGoldbach2 x =

head [n | n <- [1..], esSumaDe x n]

-- (esSumaDe x n) se verifica si x se puede escribir como la suma de n

-- primos. Por ejemplo,

-- esSumaDe 2 1 == True

-- esSumaDe 4 1 == False

-- esSumaDe 4 2 == True

-- esSumaDe 27 2 == False

-- esSumaDe 27 3 == True

esSumaDe :: Int -> Int -> Bool

esSumaDe x 1 = isPrime x

esSumaDe x n = or [esSumaDe (x-p) (n-1) | p <- takeWhile (<= x) primes]

-- 3ª definición

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

indiceGoldbach3 :: Int -> Int

indiceGoldbach3 x =

head [n | n <- [1..], esSumaDe3 x n]

esSumaDe3 :: Int -> Int -> Bool

esSumaDe3 x n = a ! (x,n) where

a = array ((2,1),(x,9)) [((i,j),f i j) | i <- [2..x], j <- [1..9]]

f i 1 = isPrime i

f i j = or [a!(i-k,j-1) | k <- takeWhile (<= i) primes]

-- 4ª definición

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

indiceGoldbach4 :: Int -> Int

indiceGoldbach4 n = v ! n where

v = array (2,n) [(i,f i) | i <- [2..n]]

f i | isPrime i = 1

| otherwise = 1 + minimum [v!(i-p) | p <- takeWhile (< (i-1)) primes]

-- Comparación de eficiencia

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

-- λ> sum (map indiceGoldbach [2..80])

-- 142

-- (2.66 secs, 1,194,330,496 bytes)

-- λ> sum (map indiceGoldbach2 [2..80])

-- 142

-- (0.01 secs, 1,689,944 bytes)

-- λ> sum (map indiceGoldbach3 [2..80])

-- 142

-- (0.03 secs, 27,319,296 bytes)

-- λ> sum (map indiceGoldbach4 [2..80])

-- 142

-- (0.03 secs, 47,823,656 bytes)

-- λ> sum (map indiceGoldbach2 [2..1000])

-- 2030

-- (0.10 secs, 200,140,264 bytes)

-- λ> sum (map indiceGoldbach3 [2..1000])

-- 2030

-- (3.10 secs, 4,687,467,664 bytes)

-- Gráfica

-- =======

graficaGoldbach :: Int -> IO ()

graficaGoldbach n =

plotList [ Key Nothing

, XRange (2,fromIntegral n)

, PNG ("Conjetura_de_Goldbach_" ++ show n ++ ".png")

]

[indiceGoldbach2 k | k <- [2..n]]

-- Comprobación de la conjetura de Goldbach

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

-- La propiedad es

prop_Goldbach :: Int -> Property

prop_Goldbach x =

x >= 2 ==> indiceGoldbach2 x < 4

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

Pensamiento

«La diferencia entre los matemáticos y los físicos es que después de que los físicos prueban un gran resultado piensan que es fantástico, pero después de que los matemáticos prueban un gran resultado piensan que es trivial.»

Lucien Szpiro.

La conjetura de Gilbreath

PorJosé A. Alonso 30-marzo-20206-abril-2020

Partiendo de los 5 primeros números primos y calculando el valor absoluto de la diferencia de cada dos números consecutivos hasta quedarse con un único número se obtiene la siguiente tabla:

   2, 3, 5, 7, 11
   1, 2, 2, 4 
   1, 0, 2
   1, 2 
   1

2, 3, 5, 7, 11

1, 2, 2, 4

1, 0, 2

1, 2

Se observa que todas las filas, salvo la inicial, comienzan con el número 1.

Repitiendo el proceso pero empezando con los 8 primeros números primos se obtiene la siguiente tabla:

   2, 3, 5, 7, 11, 13, 17, 19 
   1, 2, 2, 4,  2,  4,  2  
   1, 0, 2, 2,  2,  2 
   1, 2, 0, 0,  0 
   1, 2, 0, 0 
   1, 2, 0 
   1, 2 
   1

2, 3, 5, 7, 11, 13, 17, 19

1, 2, 2, 4, 2, 4, 2

1, 0, 2, 2, 2, 2

1, 2, 0, 0, 0

1, 2, 0, 0

1, 2, 0

1, 2

Se observa que, de nuevo, todas las filas, salvo la inicial, comienza con el número 1.

La conjetura de Gilbreath afirma que si escribimos la sucesión de números primos completa y después construimos las correspondientes sucesiones formadas por el valor absoluto de la resta de cada pareja de números consecutivos, entonces todas esas filas que obtenemos comienzan siempre por 1.

El objetivo de este ejercicio es comprobar experimentalmente dicha conjetura.

Para la representación, usaremos la simétrica de la que hemos comentado anteriormente; es decir,

    2
    3, 1
    5, 2, 1
    7, 2, 0, 1
   11, 4, 2, 2, 1
   13, 2, 2, 0, 2, 1
   17, 4, 2, 0, 0, 2, 1
   19, 2, 2, 0, 0, 0, 2, 1

3, 1

5, 2, 1

7, 2, 0, 1

11, 4, 2, 2, 1

13, 2, 2, 0, 2, 1

17, 4, 2, 0, 0, 2, 1

19, 2, 2, 0, 0, 0, 2, 1

en la que la primera columna son los números primos y el elemento de la fila i y columna j (con i, j > 1) es el valor absoluto de la diferencia de los elementos (i,j-1) e (i-1,j-1).

Definir las siguientes funciones

   siguiente           :: Integer -> [Integer] -> [Integer]
   triangulo           :: [[Integer]]
   conjeturaGilbreath  :: Int -> Bool

siguiente :: Integer -> [Integer] -> [Integer]

triangulo :: [[Integer]]

conjeturaGilbreath :: Int -> Bool

tales que

(siguiente x ys) es la línea siguiente de la ys que empieza por x en la tabla de Gilbreath; es decir, si ys es [y1,y2,…,yn], entonces (siguiente x ys) es [x,|y1-x|,|y2-|y1-x||,…]. Por ejemplo,

     siguiente  7 [5,2,1]               ==  [7,2,0,1]
     siguiente 29 [23,4,2,0,0,0,0,2,1]  ==  [29,6,2,0,0,0,0,0,2,1]

1 2	siguiente 7 [5,2,1] == [7,2,0,1] siguiente 29 [23,4,2,0,0,0,0,2,1] == [29,6,2,0,0,0,0,0,2,1]

triangulo es el triángulo de Gilbreath. Por ejemplo,

     λ> take 10 triangulo
     [[ 2],
      [ 3,1],
      [ 5,2,1],
      [ 7,2,0,1],
      [11,4,2,2,1],
      [13,2,2,0,2,1],
      [17,4,2,0,0,2,1],
      [19,2,2,0,0,0,2,1],
      [23,4,2,0,0,0,0,2,1],
      [29,6,2,0,0,0,0,0,2,1]]

λ> take 10 triangulo

[[ 2],

[ 3,1],

[ 5,2,1],

[ 7,2,0,1],

[11,4,2,2,1],

[13,2,2,0,2,1],

[17,4,2,0,0,2,1],

[19,2,2,0,0,0,2,1],

[23,4,2,0,0,0,0,2,1],

[29,6,2,0,0,0,0,0,2,1]]

(conjeturaGilbreath n) se verifica si se cumple la conjetura de Gilbreath para los n primeros números primos; es decir, en el triángulo de Gilbreath cuya primera columna son los n primeros números primos, todas las filas a partir de la segunda terminan en 1. Por ejemplo,

     λ> conjeturaGilbreath 1000
     True

1 2	λ> conjeturaGilbreath 1000 True

Soluciones

import Data.Numbers.Primes

siguiente :: Integer -> [Integer] -> [Integer]
siguiente x ys = scanl (\m n -> abs (m-n)) x ys 

triangulo :: [[Integer]]
triangulo = 
  [2] : [siguiente x ys | (x,ys) <- zip (tail primes) triangulo]

conjeturaGilbreath :: Int -> Bool
conjeturaGilbreath n = all p (tail (take n triangulo))
  where p xs = last xs == 1

import Data.Numbers.Primes

siguiente :: Integer -> [Integer] -> [Integer]

siguiente x ys = scanl (\m n -> abs (m-n)) x ys

triangulo :: [[Integer]]

triangulo =

[2] : [siguiente x ys | (x,ys) <- zip (tail primes) triangulo]

conjeturaGilbreath :: Int -> Bool

conjeturaGilbreath n = all p (tail (take n triangulo))

where p xs = last xs == 1

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>

Pensamiento

«La simplicidad es la última sofisticación.»

Leonardo da Vinci.

Cálculo de pi mediante la fracción continua de Lange

Soluciones

Otras soluciones

Cálculo de pi mediante los métodos de Gregory-Leibniz y de Beeler

Soluciones

Otras soluciones

Conjetura de Goldbach

Soluciones

Otras soluciones

Pensamiento

La conjetura de Gilbreath

Soluciones

Otras soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico