pi – Exercitium

Cálculo de pi con el producto de Wallis

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

El producto de Wallis es una expresión, descubierta por John Wallis en 1655, para representar el valor de π y que establece que:

    π     2     2     4     4     6     6     8     8
   --- = --- · --- · --- · --- · --- · --- · --- · --- ···
    2     1     3     3     5     5     7     7     9

π 2 2 4 4 6 6 8 8

--- = --- · --- · --- · --- · --- · --- · --- · --- ···

2 1 3 3 5 5 7 7 9

Definir las funciones

   factoresWallis  :: [Rational]
   productosWallis :: [Rational]
   aproximacionPi  :: Int -> Double
   errorPi         :: Double -> Int

factoresWallis :: [Rational]

productosWallis :: [Rational]

aproximacionPi :: Int -> Double

errorPi :: Double -> Int

tales que

factoresWallis es la sucesión de los factores del productos de Wallis. Por ejemplo,

     λ> take 10 factoresWallis
     [2 % 1,2 % 3,4 % 3,4 % 5,6 % 5,6 % 7,8 % 7,8 % 9,10 % 9,10 % 11]

1 2	λ> take 10 factoresWallis [2 % 1,2 % 3,4 % 3,4 % 5,6 % 5,6 % 7,8 % 7,8 % 9,10 % 9,10 % 11]

productosWallis es la sucesión de los productos de los primeros factores de Wallis. Por ejemplo,

     λ> take 7 productosWallis
     [2 % 1,4 % 3,16 % 9,64 % 45,128 % 75,256 % 175,2048 % 1225]

1 2	λ> take 7 productosWallis [2 % 1,4 % 3,16 % 9,64 % 45,128 % 75,256 % 175,2048 % 1225]

(aproximacionPi n) es la aproximación de pi obtenida multiplicando los n primeros factores de Wallis. Por ejemplo,

     aproximacionPi 20     ==  3.2137849402931895
     aproximacionPi 200    ==  3.1493784731686008
     aproximacionPi 2000   ==  3.142377365093878
     aproximacionPi 20000  ==  3.141671186534396

aproximacionPi 20 == 3.2137849402931895

aproximacionPi 200 == 3.1493784731686008

aproximacionPi 2000 == 3.142377365093878

aproximacionPi 20000 == 3.141671186534396

(errorPi x) es el menor número de factores de Wallis necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1     ==  14
     errorPi 0.01    ==  155
     errorPi 0.001   ==  1569
     errorPi 0.0001  ==  15707

errorPi 0.1 == 14

errorPi 0.01 == 155

errorPi 0.001 == 1569

errorPi 0.0001 == 15707

Soluciones

import Data.Ratio

factoresWallis :: [Rational]
factoresWallis =
  concat [[y%(y-1),  y%(y+1)] | x <- [1..], let y = 2*x]

productosWallis :: [Rational]
productosWallis = scanl1 (*) factoresWallis

aproximacionPi :: Int -> Double
aproximacionPi n =
  fromRational (2 * productosWallis !! n)

errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x =
  length (takeWhile (>=x) [abs (pi - 2 * fromRational y)
                          | y <- productosWallis])

-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n =
  2 * productosWallis2 !! n

productosWallis2 :: [Double]
productosWallis2 = scanl1 (*) factoresWallis2

factoresWallis2 :: [Double]
factoresWallis2 =
  concat [[y/(y-1),  y/(y+1)] | x <- [1..], let y = 2*x]

-- 3ª definición de errorPi
errorPi3 :: Double -> Int
errorPi3 x = head [n | n <- [1..]
                     , abs (pi - aproximacionPi2 n) < x]

-- Comparación de eficiencia
--    λ> errorPi 0.001
--    1569
--    (0.82 secs, 374,495,816 bytes)
--
--    λ> errorPi2 0.001
--    1569
--    (0.79 secs, 369,282,320 bytes)
--
--    λ> errorPi3 0.001
--    1569
--    (0.04 secs, 0 bytes)

import Data.Ratio

factoresWallis :: [Rational]

