plotLists

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>

El sesgo de Chebyshev

PorJosé A. Alonso 23-marzo-202026-marzo-2022

Un número primo distinto de 2 tiene la forma 4k + 1 o 4k + 3. Chebyshev notó en 1853 que la mayoría de las veces hay más números primos de la forma 4k + 3 que números primos de la forma 4k + 1 menores que un número dado. Esto se llama el sesgo de Chebyshev.

Definir las funciones

   distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]
   empatesRestosModulo4 :: [Integer]
   mayoria1RestosModulo4 :: [Integer]
   grafica_Chebyshev :: Int -> IO ()

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]

empatesRestosModulo4 :: [Integer]

mayoria1RestosModulo4 :: [Integer]

grafica_Chebyshev :: Int -> IO ()

tales que

distribucionPrimosModulo4 es la lista de las ternas (p,a,b) tales que p es un números primo, a es la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 y b es la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 7 distribucionPrimosModulo4
     [(2,0,0),(3,0,1),(5,1,1),(7,1,2),(11,1,3),(13,2,3),(17,3,3)]
     λ> distribucionPrimosModulo4 !! (5*10^5)
     (7368791,249888,250112)

λ> take 7 distribucionPrimosModulo4

[(2,0,0),(3,0,1),(5,1,1),(7,1,2),(11,1,3),(13,2,3),(17,3,3)]

λ> distribucionPrimosModulo4 !! (5*10^5)

(7368791,249888,250112)

empatesRestosModulo4 es la lista de los primos p tales que la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 es igual a la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 10 empatesRestosModulo4
     [2,5,17,41,461,26833,26849,26863,26881,26893]
     λ> length (takeWhile (<= 10^6) empatesRestosModulo4)
     112

λ> take 10 empatesRestosModulo4

[2,5,17,41,461,26833,26849,26863,26881,26893]

λ> length (takeWhile (<= 10^6) empatesRestosModulo4)

112

mayoria1RestosModulo4 es la lista de los primos p tales que la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 es mayor que la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo,

     λ> take 10 mayoria1RestosModulo4
     [26861,616841,616849,616877,616897,616909,616933,616943,616951,616961]
     λ> length (takeWhile (<= 10^6) mayoria1RestosModulo4)
     239

λ> take 10 mayoria1RestosModulo4

[26861,616841,616849,616877,616897,616909,616933,616943,616951,616961]

λ> length (takeWhile (<= 10^6) mayoria1RestosModulo4)

239

(graficaChebyshev n) dibuja la gráfica de los puntos (p,b-a) donde p es uno de los n primeros primos impares, a es la cantidad de primos menores o iguales que p congruentes con 1 módulo 4 y b es la cantidad de primos menores o iguales que p congruentes con 3 módulo 4. Por ejemplo, (graficaChebyshev 5000) dibuja la figura

Soluciones

[schedule expon=’2020-03-30′ expat=»06:00″]

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

Pensamiento

«El valor de un problema no es tanto el de encontrar la respuesta como el de las ideas e intentos que obliga su resolución.»

Israel Nathan Herstein.

[/schedule]

[schedule on=’2020-03-30′ at=»06:00″]

import Data.Numbers.Primes
import Graphics.Gnuplot.Simple

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]
distribucionPrimosModulo4 = (2,0,0) : aux (tail primes) (0, 0)
  where aux (p:ps) (a,b)
          | p `mod` 4 == 1 = (p,a+1,b) : aux ps (a+1,b)
          | otherwise      = (p,a,b+1) : aux ps (a,b+1)

empatesRestosModulo4 :: [Integer]
empatesRestosModulo4 =
  [p | (p,a,b) <- distribucionPrimosModulo4
     , a == b]
  
mayoria1RestosModulo4 :: [Integer]
mayoria1RestosModulo4 =
  [p | (p,a,b) <- distribucionPrimosModulo4
     , a > b]
  
grafica :: Int -> IO ()
grafica n = 
  plotLists [Key Nothing]
            [ [(p, a) | (p,a,b) <- xs]
            , [(p, b) | (p,a,b) <- xs]
            , [(p, b-a) | (p,a,b) <- xs]]
  where xs = take n (tail (distribucionPrimosModulo4))

