Menu Close

Cálculo de dígitos de pi y su distribución

Se pueden generar los dígitos de Pi, como se explica en el artículo Unbounded spigot algorithms for the digits of pi c0on la función digitosPi definida por

   digitosPi :: [Integer]
   digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) = 
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

Por ejemplo,

   λ> take 25 digitosPi
   [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

La distribución de los primeros 25 dígitos de pi es [0,2,3,5,3,3,3,1,2,3] ya que el 0 no aparece, el 1 ocurre 2 veces, el 3 ocurre 3 veces, el 4 ocurre 5 veces, …

Usando digitosPi, definir las siguientes funciones

   distribucionDigitosPi :: Int -> [Int]
   frecuenciaDigitosPi   :: Int -> [Double]

tales que

  • (distribucionDigitosPi n) es la distribución de los n primeros dígitos de pi. Por ejemplo,
     λ> distribucionDigitosPi 10
     [0,2,1,2,1,2,1,0,0,1]
     λ> distribucionDigitosPi 100
     [8,8,12,12,10,8,9,8,12,13]
     λ> distribucionDigitosPi 1000
     [93,116,103,103,93,97,94,95,101,105]
     λ> distribucionDigitosPi 5000
     [466,531,496,460,508,525,513,488,492,521]
  • (frecuenciaDigitosPi n) es la frecuencia de los n primeros dígitos de pi. Por ejemplo,
   λ> frecuenciaDigitosPi 10
   [0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]
   λ> frecuenciaDigitosPi 100
   [8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]
   λ> frecuenciaDigitosPi 1000
   [9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]
   λ> frecuenciaDigitosPi 5000
   [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

Soluciones

import Data.Array
import Data.List (group, sort)
 
digitosPi :: [Integer]
digitosPi = g(1,0,1,1,3,3) where
  g (q,r,t,k,n,l) = 
    if 4*q+r-t < n*t
    then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
    else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
 
-- 1ª definición
-- =============
 
distribucionDigitosPi :: Int -> [Int]
distribucionDigitosPi n =
  elems (accumArray (+) 0 (0,9) [ (i,1)
                                | i <- take n digitosPi]) 
 
frecuenciaDigitosPi :: Int -> [Double]
frecuenciaDigitosPi n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]
  where m = fromIntegral n
 
-- 2ª definición
-- =============
 
distribucionDigitosPi2 :: Int -> [Int]
distribucionDigitosPi2 n =
  [length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]
 
frecuenciaDigitosPi2 :: Int -> [Double]
frecuenciaDigitosPi2 n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]
  where m = fromIntegral n
 
-- Comparación de eficiencia
-- =========================
 
--    λ> last (take 5000 digitosPi)
--    2
--    (4.47 secs, 3,927,848,448 bytes)
--    λ> frecuenciaDigitosPi 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.01 secs, 0 bytes)
--    λ> frecuenciaDigitosPi2 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.02 secs, 0 bytes)

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

4 soluciones de “Cálculo de dígitos de pi y su distribución

  1. Enrique Zubiría
    digitosPi :: [Integer]
    digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) =
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
     
    distribucionDigitosPi :: Int -> [Int]
    distribucionDigitosPi n = map f [0..9]
      where digPi = take n digitosPi
            f x   = length $ filter (x==) digPi
     
    frecuenciaDigitosPi :: Int -> [Double]
    frecuenciaDigitosPi n =
      map (x -> 100*x/(fromIntegral n))
          (map fromIntegral (distribucionDigitosPi n))
  2. juabaerui
    digitosPi :: [Integer]
    digitosPi = g(1,0,1,1,3,3) where
      g (q,r,t,k,n,l) = 
        if 4*q+r-t < n*t
        then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
        else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
     
    distribucionDigitosPi :: Int -> [Int]
    distribucionDigitosPi n = [a 0, a 1, a 2, a 3, a 4, a 5, a 6, a 7, a 8, a 9]
      where a x = apariciones x (take n digitosPi)
            apariciones k []     = 0
            apariciones k (x:xs) | k == x    = 1 + apariciones k xs
                                 | otherwise = apariciones k xs
     
    frecuenciaDigitosPi :: Int -> [Double]
    frecuenciaDigitosPi n = map (aproximandecimales 2) (map (*(100/k)) ys)
      where k  = fromIntegral n :: Double
            ys = map fromIntegral (distribucionDigitosPi n)
     
    aproximandecimales :: Int -> Double -> Double
    aproximandecimales n x = (fromIntegral (round (x*10^n))) / 10^n
  3. rebgongor
    digitosPi :: [Integer]
    digitosPi = g(1,0,1,1,3,3) where
         g (q,r,t,k,n,l) = 
           if 4*q+r-t < n*t
           then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
           else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
     
    distribucionDigitosPi :: Int -> [Int]
    distribucionDigitosPi n =
      map (x -> length $ filter (==x) $ take n digitosPi) [0..9]
     
    frecuenciaDigitosPi :: Int -> [Double]
    frecuenciaDigitosPi n =
      map (*100) [fromIntegral x/fromIntegral n | x <- distribucionDigitosPi n]
  4. javjimord
    digitosPi :: [Integer]
    digitosPi = g(1,0,1,1,3,3) where
         g (q,r,t,k,n,l) = 
           if 4*q+r-t < n*t
           then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
           else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
     
    distribucionDigitosPi :: Int -> [Int]
    distribucionDigitosPi n =
      [ocurrencias x (take n digitosPi) | x <- [0..9]]
     
    ocurrencias :: Integer -> [Integer] -> Int
    ocurrencias x xs
      | x `notElem` xs = 0
      | (head xs) == x = 1 + ocurrencias x (tail xs)
      | otherwise      = ocurrencias x (tail xs)
     
    frecuenciaDigitosPi :: Int -> [Double]
    frecuenciaDigitosPi n =
      [100*x / fromIntegral n | x <- map fromIntegral (distribucionDigitosPi n)]

Leave a Reply

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