factoresWallis =

concat [[y%(y-1), y%(y+1)] | x <- [1..], let y = 2*x]

productosWallis :: [Rational]

productosWallis = scanl1 (*) factoresWallis

aproximacionPi :: Int -> Double

aproximacionPi n =

fromRational (2 * productosWallis !! n)

errorPi :: Double -> Int

errorPi x = head [n | n <- [1..]

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x =

length (takeWhile (>=x) [abs (pi - 2 * fromRational y)

| y <- productosWallis])

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

aproximacionPi2 n =

2 * productosWallis2 !! n

productosWallis2 :: [Double]

productosWallis2 = scanl1 (*) factoresWallis2

factoresWallis2 :: [Double]

factoresWallis2 =

concat [[y/(y-1), y/(y+1)] | x <- [1..], let y = 2*x]

-- 3ª definición de errorPi

errorPi3 :: Double -> Int

errorPi3 x = head [n | n <- [1..]

, abs (pi - aproximacionPi2 n) < x]

-- Comparación de eficiencia

-- λ> errorPi 0.001

-- 1569

-- (0.82 secs, 374,495,816 bytes)

-- λ> errorPi2 0.001

-- 1569

-- (0.79 secs, 369,282,320 bytes)

-- λ> errorPi3 0.001

-- 1569

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

Pensamiento

«¿Por qué son hermosos los números? Es como preguntar por qué es bella la Novena Sinfonía de Beethoven. Si no ves por qué, alguien no puede decírtelo. Yo sé que los números son hermosos. Si no son hermosos, nada lo es.»

Paul Erdös.

Cálculo de pi usando la fórmula de Vieta

PorJosé A. Alonso 16-diciembre-201923-diciembre-2019

La fórmula de Vieta para el cálculo de pi es la siguiente

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

1 2	aproximacionPi :: Int -> Double errorPi :: Double -> Int

tales que

(aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,

     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795

aproximacionPi 5 == 3.140331156954753

aproximacionPi 10 == 3.1415914215112

aproximacionPi 15 == 3.141592652386592

aproximacionPi 20 == 3.1415926535886207

aproximacionPi 25 == 3.141592653589795

(errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

errorPi 0.1 == 2

errorPi 0.01 == 4

errorPi 0.001 == 6

errorPi 0.0001 == 7

errorPi 1e-4 == 7

errorPi 1e-14 == 24

pi == 3.141592653589793

aproximacionPi 24 == 3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

-- 1ª definición de aproximacionPi

aproximacionPi :: Int -> Double

aproximacionPi n = product [2 / aux x | x <- [0..n]]

where

aux 0 = 1

aux 1 = sqrt 2

aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = product [2/x | x <- 1 : xs]

where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi

aproximacionPi3 :: Int -> Double

aproximacionPi3 n = product (2 : take n (map (2/) xs))

where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi

errorPi :: Double -> Int

errorPi x = head [n | n <- [1..]

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x = until aceptable (+1) 1

where aceptable n = abs (pi - aproximacionPi n) < x

Pensamiento

El tiempo que la barba me platea,
cavó mis ojos y agrandó mi frente,
va siendo en mi recuerdo transparente,
y mientras más al fondo, más clarea.

Antonio Machado

Aproximación entre pi y e

PorJosé A. Alonso 7-diciembre-201817-enero-2019

El día 11 de noviembre, se publicó en la cuenta de Twitter de Fermat’s Library la siguiente curiosa identidad que relaciona los números e y pi:

Definir las siguientes funciones:

   sumaTerminos :: Int -> Double
   aproximacion :: Double -> Int

1 2	sumaTerminos :: Int -> Double aproximacion :: Double -> Int

tales que

(sumaTerminos n) es la suma de los primeros n términos de la serie 1/(π²+ 1) + 1/(4π²+1) + 1/(9π²+1) + 1/(16π²+ ) + … Por ejemplo,

     sumaTerminos 10     ==  0.14687821811081034
     sumaTerminos 100    ==  0.15550948345688423
     sumaTerminos 1000   ==  0.15641637221314514
     sumaTerminos 10000  ==  0.15650751113789382

sumaTerminos 10 == 0.14687821811081034

sumaTerminos 100 == 0.15550948345688423

sumaTerminos 1000 == 0.15641637221314514

sumaTerminos 10000 == 0.15650751113789382

(aproximación x) es el menor número de términos que hay que sumar de la serie anterior para que se diferencie (en valor absoluto) de 1/(e²-1) menos que x. Por ejemplo,

     aproximacion 0.1     ==  1
     aproximacion 0.01    ==  10
     aproximacion 0.001   ==  101
     aproximacion 0.0001  ==  1013

aproximacion 0.1 == 1

aproximacion 0.01 == 10

aproximacion 0.001 == 101

aproximacion 0.0001 == 1013

Soluciones

-- 1ª definición de sumaTerminos
sumaTerminos :: Int -> Double
sumaTerminos n =
  sum [1 / (((x ^ 2) * (pi ^ 2)) + 1) | x <- [1 .. fromIntegral n]]

-- 2ª definición de sumaTerminos
sumaTerminos2 :: Int -> Double
sumaTerminos2 0 = 0
sumaTerminos2 n = 1 / (m^2 * pi^2 + 1) + sumaTerminos2 (n-1)
  where m = fromIntegral n

-- Definición de aproximacion
aproximacion :: Double -> Int
aproximacion x =
  head [n | n <- [0..]
          , abs (sumaTerminos n - 1 / (e^2 - 1)) < x]
  where e = exp 1

-- 1ª definición de sumaTerminos

sumaTerminos :: Int -> Double

sumaTerminos n =

sum [1 / (((x ^ 2) * (pi ^ 2)) + 1) | x <- [1 .. fromIntegral n]]

-- 2ª definición de sumaTerminos

sumaTerminos2 :: Int -> Double

sumaTerminos2 0 = 0

sumaTerminos2 n = 1 / (m^2 * pi^2 + 1) + sumaTerminos2 (n-1)

where m = fromIntegral n

-- Definición de aproximacion

aproximacion :: Double -> Int

aproximacion x =

head [n | n <- [0..]

, abs (sumaTerminos n - 1 / (e^2 - 1)) < x]

where e = exp 1

Pensamiento

«Sólo sé que no se nada» contenía la jactancia de un excesivo saber, puesto que olvidó añadir: y aun de esto mismo no estoy completamente seguro.

Antonio Machado

Fractal hexagonal

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

Escribir, usando CodeWorld, un programa para dibujar el fractal hexagonal que se muestra en la siguiente animación

Las 4 primeras fases de la animación son

Fase 0:
Fase 1:
Fase 2:
Fase 3:

Nota: Este ejercicio ha sido propuesto por Agustín Martín Aguera.

Soluciones

import CodeWorld

main :: IO()
main = animationOf (hexagono . s)

hexagono :: Int -> Picture
hexagono 0 =
  colored red $ solidPolygon [(9,0),(4.5,c),(-4.5,c),(-9,0),(-4.5,-c),(4.5,-c)]
  where c = 9 * sin (pi / 3)
hexagono n =
  pictures (hex : take 6 (iterate (rotated (pi / 3)) (translated 6 0 hex)))
  where hex = scaled (1/3) (1/3) $ hexagono (n-1)

s :: Double -> Int
s t = mod (floor t) 5

import CodeWorld

main :: IO()

main = animationOf (hexagono . s)

hexagono :: Int -> Picture

hexagono 0 =

colored red $ solidPolygon [(9,0),(4.5,c),(-4.5,c),(-9,0),(-4.5,-c),(4.5,-c)]

where c = 9 * sin (pi / 3)

hexagono n =

pictures (hex : take 6 (iterate (rotated (pi / 3)) (translated 6 0 hex)))

where hex = scaled (1/3) (1/3) $ hexagono (n-1)

s :: Double -> Int

s t = mod (floor t) 5

Número de dígitos del factorial

PorJosé A. Alonso 23-marzo-201731-marzo-2017

Definir las funciones

   nDigitosFact :: Integer -> Integer
   graficas     :: [Integer] -> IO ()

1 2	nDigitosFact :: Integer -> Integer graficas :: [Integer] -> IO ()

tales que

(nDigitosFact n) es el número de dígitos de n!. Por ejemplo,

     nDigitosFact 0        ==  1
     nDigitosFact 4        ==  2
     nDigitosFact 5        ==  3
     nDigitosFact 10       ==  7
     nDigitosFact 100      ==  158
     nDigitosFact 1000     ==  2568
     nDigitosFact 10000    ==  35660
     nDigitosFact 100000   ==  456574
     nDigitosFact 1000000  ==  5565709

nDigitosFact 0 == 1

nDigitosFact 4 == 2

nDigitosFact 5 == 3

nDigitosFact 10 == 7

nDigitosFact 100 == 158

nDigitosFact 1000 == 2568

nDigitosFact 10000 == 35660

nDigitosFact 100000 == 456574

nDigitosFact 1000000 == 5565709

(graficas xs) dibuja las gráficas de los números de dígitos del factorial de k (para k en xs) y de la recta y = 5.5 x. Por ejemplo, (graficas [0,500..10^6]) dibuja

Nota: Este ejercicio está basado en el problema How many digits? de Kattis en donde se impone la restricción de calcular, en menos de 1 segundo, el número de dígitos de los factoriales de 10.000 números del rango [0,1.000.000].

Se puede simular como sigue

   λ> import System.Random 
   λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]
   λ> maximum (map nDigitosFact xs)
   5561492
   λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]
   λ> maximum (map nDigitosFact xs)
   5563303

λ> import System.Random

λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]

λ> maximum (map nDigitosFact xs)

5561492

λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]