graficaChebyshev :: Int -> IO ()
graficaChebyshev n = 
  plotList [ Key Nothing
           , PNG "El_sesgo_de_Chebyshev.png"
           ]
           [(p, b-a) | (p,a,b) <- xs]
  where xs = take n (tail (distribucionPrimosModulo4))

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

distribucionPrimosModulo4 :: [(Integer, Integer, Integer)]

distribucionPrimosModulo4 = (2,0,0) : aux (tail primes) (0, 0)

where aux (p:ps) (a,b)

| p `mod` 4 == 1 = (p,a+1,b) : aux ps (a+1,b)

| otherwise = (p,a,b+1) : aux ps (a,b+1)

empatesRestosModulo4 :: [Integer]

empatesRestosModulo4 =

[p | (p,a,b) <- distribucionPrimosModulo4

, a == b]

mayoria1RestosModulo4 :: [Integer]

mayoria1RestosModulo4 =

[p | (p,a,b) <- distribucionPrimosModulo4

, a > b]

grafica :: Int -> IO ()

grafica n =

plotLists [Key Nothing]

[ [(p, a) | (p,a,b) <- xs]

, [(p, b) | (p,a,b) <- xs]

, [(p, b-a) | (p,a,b) <- xs]]

where xs = take n (tail (distribucionPrimosModulo4))

graficaChebyshev :: Int -> IO ()

graficaChebyshev n =

plotList [ Key Nothing

, PNG "El_sesgo_de_Chebyshev.png"

]

[(p, b-a) | (p,a,b) <- xs]

where xs = take n (tail (distribucionPrimosModulo4))

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>

[/schedule]

La conjetura de Mertens

PorJosé A. Alonso 27-enero-20203-febrero-2020

Un número entero n es libre de cuadrados si no existe un número primo p tal que p² divide a n; es decir, los factores primos de n son todos distintos.

La función de Möbius μ(n) está definida para todos los enteros positivos como sigue:

μ(n) = 1 si n es libre de cuadrados y tiene un número par de factores primos.
μ(n) = -1 si n es libre de cuadrados y tiene un número impar de factores primos.
μ(n) = 0 si n no es libre de cuadrados.

Sus primeros valores son 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, …

La función de Mertens M(n) está definida para todos los enteros positivos como la suma de μ(k) para 1 ≤ k ≤ n. Sus primeros valores son 1, 0, -1, -1, -2, -1, -2, -2, …

La conjetura de Mertens afirma que

Para todo entero x mayor que 1, el valor absoluto de la función de Mertens en x es menor que la raíz cuadrada de x.

La conjetura fue planteada por Franz Mertens en 1897. Riele Odlyzko, demostraronen 1985 que la conjetura de Mertens deja de ser cierta más o menos a partir de $10^{10^{64}}$ , cifra que luego de algunos refinamientos se redujo a $10^{10^{40}}$ .

Definir las funciones

   mobius :: Integer -> Integer
   mertens :: Integer -> Integer
   graficaMertens :: Integer -> IO ()

mobius :: Integer -> Integer

mertens :: Integer -> Integer

graficaMertens :: Integer -> IO ()

tales que

(mobius n) es el valor de la función de Möbius en n. Por ejemplo,

     mobius 6   ==  1
     mobius 30  ==  -1
     mobius 12  ==  0

mobius 6 == 1

mobius 30 == -1

mobius 12 == 0

(mertens n) es el valor de la función de Mertens en n. Por ejemplo,

     mertens 1     ==  1
     mertens 2     ==  0
     mertens 3     ==  -1
     mertens 5     ==  -2
     mertens 661   ==  -11
     mertens 1403  ==  11

mertens 1 == 1

mertens 2 == 0

mertens 3 == -1

mertens 5 == -2

mertens 661 == -11

mertens 1403 == 11

