Menu Close

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

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>

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

  1. rebgongor
    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
    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
    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]]

Leave a Reply

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