λ> maximum (map nDigitosFact xs)

5563303

Soluciones

import Data.List (genericLength, genericIndex)
import Graphics.Gnuplot.Simple

-- 1ª definición
nDigitosFact1 :: Integer -> Integer
nDigitosFact1 n =
  genericLength (show (product [1..n]))

-- 2ª definición (usando logaritmos)
-- =================================

-- Basado en
--    nDigitos (n!) = 1 + floor (log (n!))
--                  = 1 + floor (log (1*2*3*...*n))
--                  = 1 + floor (log(1) + log(2) + log(3) +...+ log(n))

nDigitosFact2 :: Integer -> Integer
nDigitosFact2 n =
  1 + floor (sum [logBase 10 (fromIntegral k) | k <- [1..n]])

-- 3ª definición (usando logaritmos y programación dinámica)
-- =========================================================

nDigitosFact3 :: Integer -> Integer
nDigitosFact3 0 = 1
nDigitosFact3 n = 1 + floor ((sumLogs `genericIndex` (n-1)) / log 10)

logs :: [Double]
logs = map log [1..]

sumLogs :: [Double]
sumLogs = scanl1 (+) logs

-- 4ª definición (Usando la fórmula de Kamenetsky)
-- ===============================================

-- La fórmula se encuentra en https://oeis.org/A034886

