Menu Close

Aproximación del número pi

Una forma de aproximar el número π es usando la siguiente igualdad:

    π         1     1·2     1·2·3     1·2·3·4     
   --- = 1 + --- + ----- + ------- + --------- + ....
    2         3     3·5     3·5·7     3·5·7·9

Es decir, la serie cuyo término general n-ésimo es el cociente entre el producto de los primeros n números y los primeros n números impares:

               Π i   
   s(n) =  -----------
            Π (2*i+1)

Definir la función

   aproximaPi :: Integer -> Double

tal que (aproximaPi n) es la aproximación del número π calculada con la serie anterior hasta el término n-ésimo. Por ejemplo,

   aproximaPi 10     == 3.1411060206
   aproximaPi 20     == 3.1415922987403397
   aproximaPi 30     == 3.1415926533011596
   aproximaPi 40     == 3.1415926535895466
   aproximaPi 50     == 3.141592653589793
   aproximaPi (10^4) == 3.141592653589793
   pi                == 3.141592653589793

Soluciones

import Data.Ratio ((%))
import Data.List (genericTake)
import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)
 
-- 1ª solución
-- ===========
 
aproximaPi1 :: Integer -> Double
aproximaPi1 n = 
  fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])
 
-- 2ª solución
-- ===========
 
aproximaPi2 :: Integer -> Double
aproximaPi2 0 = 2
aproximaPi2 n = 
  aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])
 
-- 3ª solución
-- ===========
 
aproximaPi3 :: Integer -> Double
aproximaPi3 n = 
  fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))
 
-- numeradores es la sucesión de los numeradores. Por ejemplo,
--    λ> take 10 numeradores
--    [1,2,6,24,120,720,5040,40320,362880,3628800]
numeradores :: [Integer]
numeradores = scanl (*) 1 [2..]
 
-- denominadores es la sucesión de los denominadores. Por ejemplo,
--    λ> take 10 denominadores
--    [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]
denominadores :: [Integer]
denominadores = scanl (*) 3 [5, 7..]
 
-- 4ª solución
-- ===========
 
aproximaPi4 :: Integer -> Double
aproximaPi4 n = 
  read (x : "." ++ xs)
  where (x:xs) = show (aproximaPi4' n)
 
aproximaPi4' :: Integer -> Integer
aproximaPi4' n = 
  2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores))) 
  where p = 10^n
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_aproximaPi :: Property
prop_aproximaPi =
  forAll (arbitrary `suchThat` (> 3)) $ \n ->
  all (=~ aproximaPi1 n)
      [aproximaPi2 n,
       aproximaPi3 n,
       aproximaPi4 n]
 
(=~) :: Double -> Double -> Bool
x =~ y = abs (x - y) < 0.001
 
-- La comprobación es
--    λ> quickCheck prop_aproximaPi
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> aproximaPi1 3000 
--    3.141592653589793
--    (4.96 secs, 27,681,824,408 bytes)
--    λ> aproximaPi2 3000 
--    3.1415926535897922
--    (3.00 secs, 20,496,194,496 bytes)
--    λ> aproximaPi3 3000 
--    3.141592653589793
--    (3.13 secs, 13,439,528,432 bytes)
--    λ> aproximaPi4 3000 
--    3.141592653589793
--    (0.09 secs, 23,142,144 bytes)

El código se encuentra en GitHub.

Posted in Ejercicio

Escribe tu solución

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