Menu Close

Cálculo de pi usando el producto de Wallis

El producto de Wallis es una expresión, descubierta por John Wallis en 1655, para representar el valor de π y que establece que:

    π     2     2     4     4     6     6     8     8
   --- = --- · --- · --- · --- · --- · --- · --- · --- ···
    2     1     3     3     5     5     7     7     9

Definir las funciones

   factoresWallis  :: [Rational]
   productosWallis :: [Rational]
   aproximacionPi  :: Int -> Double
   errorPi         :: Double -> Int

tales que

  • factoresWallis es la sucesión de los factores del productos de Wallis. Por ejemplo,
     λ> take 10 factoresWallis
     [2 % 1,2 % 3,4 % 3,4 % 5,6 % 5,6 % 7,8 % 7,8 % 9,10 % 9,10 % 11]
  • productosWallis es la sucesión de los productos de los primeros factores de Wallis. Por ejemplo,
     λ> take 7 productosWallis
     [2 % 1,4 % 3,16 % 9,64 % 45,128 % 75,256 % 175,2048 % 1225]
  • (aproximacionPi n) es la aproximación de pi obtenida multiplicando los n primeros factores de Wallis. Por ejemplo,
     aproximacionPi 20     ==  3.2137849402931895
     aproximacionPi 200    ==  3.1493784731686008
     aproximacionPi 2000   ==  3.142377365093878
     aproximacionPi 20000  ==  3.141671186534396
  • (errorPi x) es el menor número de factores de Wallis necesarios para obtener pi con un error menor que x. Por ejemplo,
     errorPi 0.1     ==  14
     errorPi 0.01    ==  155
     errorPi 0.001   ==  1569
     errorPi 0.0001  ==  15707

Soluciones

import Data.Ratio
 
factoresWallis :: [Rational]
factoresWallis =
  concat [[y%(y-1),  y%(y+1)] | x <- [1..], let y = 2*x]
 
productosWallis :: [Rational]
productosWallis = scanl1 (*) factoresWallis
 
aproximacionPi :: Int -> Double
aproximacionPi n =
  fromRational (2 * productosWallis !! n)
 
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]
 
-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x =
  length (takeWhile (>=x) [abs (pi - 2 * fromRational y)
                          | y <- productosWallis])
 
-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n =
  2 * productosWallis2 !! n
 
productosWallis2 :: [Double]
productosWallis2 = scanl1 (*) factoresWallis2
 
factoresWallis2 :: [Double]
factoresWallis2 =
  concat [[y/(y-1),  y/(y+1)] | x <- [1..], let y = 2*x]
 
-- 3ª definición de errorPi
errorPi3 :: Double -> Int
errorPi3 x = head [n | n <- [1..]
                     , abs (pi - aproximacionPi2 n) < x]
 
-- Comparación de eficiencia
--    λ> errorPi 0.001
--    1569
--    (0.82 secs, 374,495,816 bytes)
--
--    λ> errorPi2 0.001
--    1569
--    (0.79 secs, 369,282,320 bytes)
--
--    λ> errorPi3 0.001
--    1569
--    (0.04 secs, 0 bytes)
Medio

5 soluciones de “Cálculo de pi usando el producto de Wallis

  1. enrnarbej
    factoresWallis  :: [Rational]
    factoresWallis = zipWith (/) parR imparR
      where
        parR   = concat [[x,x] | x <- [2,4..]]
        imparR = 1 : concat [[x,x] | x <- [3,5..]]
     
    productosWallis :: [Rational]
    productosWallis = scanl (*) (2/1) (tail factoresWallis)
     
    aproximacionPi :: Int -> Double
    aproximacionPi n = 2 * fromRational (productosWallis !! n)
     
    errorPi  :: Double -> Int
    errorPi e =
      length $ takeWhile (k -> abs (pi - 2*fromRational k) >= e)
                         productosWallis
  2. albcercid
    factoresWallis  :: [Rational]
    factoresWallis = [a % b | (a,b) <- zip (r [2,4..]) (1:r [3,5..])]
      where r (x:xs) = x:x:r xs
     
    productosWallis :: [Rational]
    productosWallis = aux (2 % 1) (tail factoresWallis)
      where aux x (y:ys) = x:aux (x*y) ys
     
    aproximacionPi  :: Int -> Double
    aproximacionPi n = 2*(fromRational $ productosWallis!!n)
     
    aproximacionPi2  :: Int -> Double
    aproximacionPi2 n = 2*aux (-1) 1 factoresWallis n
      where aux n b (x:xs) y | y == n = b
                             | otherwise = aux (n+1) (b*fromRational x) xs y
     
    errorPi :: Double -> Int
    errorPi n = aux n factWallis2 (-1) 2
      where aux a (x:xs) b c | abs (pi-c) < a = b
                             | otherwise = aux a xs (b+1) (c*x)
            factWallis2 = [ a/b | (a,b) <- zip (r [2,4..]) (1:r [3,5..])]
            r (x:xs) = x:x:r xs
  3. Juanjo Ortega (juaorture)
    import Data.Ratio
     
    factoresWallis :: [Rational]
    factoresWallis = [ a % b | (a,b) <- xs ]
      where xs = zip (concatMap (replicate 2) [2,4..])
                     (1 : concatMap (replicate 2) [3,5..])
     
    productosWallis :: [Rational]
    productosWallis = aux factoresWallis
      where aux (x:y:xs) = x : aux (x*y:xs)
     
    aproximacionPi :: Int -> Double
    aproximacionPi n =  2 * fromRational x
      where x = product $ take (n+1) $ factoresWallis
     
    aproximacionPi' :: Int -> Double
    aproximacionPi' n = 2 * fromRational (productosWallis!!n)
     
    errorPi :: Double -> Int
    errorPi x =
      head [ n | n <- [1..]
               , abs (pi - aproximacionPi' n) <= x]
  4. eliguivil
    factoresWallis :: [Rational]
    factoresWallis = [n % d | (n,d) <- zip xs ys]
      where
        xs = concat [[n,n] | n <- [2,4..]]
        ys = 1: concat [[n,n] | n <- [3,5..]]
     
    productosWallis :: [Rational]
    productosWallis = scanl1 (*) factoresWallis
     
    aproximacionPi :: Int -> Double
    aproximacionPi n = 2*fromRational (productosWallis !! n)
     
    errorPi :: Double -> Int
    errorPi e = aux e 2
      where
        aux e n | abs (pi- (aproximacionPi n)) < e = n
                | otherwise = aux e (n+1)
  5. cescarde
    factoresWallis :: [Rational]
    factoresWallis = intercala [2*x/(2*x-1) | x <- [1..]] [2*x/(2*x+1) | x <- [1..]]
     
    intercala :: [a] -> [a] -> [a]
    intercala [] [] = []
    intercala [] ys = ys
    intercala xs [] = xs
    intercala (x:xs) (y:ys) = x : y : (intercala xs ys)
     
    productosWallis :: [Rational]
    productosWallis = [product (take n factoresWallis) | n <- [1..]]
     
    aproximacionPi  :: Int -> Double
    aproximacionPi n = fromRational $ product $ take (n+1) (double factoresWallis)
                   where double (x:xs) = 2*x : xs
     
    errorPi :: Double -> Int
    errorPi x = head [n | n <- [1..], abs (pi- aproximacionPi n) < x]

Escribe tu solución

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