nDigitosFact4 :: Integer -> Integer
nDigitosFact4 0 = 1
nDigitosFact4 1 = 1
nDigitosFact4 n =
  1 + floor (m * logBase 10 (m/e) + logBase 10 (2*pi*m)/2)
  where m = fromIntegral n
        e = exp 1

--    λ> nDigitosFact1 (4*10^4)
--    166714
--    (5.61 secs, 1,649,912,864 bytes)
--    λ> nDigitosFact2 (4*10^4)
--    166714
--    (0.10 secs, 13,741,360 bytes)
--
--    λ> nDigitosFact2 (7*10^5)
--    3787566
--    (1.86 secs, 187,666,328 bytes)
--    λ> nDigitosFact3 (7*10^5)
--    3787566
--    (0.02 secs, 0 bytes)
--    λ> nDigitosFact3 (7*10^5)
--    3787566
--    (1.01 secs, 238,682,064 bytes)
--    λ> nDigitosFact4 (7*10^5)
--    3787566
--    (0.01 secs, 0 bytes)
-- 
--    λ> import System.Random 
--    λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^2]]
--    λ> maximum (map nDigitosFact3 xs)
--    5565097
--    (11.19 secs, 7,336,595,440 bytes)
--    λ> maximum (map nDigitosFact4 xs)
--    5565097
--    (0.01 secs, 0 bytes)