(graficaMertens n) dibuja la gráfica de la función de Mertens, la raíz cuadrada y el opuestos de la raíz cuadrada para los n primeros n enteros positivos. Por ejemplo, (graficaMertens 1000) dibuja

Comprobar con QuickCheck la conjetura de Mertens.

Nota: El ejercicio está basado en La conjetura de Merterns y su relación con un número tan raro como extremada y colosalmente grande publicado por @Alvy la semana pasada en Microsiervos.

Soluciones

import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck
import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer
mobius n | tieneRepetidos xs = 0
         | otherwise         = (-1)^(length xs)
  where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool
tieneRepetidos xs =
  or [x == y | (x,y) <- zip xs (tail xs)]
        
mertens :: Integer -> Integer
mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens
-- ============================

graficaMertens :: Integer -> IO ()
graficaMertens n = do
  plotLists [ Key Nothing
            , Title "Conjetura de Mertens"
            , PNG "La_conjetura_de_Mertens.png"
            ]
            [ [mertens k | k <- [1..n]]
            , raices
            , map negate raices
            ]
    
  where
    raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens
-- ====================

-- La conjetura es
conjeturaDeMertens :: Integer -> Property
conjeturaDeMertens n =
  n > 1
  ==>
  abs (mertens n) < ceiling (sqrt n')
  where n' = fromIntegral n

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer

mobius n | tieneRepetidos xs = 0

| otherwise = (-1)^(length xs)

where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool

tieneRepetidos xs =

or [x == y | (x,y) <- zip xs (tail xs)]

mertens :: Integer -> Integer

mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens

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

graficaMertens :: Integer -> IO ()

graficaMertens n = do

plotLists [ Key Nothing

, Title "Conjetura de Mertens"

, PNG "La_conjetura_de_Mertens.png"

]

[ [mertens k | k <- [1..n]]

, raices

, map negate raices

]

where

raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens

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

-- La conjetura es

conjeturaDeMertens :: Integer -> Property

conjeturaDeMertens n =

n > 1

==>

abs (mertens n) < ceiling (sqrt n')

where n' = fromIntegral n

-- La comprobación es

-- λ> quickCheck conjeturaDeMertens

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

«El control de la complejidad es la esencia de la programación informática.»

Brian Kernighan.

Distribución de diferencias de dígitos consecutivos de pi

PorJosé A. Alonso 4-abril-201911-abril-2019

Usando la librería Data.Number.CReal, que se instala con

   cabal install number

1	cabal install number

se pueden calcular el número pi con la precisión que se desee. Por ejemplo,

   λ> import Data.Number.CReal
   λ> showCReal 60 pi
   "3.141592653589793238462643383279502884197169399375105820974945"

λ> import Data.Number.CReal

λ> showCReal 60 pi

"3.141592653589793238462643383279502884197169399375105820974945"

importa la librería y calcula el número pi con 60 decimales.

La distribución de las diferencias de los dígitos consecutivos para los 18 primeros n dígitos de pi se calcula como sigue: los primeros 18 dígitos de pi son

   3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

1	3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

Las diferencias de sus elementos consecutivos es

   2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

1	2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

y la distribución de sus frecuencias en el intervalo [-9,9] es

   0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

1	0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

es decir, el desde el -9 a -5 no aparecen, el -4 aparece 3 veces, el -2 aparece 2 veces y así sucesivamente.

Definir las funciones

   distribucionDDCpi :: Int -> [Int]
   graficas :: [Int] -> FilePath -> IO ()

1 2	distribucionDDCpi :: Int -> [Int] graficas :: [Int] -> FilePath -> IO ()

tales que

(distribucionDDCpi n) es la distribución de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi. Por ejemplo,

     λ> distribucionDDCpi 18
     [0,0,0,0,0,3,2,2,2,0,2,3,1,0,0,1,1,0,0]
     λ> distribucionDDCpi 100
     [1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]
     λ> distribucionDDCpi 200
     [3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]
     λ> distribucionDDCpi 1000
     [16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]
     λ> distribucionDDCpi 5000
     [67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

λ> distribucionDDCpi 18

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

λ> distribucionDDCpi 100

[1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]

λ> distribucionDDCpi 200

[3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]

λ> distribucionDDCpi 1000

[16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]

λ> distribucionDDCpi 5000

[67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

(graficas ns f) dibuja en el fichero f las gráficas de las distribuciones de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi, para n en ns. Por ejemplo, al evaluar (graficas [100,250..4000] «distribucionDDCpi.png» se escribe en el fichero «distribucionDDCpi.png» la siguiente gráfica

Soluciones

import Data.Number.CReal
import Graphics.Gnuplot.Simple
import Data.Array

--    λ> digitosPi 18
--    [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]
digitosPi :: Int -> [Int]
digitosPi n = init [read [c] | c <- (x:xs)]
  where (x:_:xs) = showCReal n pi

--    λ> diferenciasConsecutivos (digitosPi 18)
--    [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]
diferenciasConsecutivos :: Num a => [a] -> [a]
diferenciasConsecutivos xs =
  zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]
distribucionDDCpi =
  distribucion . diferenciasConsecutivos . digitosPi
  where distribucion xs =
          elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()
graficas ns f = 
  plotLists [Key Nothing, PNG f]
            [puntos n | n <- ns]
  where puntos :: Int -> [(Int,Int)]
        puntos n = zip [-9..9] (distribucionDDCpi n)

import Data.Number.CReal

import Graphics.Gnuplot.Simple

import Data.Array

-- λ> digitosPi 18

-- [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]

digitosPi :: Int -> [Int]

digitosPi n = init [read [c] | c <- (x:xs)]

where (x:_:xs) = showCReal n pi

-- λ> diferenciasConsecutivos (digitosPi 18)

-- [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]

diferenciasConsecutivos :: Num a => [a] -> [a]

diferenciasConsecutivos xs =

zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]

distribucionDDCpi =

distribucion . diferenciasConsecutivos . digitosPi

where distribucion xs =

elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()

graficas ns f =

plotLists [Key Nothing, PNG f]

[puntos n | n <- ns]

where puntos :: Int -> [(Int,Int)]

puntos n = zip [-9..9] (distribucionDDCpi n)

Pensamiento

Doy consejo, a fuer de viejo:
nunca sigas mi consejo.

Antonio Machado

Números apocalípticos

PorJosé A. Alonso 8-enero-201819-febrero-2018

Un número apocalíptico es aquel número natural n tal que 2^n contiene la secuencia 666.

Definir las funciones

   esApocaliptico           :: Integer -> Bool
   apocalipticos            :: [Integer]
   mayorNoApocalipticoMenor :: Integer -> Maybe Integer
   grafica                  :: Integer -> IO ()

esApocaliptico :: Integer -> Bool

apocalipticos :: [Integer]

mayorNoApocalipticoMenor :: Integer -> Maybe Integer

grafica :: Integer -> IO ()

tales que

(esApocaliptico n) se verifica si n es un número apocalíptico. Por ejemplo,

     esApocaliptico 666    ==  True
     esApocaliptico 29784  ==  False

1 2	esApocaliptico 666 == True esApocaliptico 29784 == False

apocalipticos es la lista de los números apocalípticos. Por ejemplo,

     take 9 apocalipticos  ==  [157,192,218,220,222,224,226,243,245]
     apocalipticos !! 55   ==  666

1 2	take 9 apocalipticos == [157,192,218,220,222,224,226,243,245] apocalipticos !! 55 == 666

(mayorNoApocalipticoMenor n) es justo el mayor número no apocalíptico menor que n. Por ejemplo,

     mayorNoApocalipticoMenor  40000  ==  Just 29784
     mayorNoApocalipticoMenor  29784  ==  Just 26667

1 2	mayorNoApocalipticoMenor 40000 == Just 29784 mayorNoApocalipticoMenor 29784 == Just 26667

(grafica n) dibuja las gráficas de los n primeros términos de la sucesión de los números apocalípticos junto con los de la sucesión a(n) = 3715+n. Por ejemplo, (grafica 3000) dibuja

y (grafica 30000) dibuja

Nota: Este ejercicio ha sido propuesto por Ángel Ruiz Campos.

Soluciones

import Data.List (isInfixOf, find, genericTake)
import Graphics.Gnuplot.Simple

esApocaliptico :: Integer -> Bool
esApocaliptico = isInfixOf "666" . show . (2^)

apocalipticos :: [Integer]
apocalipticos = filter esApocaliptico [1..]

mayorNoApocalipticoMenor :: Integer -> Maybe Integer
mayorNoApocalipticoMenor n = find (not . esApocaliptico) [n-1,n-2..1]

grafica :: Integer -> IO ()
grafica n = 
  plotLists [ Key Nothing
            , PNG ("Numeros_apocalipticos_" ++ show n ++ ".png")
            ]
            [ genericTake n apocalipticos
            , [3715..3715+n-1] ]

import Data.List (isInfixOf, find, genericTake)

import Graphics.Gnuplot.Simple

esApocaliptico :: Integer -> Bool

esApocaliptico = isInfixOf "666" . show . (2^)

apocalipticos :: [Integer]

apocalipticos = filter esApocaliptico [1..]

mayorNoApocalipticoMenor :: Integer -> Maybe Integer

mayorNoApocalipticoMenor n = find (not . esApocaliptico) [n-1,n-2..1]

grafica :: Integer -> IO ()

grafica n =

plotLists [ Key Nothing

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

]

[ genericTake n apocalipticos

, [3715..3715+n-1] ]

Distancias entre primos consecutivos

PorJosé A. Alonso 4-abril-201711-abril-2017

Los 15 primeros números primos son

   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47

1	2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47

Las distancias entre los elementos consecutivos son

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

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

La distribución de las distancias es

   (1,1), (2,6), (4,5), (6,2)

1	(1,1), (2,6), (4,5), (6,2)

(es decir, el 1 aparece una vez, el 2 aparece 6 veces, etc.) La frecuencia de las distancias es

   (1,7.142857), (2,42.857143), (4,35.714287), (6,14.285714)

1	(1,7.142857), (2,42.857143), (4,35.714287), (6,14.285714)

(es decir, el 1 aparece el 7.142857%, el 2 el 42.857143% etc.)

Definir las funciones

  cuentaDistancias        :: Int -> [(Int,Int)]
  frecuenciasDistancias   :: Int -> [(Int,Float)]
  graficas                :: [Int] -> IO ()
  distanciasMasFrecuentes :: Int -> [Int]

cuentaDistancias :: Int -> [(Int,Int)]

frecuenciasDistancias :: Int -> [(Int,Float)]

graficas :: [Int] -> IO ()

distanciasMasFrecuentes :: Int -> [Int]

tales que

(cuentaDistancias n) es la distribución de distancias entre los n primeros primos consecutivos. Por ejemplo,

     λ> cuentaDistancias 15
     [(1,1),(2,6),(4,5),(6,2)]
     λ> cuentaDistancias 100
     [(1,1),(2,25),(4,26),(6,25),(8,7),(10,7),(12,4),(14,3),(18,1)]

λ> cuentaDistancias 15

[(1,1),(2,6),(4,5),(6,2)]

λ> cuentaDistancias 100

[(1,1),(2,25),(4,26),(6,25),(8,7),(10,7),(12,4),(14,3),(18,1)]

(frecuenciasDistancias n) es la frecuencia de distancias entre los n primeros primos consecutivos. Por ejemplo,

     λ> frecuenciasDistancias 15
     [(1,7.142857),(2,42.857143),(4,35.714287),(6,14.285714)]
     λ> frecuenciasDistancias 30
     [(1,3.4482758),(2,34.482758),(4,34.482758),(6,24.137932),(8,3.4482758)]

λ> frecuenciasDistancias 15

[(1,7.142857),(2,42.857143),(4,35.714287),(6,14.285714)]

λ> frecuenciasDistancias 30

[(1,3.4482758),(2,34.482758),(4,34.482758),(6,24.137932),(8,3.4482758)]

(graficas ns) dibuja las gráficas de (frecuenciasDistancias k) para k en ns. Por ejemplo, (graficas [10,20,30]) dibuja

(graficas [1000,2000,3000]) dibuja

y (graficas [100000,200000,300000]) dibuja
(distanciasMasFrecuentes n) es la lista de las distancias más frecuentes entre los elementos consecutivos de la lista de los n primeros primos. Por ejemplo,

     distanciasMasFrecuentes 15     ==  [2]
     distanciasMasFrecuentes 26     ==  [2,4]
     distanciasMasFrecuentes 32     ==  [4]
     distanciasMasFrecuentes 41     ==  [2,4,6]
     distanciasMasFrecuentes 77     ==  [6]
     distanciasMasFrecuentes 160    ==  [4,6]
     distanciasMasFrecuentes (10^6) == [6]

distanciasMasFrecuentes 15 == [2]

distanciasMasFrecuentes 26 == [2,4]

distanciasMasFrecuentes 32 == [4]

distanciasMasFrecuentes 41 == [2,4,6]

distanciasMasFrecuentes 77 == [6]

distanciasMasFrecuentes 160 == [4,6]

distanciasMasFrecuentes (10^6) == [6]

Comprobar con QuickCheck si para todo n > 160 se verifica que (distanciasMasFrecuentes n) es [6].

Soluciones

import Data.Numbers.Primes
import qualified Data.Map as M 
import Graphics.Gnuplot.Simple
import Test.QuickCheck

cuentaDistancias :: Int -> [(Int,Int)]
cuentaDistancias = M.toList . dicDistancias

dicDistancias :: Int -> M.Map Int Int
dicDistancias n =
  M.fromListWith (+) (zip (distancias n) (repeat 1))

distancias :: Int -> [Int]
distancias n =
  zipWith (-) (tail xs) xs
  where xs = take n primes

frecuenciasDistancias :: Int -> [(Int,Float)]
frecuenciasDistancias n =
  [(k,(100 * fromIntegral x) / n1) | (k,x) <- cuentaDistancias n]
  where n1 = fromIntegral (n-1)

graficas :: [Int] -> IO ()
graficas ns =
  plotLists [Key Nothing]
            (map frecuenciasDistancias ns)

distanciasMasFrecuentes :: Int -> [Int]
distanciasMasFrecuentes n =
  M.keys (M.filter (==m) d)
  where d = dicDistancias n
        m = maximum (M.elems d)

-- La propiedad es
prop_distanciasMasFrecuentes :: Int -> Bool
prop_distanciasMasFrecuentes n =
  distanciasMasFrecuentes (161 + abs n) == [6]

import Data.Numbers.Primes

import qualified Data.Map as M

import Graphics.Gnuplot.Simple

import Test.QuickCheck

cuentaDistancias :: Int -> [(Int,Int)]

cuentaDistancias = M.toList . dicDistancias

dicDistancias :: Int -> M.Map Int Int

dicDistancias n =

M.fromListWith (+) (zip (distancias n) (repeat 1))

distancias :: Int -> [Int]

distancias n =

zipWith (-) (tail xs) xs

where xs = take n primes

frecuenciasDistancias :: Int -> [(Int,Float)]

frecuenciasDistancias n =

[(k,(100 * fromIntegral x) / n1) | (k,x) <- cuentaDistancias n]

where n1 = fromIntegral (n-1)

graficas :: [Int] -> IO ()

graficas ns =

plotLists [Key Nothing]

(map frecuenciasDistancias ns)

distanciasMasFrecuentes :: Int -> [Int]

distanciasMasFrecuentes n =

M.keys (M.filter (==m) d)

where d = dicDistancias n

m = maximum (M.elems d)

-- La propiedad es

prop_distanciasMasFrecuentes :: Int -> Bool

prop_distanciasMasFrecuentes n =

distanciasMasFrecuentes (161 + abs n) == [6]

Distribución de diferencias de dígitos consecutivos de pi

PorJosé A. Alonso 14-marzo-201725-abril-2021

La distribución de las diferencias de los dígitos consecutivos para los 18 primeros dígitos de pi se calcula como sigue: los primeros 18 dígitos de pi son

   3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

1	3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

Las diferencias de sus elementos consecutivos es

   2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

1	2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

y la distribución de sus frecuencias en el intervalo [-9,9] es

   0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

1	0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

es decir, el desde el -9 a -5 no aparecen, el -4 aparece 3 veces, el -2 aparece 2 veces y así sucesivamente.

Definir las funciones

   distribucionDDCpi :: Int -> [Int]
   graficas :: [Int] -> FilePath -> IO ()

1 2	distribucionDDCpi :: Int -> [Int] graficas :: [Int] -> FilePath -> IO ()

tales que

(distribucionDDCpi n) es la distribución de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi. Por ejemplo,

   λ> distribucionDDCpi 18
   [0,0,0,0,0,3,2,2,2,0,2,3,1,0,0,1,1,0,0]
   λ> distribucionDDCpi 100
   [1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]
   λ> distribucionDDCpi 200
   [3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]
   λ> distribucionDDCpi 1000
   [16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]
   λ> distribucionDDCpi 5000
   [67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

λ> distribucionDDCpi 18

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

λ> distribucionDDCpi 100

[1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]

λ> distribucionDDCpi 200

[3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]

λ> distribucionDDCpi 1000

[16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]

λ> distribucionDDCpi 5000

[67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

(graficas ns f) dibuja en el fichero f las gráficas de las distribuciones de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi, para n en ns. Por ejemplo, al evaluar (graficas [100,250..4000] «distribucionDDCpi.png» se escribe en el fichero «distribucionDDCpi.png» la siguiente gráfica

Nota: Se puede usar la librería Data.Number.CReal.

Soluciones

import Data.Number.CReal
import Graphics.Gnuplot.Simple
import Data.Array

--    λ> digitosPi 18
--    [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]
digitosPi :: Int -> [Int]
digitosPi n = init [read [c] | c <- (x:xs)]
  where (x:_:xs) = showCReal n pi

--    λ> diferenciasConsecutivos (digitosPi 18)
--    [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]
diferenciasConsecutivos :: Num a => [a] -> [a]
diferenciasConsecutivos xs =
  zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]
distribucionDDCpi =
  distribucion . diferenciasConsecutivos . digitosPi
  where distribucion xs =
          elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()
graficas ns f = 
  plotLists [Key Nothing, PNG f]
            [puntos n | n <- ns]
  where puntos :: Int -> [(Int,Int)]
        puntos n = zip [-9..9] (distribucionDDCpi n)

import Data.Number.CReal

import Graphics.Gnuplot.Simple

import Data.Array

-- λ> digitosPi 18

-- [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]

digitosPi :: Int -> [Int]

digitosPi n = init [read [c] | c <- (x:xs)]

where (x:_:xs) = showCReal n pi

-- λ> diferenciasConsecutivos (digitosPi 18)

-- [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]

diferenciasConsecutivos :: Num a => [a] -> [a]

diferenciasConsecutivos xs =

zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]

distribucionDDCpi =

distribucion . diferenciasConsecutivos . digitosPi

where distribucion xs =

elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()

graficas ns f =

plotLists [Key Nothing, PNG f]

[puntos n | n <- ns]

where puntos :: Int -> [(Int,Int)]

puntos n = zip [-9..9] (distribucionDDCpi n)

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

PorJosé A. Alonso 23-febrero-201726-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.

Nota: Este ejercicio ha sido propuesto por Enrique Naranjo.

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

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

Soluciones

Otras soluciones

El sesgo de Chebyshev

Soluciones

Pensamiento

Otras soluciones

La conjetura de Mertens

Soluciones

Otras soluciones

Pensamiento

Distribución de diferencias de dígitos consecutivos de pi

Soluciones

Pensamiento

Números apocalípticos

Soluciones

Distancias entre primos consecutivos

Soluciones

Distribución de diferencias de dígitos consecutivos de pi

Soluciones

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

Soluciones

Entradas recientes

Buscar

Correo electrónico