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

y se obtiene a partir de la serie armónica

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

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

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

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

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 Comentarios

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

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

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

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

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

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

Responder

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

Soluciones

Otras soluciones

3 Comentarios

Escribe tu soluciónCancel reply

Entradas recientes

Buscar

Correo electrónico

Soluciones

Otras soluciones

Se puede imprimir o compartir con

3 Comentarios

Escribe tu soluciónCancel reply

Entradas recientes

Buscar

Correo electrónico