graficas :: [Integer] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
              [p1, p2]
    where p1, p2 :: [(Float, Float)]
          p1 = [(fi k, fi (nDigitosFact4 k)) | k <- xs]
          p2 = [(k',5.5*k') | k <- xs, let k' = fi k]
          fi = fromIntegral

import Data.List (genericLength, genericIndex)

import Graphics.Gnuplot.Simple

-- 1ª definición

nDigitosFact1 :: Integer -> Integer

nDigitosFact1 n =

genericLength (show (product [1..n]))

-- 2ª definición (usando logaritmos)

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

-- Basado en

-- nDigitos (n!) = 1 + floor (log (n!))

-- = 1 + floor (log (1*2*3*...*n))

-- = 1 + floor (log(1) + log(2) + log(3) +...+ log(n))

nDigitosFact2 :: Integer -> Integer

nDigitosFact2 n =

1 + floor (sum [logBase 10 (fromIntegral k) | k <- [1..n]])

-- 3ª definición (usando logaritmos y programación dinámica)

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

nDigitosFact3 :: Integer -> Integer

nDigitosFact3 0 = 1

nDigitosFact3 n = 1 + floor ((sumLogs `genericIndex` (n-1)) / log 10)

logs :: [Double]

logs = map log [1..]

sumLogs :: [Double]

sumLogs = scanl1 (+) logs

-- 4ª definición (Usando la fórmula de Kamenetsky)

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

-- La fórmula se encuentra en https://oeis.org/A034886

nDigitosFact4 :: Integer -> Integer

nDigitosFact4 0 = 1

nDigitosFact4 1 = 1

nDigitosFact4 n =

1 + floor (m * logBase 10 (m/e) + logBase 10 (2*pi*m)/2)

where m = fromIntegral n

e = exp 1

-- λ> nDigitosFact1 (4*10^4)

-- 166714

-- (5.61 secs, 1,649,912,864 bytes)

-- λ> nDigitosFact2 (4*10^4)

-- 166714

-- (0.10 secs, 13,741,360 bytes)

-- λ> nDigitosFact2 (7*10^5)

-- 3787566

-- (1.86 secs, 187,666,328 bytes)

-- λ> nDigitosFact3 (7*10^5)

-- 3787566

-- (0.02 secs, 0 bytes)

-- λ> nDigitosFact3 (7*10^5)

-- 3787566

-- (1.01 secs, 238,682,064 bytes)

-- λ> nDigitosFact4 (7*10^5)

-- 3787566

-- (0.01 secs, 0 bytes)

-- λ> import System.Random

-- λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^2]]

-- λ> maximum (map nDigitosFact3 xs)

-- 5565097

-- (11.19 secs, 7,336,595,440 bytes)

-- λ> maximum (map nDigitosFact4 xs)

-- 5565097

-- (0.01 secs, 0 bytes)

graficas :: [Integer] -> IO ()

graficas xs =

plotLists [Key Nothing]

[p1, p2]

where p1, p2 :: [(Float, Float)]

p1 = [(fi k, fi (nDigitosFact4 k)) | k <- xs]

p2 = [(k',5.5*k') | k <- xs, let k' = fi k]

fi = fromIntegral

pi

Cálculo de pi con el producto de Wallis

Soluciones

Otras soluciones

Pensamiento

Cálculo de pi usando la fórmula de Vieta

Soluciones

Pensamiento

Aproximación entre pi y e

Soluciones

Pensamiento

Fractal hexagonal

Soluciones

Número de dígitos del factorial

Soluciones

Entradas recientes

Buscar

Correo electrónico