Ir al contenido principal

Cálculo de pi mediante la fórmula de Brouncker

El mes de marzo es el mes de pi, ya que el 14 de marzo (3/14) es el día de pi. Con ese motivo, el pasado 3 de marzo se publicó en Twitter un mensaje con la fórmula de Brouncker para el cálculo de pi Cálculo de pi mediante la fórmula de Brouncker

La primeras aproximaciones son

     a(1) = 4                                  =  4
     a(2) = 4/(1+1^2)                          =  2.0
     a(3) = 4/(1+1^2/(2+3^2))                  =  3.666666666666667
     a(4) = 4/(1+1^2/(2+3^2/(2+5^2)))          =  2.8
     a(5) = 4/(1+1^2/(2+3^2/(2+5^2/(2+7^2))))  =  3.395238095238095

Definir las funciones

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

tales que

  • (aproximacionPi n) es la n-ésima aproximación de pi con la fórmula de Brouncker. Por ejemplo,
     aproximacionPi 1      ==  4.0
     aproximacionPi 2      ==  2.0
     aproximacionPi 3      ==  3.666666666666667
     aproximacionPi 4      ==  2.8
     aproximacionPi 5      ==  3.395238095238095
     aproximacionPi 10     ==  3.0301437124966535
     aproximacionPi 1000   ==  3.1405916523380406
     aproximacionPi 1001   ==  3.142592653839793
     aproximacionPi 10000  ==  3.141492643588543
     aproximacionPi 10001  ==  3.1416926535900433
     pi                    ==  3.141592653589793
  • (grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi para k en xs. Por ejemplo, (grafica [10..100]) dibuja Cálculo de pi mediante la fórmula de Brouncker

Soluciones

import Graphics.Gnuplot.Simple (Attribute (Key, PNG), plotList)

-- 1ª solución
-- ===========

aproximacionPi :: Int -> Double
aproximacionPi 1 = 4
aproximacionPi n = 4/(1 + 1/(aux 1 3))
  where aux a b | a == n-1  = 1
                | otherwise = 2 + (b^2)/(aux (a+1) (b+2))

-- El cálculo es
--    aproximacionPi 2
--      = 4/(1 + 1/(aux 1 3))
--      = 4/(1 + 1/1)
--      = 2.0
--
--    aproximacionPi 3
--      = 4/(1 + 1/(aux 1 3))
--      = 4/(1 + 1/(2 + 3^2/(aux 2 5))
--      = 4/(1 + 1/(2 + 3^2/1))
--      = 3.666666666666667
--
--    aproximacionPi 4
--      = 4/(1 + 1/(aux 1 3))
--      = 4/(1 + 1/(2 + 3^2/(aux 2 5))
--      = 4/(1 + 1/(2 + 3^2/(2 + 5^2/(aux 3 7))))
--      = 4/(1 + 1/(2 + 3^2/(2 + 5^2/1)))
--      = 2.8

-- 2ª solución
-- ===========

aproximacionPi2 :: Int -> Double
aproximacionPi2 n =
  aproximacionFC n fraccionPi

-- fraccionPi es la representación de la fracción continua de pi como un
-- par de listas infinitas.
fraccionPi :: [(Integer, Integer)]
fraccionPi = zip (0 : 1 : [2,2..]) (4 : map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción
-- continua fc (como un par de listas).
aproximacionFC :: Int -> [(Integer, Integer)] -> Double
aproximacionFC n =
  foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

-- Gráfica
-- =======

grafica :: [Int] -> IO ()
grafica xs =
  plotList [ Key Nothing
           , PNG "Calculo_de_pi_mediante_la_formula_de_Brouncker_1.png"
           ]
           [(k,aproximacionPi k) | k <- xs]