Menu Close

Cálculo de pi mediante la variante de Euler de la serie armónica

José A. Alonso,

En el artículo El desarrollo más bello de Pi como suma infinita, Miguel Ángel Morales comenta el desarrollo de pi publicado por Leonhard Euler en su libro «Introductio in Analysis Infinitorum» (1748).

El desarrollo es el siguiente
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_1
y se obtiene a partir de la serie armónica
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_2
modificando sólo el signo de algunos términos según el siguiente criterio:

  • Dejamos un + cuando el denominador de la fracción sea un 2 o un primo de la forma 4m-1.
  • Cambiamos a – si el denominador de la fracción es un primo de la forma 4m+1.
  • Si el número es compuesto ponemos el signo que quede al multiplicar los signos correspondientes a cada factor.

Por ejemplo,

  • la de denominador 3 = 4×1-1 lleva un +,
  • la de denominador 5 = 4×1+1 lleva un -,
  • la de denominador 13 = 4×3+1 lleva un -,
  • la de denominador 6 = 2×3 lleva un + (porque los dos llevan un +),
  • la de denominador 10 = 2×5 lleva un – (porque el 2 lleva un + y el 5 lleva un -) y
  • la de denominador 50 = 5x5x2 lleva un + (un – por el primer 5, otro – por el segundo 5 y un + por el 2).

Definir las funciones

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

tales que

  • (aproximacionPi n) es la aproximación de pi obtenida sumando los n primeros términos de la serie de Euler. Por ejemplo.
     aproximacionPi 1        ==  1.0
     aproximacionPi 10       ==  2.3289682539682537
     aproximacionPi 100      ==  2.934318000847734
     aproximacionPi 1000     ==  3.0603246224585128
     aproximacionPi 10000    ==  3.1105295744825403
     aproximacionPi 100000   ==  3.134308801935256
     aproximacionPi 1000000  ==  3.1395057903490806
  • (grafica n) dibuja la gráfica de las aproximaciones de pi usando k sumando donde k toma los valores de la lista [100,110..n]. Por ejemplo, al evaluar (grafica 4000) se obtiene
    Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_3.png

Soluciones

import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
 
-- 1ª definición
-- =============
 
aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [1 / fromIntegral (k * signo k) | k <- [1..n]] 
 
signoPrimo :: Int -> Int
signoPrimo 2 = 1
signoPrimo p | p `mod` 4 == 3 = 1
             | otherwise      = -1
 
signo :: Int -> Int
signo n | isPrime n = signoPrimo n
        | otherwise = product (map signoPrimo (primeFactors n))
 
-- 2ª definición
-- =============
 
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = serieEuler !! (n-1)
 
serieEuler :: [Double]
serieEuler =
  scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]
 
-- Definición de grafica
-- =====================
 
grafica :: Int -> IO ()
grafica n = 
    plotList [Key Nothing]
             [(k,aproximacionPi2 k) | k <- [100,110..n]]

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>
Posted in Inicial
Tagged , , , , , , , , , ,

3 Comments

  1. rebgongor
    Responder 
    import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
 
aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [fromIntegral (signo k) / fromIntegral k | k <- [1..n]]
 
signo :: Int -> Int
signo n
  | (not . isPrime) n         = product $ map signo $ primeFactors n
  | primoNegativo n           = -1
  | n == 2 || primoPositivo n = 1
  where primoPositivo n = elem n [4*k-1 | k <- [1..fromIntegral (n-1)]]
        primoNegativo n = elem n [4*k+1 | k <- [1..fromIntegral (n-1)]]
 
grafica :: Int -> IO ()
grafica n =
  plotLists [Key Nothing, PNG "graficaEuler.png"]
            [map aproximacionPi [100,110..n]]
  2. Enrique Zubiría
    Responder 
    import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
 
factores :: Int -> [Int]
factores n
  | isPrime n = [n]
  | otherwise = primerFactor : factores (n `div` primerFactor)
  where primerFactor = head [m | m <- [2..n], 0 == mod n m]
 
signo :: Int -> Double
signo n
  | n == 1 = 1
  | n == 2 = 1
  | isPrime n && (mod (n+1) 4) == 0 = 1
  | isPrime n && (mod (n-1) 4) == 0 = (-1)
  | otherwise = product (map signo (factores n))
 
aproximacionPi :: Int -> Double
aproximacionPi n = sum [(signo m)/(fromIntegral m) | m <- [1..n]]
 
grafica        :: Int -> IO ()
grafica n = do
  plotList [ Key Nothing
           , PNG "calculoPiEulerSerieArmonica.png"
           ]
           (map aproximacionPi [100,110..n])
  3. javjimord
    Responder 
    import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple
 
aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [fromIntegral m * (1/fromIntegral x) | (m,x) <- zip ys xs]
  where ys = [product (map signo $ primeFactors x) | x <- [1..n]]
        xs = [1..n]
 
signo :: Int -> Int
signo p
  | p `mod` 4 == 1 = (-1)
  | otherwise      = 1
 
grafica :: Int -> IO ()
grafica n =
  plotList [Key Nothing]
           [(x,aproximacionPi x) | x <- [100,110..n]]

Escribe tu solución

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.