Menu Close

Etiqueta: abs

Cálculo de pi mediante la serie de Nilakantha

Una serie infinita para el cálculo de pi, publicada por Nilakantha en el siglo XV, es
Calculo_de_pi_mediante_la_serie_de_Nilakantha

Definir las funciones

   aproximacionPi :: Int -> Double
   tabla          :: FilePath -> [Int] -> IO ()

tales que

  • (aproximacionPi n) es la n-ésima aproximación de pi obtenido sumando los n primeros términos de la serie de Nilakantha. Por ejemplo,
     aproximacionPi 0        ==  3.0
     aproximacionPi 1        ==  3.1666666666666665
     aproximacionPi 2        ==  3.1333333333333333
     aproximacionPi 3        ==  3.145238095238095
     aproximacionPi 4        ==  3.1396825396825396
     aproximacionPi 5        ==  3.1427128427128426
     aproximacionPi 10       ==  3.1414067184965018
     aproximacionPi 100      ==  3.1415924109719824
     aproximacionPi 1000     ==  3.141592653340544
     aproximacionPi 10000    ==  3.141592653589538
     aproximacionPi 100000   ==  3.1415926535897865
     aproximacionPi 1000000  ==  3.141592653589787
     pi                      ==  3.141592653589793
  • (tabla f ns) escribe en el fichero f las n-ésimas aproximaciones de pi, donde n toma los valores de la lista ns, junto con sus errores. Por ejemplo, al evaluar la expresión
     tabla "AproximacionesPi.txt" [0,10..100]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|   10 | 3.141406718497 | 0.000185935093 |
|   20 | 3.141565734659 | 0.000026918931 |
|   30 | 3.141584272675 | 0.000008380915 |
|   40 | 3.141589028941 | 0.000003624649 |
|   50 | 3.141590769850 | 0.000001883740 |
|   60 | 3.141591552546 | 0.000001101044 |
|   70 | 3.141591955265 | 0.000000698325 |
|   80 | 3.141592183260 | 0.000000470330 |
|   90 | 3.141592321886 | 0.000000331704 |
|  100 | 3.141592410972 | 0.000000242618 |
+------+----------------+----------------+

al evaluar la expresión

     tabla "AproximacionesPi.txt" [0,500..5000]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|  500 | 3.141592651602 | 0.000000001988 |
| 1000 | 3.141592653341 | 0.000000000249 |
| 1500 | 3.141592653516 | 0.000000000074 |
| 2000 | 3.141592653559 | 0.000000000031 |
| 2500 | 3.141592653574 | 0.000000000016 |
| 3000 | 3.141592653581 | 0.000000000009 |
| 3500 | 3.141592653584 | 0.000000000006 |
| 4000 | 3.141592653586 | 0.000000000004 |
| 4500 | 3.141592653587 | 0.000000000003 |
| 5000 | 3.141592653588 | 0.000000000002 |
+------+----------------+----------------+

Soluciones

import Text.Printf
 
-- 1ª solución
-- ===========
 
aproximacionPi :: Int -> Double
aproximacionPi n = serieNilakantha !! n
 
serieNilakantha :: [Double]
serieNilakantha = scanl1 (+) terminosNilakantha
 
terminosNilakantha :: [Double]
terminosNilakantha = zipWith (/) numeradores denominadores
  where numeradores   = 3 : cycle [4,-4]
        denominadores = 1 : [n*(n+1)*(n+2) | n <- [2,4..]]
 
-- 2ª solución
-- ===========
 
aproximacionPi2 :: Int -> Double
aproximacionPi2 = aux 3 2 1
  where aux x _ _ 0 = x
        aux x y z m =
          aux (x+4/product[y..y+2]*z) (y+2) (negate z) (m-1)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> aproximacionPi (2*10^5)
--    3.141592653589787
--    (0.82 secs, 145,964,728 bytes)
--    λ> aproximacionPi2 (2*10^5)
--    3.141592653589787
--    (2.27 secs, 432,463,496 bytes)
--    λ> aproximacionPi (3*10^5)
--    3.141592653589787
--    (0.34 secs, 73,056,488 bytes)
--    λ> aproximacionPi2 (3*10^5)
--    3.141592653589787
--    (3.24 secs, 648,603,824 bytes)
 
-- Definicioń de tabla
-- ===================
 
tabla :: FilePath -> [Int] -> IO ()
tabla f ns = do
  writeFile f (tablaAux ns)
 
tablaAux :: [Int] -> String
tablaAux ns =
     linea
  ++ cabecera
  ++ linea
  ++ concat [printf "| %4d | %.12f | %.12f |\n" n a e
            | n <- ns
            , let a = aproximacionPi n
            , let e = abs (pi - a)]
  ++ linea
 
linea :: String
linea = "+------+----------------+----------------+\n"
 
cabecera :: String
cabecera = "| n    | Aproximación   | Error          |\n"

La conjetura de Levy

Hyman Levy observó que

    7 = 3 + 2 x 2
    9 = 3 + 2 x 3 =  5 + 2 x 2
   11 = 5 + 2 x 3 =  7 + 2 x 2
   13 = 3 + 2 x 5 =  7 + 2 x 3
   15 = 3 + 2 x 5 = 11 + 2 x 2
   17 = 3 + 2 x 7 =  7 + 2 x 5 = 11 + 2 x 3 = 13 + 2 x 2
   19 = 5 + 2 x 7 = 13 + 2 x 3

y conjeturó que todos los número impares mayores o iguales que 7 se pueden escribir como la suma de un primo y el doble de un primo. El objetivo de los siguientes ejercicios es comprobar la conjetura de Levy.

Definir las siguientes funciones

   descomposicionesLevy :: Integer -> [(Integer,Integer)]
   graficaLevy          :: Integer -> IO ()

tales que

  • (descomposicionesLevy x) es la lista de pares de primos (p,q) tales que x = p + 2q. Por ejemplo,
     descomposicionesLevy  7  ==  [(3,2)]
     descomposicionesLevy  9  ==  [(3,3),(5,2)]
     descomposicionesLevy 17  ==  [(3,7),(7,5),(11,3),(13,2)]
  • (graficaLevy n) dibuja los puntos (x,y) tales que x pertenece a [7,9..7+2x(n-1)] e y es el número de descomposiciones de Levy de x. Por ejemplo, (graficaLevy 200) dibuja
    La_conjetura_de_Levy-200

Comprobar con QuickCheck la conjetura de Levy.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck
import Graphics.Gnuplot.Simple
 
descomposicionesLevy :: Integer -> [(Integer,Integer)]
descomposicionesLevy x =
  [(p,q) | p <- takeWhile (< x) (tail primes)
         , let q = (x - p) `div` 2
         , isPrime q]
 
-- graficaLevy 300
graficaLevy :: Integer -> IO ()
graficaLevy n =
  plotList [ Key Nothing
           , XRange (7,fromIntegral (7+2*(n-1)))
           , PNG ("La_conjetura_de_Levy-" ++ show n ++ ".png")
           ]
           [(x, length (descomposicionesLevy x)) | x <- [7,9..7+2*(n-1)]] 
 
-- La propiedad es
prop_Levy :: Integer -> Bool
prop_Levy x =
  not (null (descomposicionesLevy (7 + 2 * abs x)))
 
-- La comprobación es
--    λ> quickCheck prop_Levy
--    +++ OK, passed 100 tests.

La conjetura de Gilbreath

Partiendo de los 5 primeros números primos y calculando el valor absoluto de la diferencia de cada dos números consecutivos hasta quedarse con un único número se obtiene la siguiente tabla:

   2, 3, 5, 7, 11
   1, 2, 2, 4 
   1, 0, 2
   1, 2 
   1

Se observa que todas las filas, salvo la inicial, comienzan con el número 1.

Repitiendo el proceso pero empezando con los 8 primeros números primos se obtiene la siguiente tabla:

   2, 3, 5, 7, 11, 13, 17, 19 
   1, 2, 2, 4,  2,  4,  2  
   1, 0, 2, 2,  2,  2 
   1, 2, 0, 0,  0 
   1, 2, 0, 0 
   1, 2, 0 
   1, 2 
   1

Se observa que, de nuevo, todas las filas, salvo la inicial, comienza con el número 1.

La conjetura de Gilbreath afirma que si escribimos la sucesión de números primos completa y después construimos las correspondientes sucesiones formadas por el valor absoluto de la resta de cada pareja de números consecutivos, entonces todas esas filas que obtenemos comienzan siempre por 1.

El objetivo de este ejercicio es comprobar experimentalmente dicha conjetura.

Para la representación, usaremos la simétrica de la que hemos comentado anteriormente; es decir,

    2
    3, 1
    5, 2, 1
    7, 2, 0, 1
   11, 4, 2, 2, 1
   13, 2, 2, 0, 2, 1
   17, 4, 2, 0, 0, 2, 1
   19, 2, 2, 0, 0, 0, 2, 1

en la que la primera columna son los números primos y el elemento de la fila i y columna j (con i, j > 1) es el valor absoluto de la diferencia de los elementos (i,j-1) e (i-1,j-1).

Definir las siguientes funciones

   siguiente           :: Integer -> [Integer] -> [Integer]
   triangulo           :: [[Integer]]
   conjetura_Gilbreath :: Int -> Bool

tales que

  • (siguiente x ys) es la línea siguiente de la ys que empieza por x en la tabla de Gilbreath; es decir, si ys es [y1,y2,…,yn], entonces (siguiente x ys) es [x,|y1-x|,|y2-|y1-x||,…] Por ejemplo,
     siguiente  7 [5,2,1]               ==  [7,2,0,1]
     siguiente 29 [23,4,2,0,0,0,0,2,1]  ==  [29,6,2,0,0,0,0,0,2,1]
  • triangulo es el triángulo de Gilbreath. Por ejemplo,
     ghci> take 10 triangulo
     [[ 2],
      [ 3,1],
      [ 5,2,1],
      [ 7,2,0,1],
      [11,4,2,2,1],
      [13,2,2,0,2,1],
      [17,4,2,0,0,2,1],
      [19,2,2,0,0,0,2,1],
      [23,4,2,0,0,0,0,2,1],
      [29,6,2,0,0,0,0,0,2,1]]
  • (conjeturaGilbreath n) se verifica si se cumple la conjetura de Gilbreath para los n primeros números primos; es decir, en el triángulo de Gilbreath cuya primera columna son los n primeros números primos, todas las filas a partir de la segunda terminan en 1. Por ejemplo,
     ghci> conjeturaGilbreath 1000
     True

Soluciones

import Data.Numbers.Primes
 
siguiente :: Integer -> [Integer] -> [Integer]
siguiente x ys = scanl (\m n -> abs (m-n)) x ys 
 
triangulo :: [[Integer]]
triangulo = 
    [2] : [siguiente x ys | (x,ys) <- zip (tail primes) triangulo]
 
conjeturaGilbreath :: Int -> Bool
conjeturaGilbreath n = all p (tail (take n triangulo))
  where p xs = last xs == 1

Máxima distancia en árbol

Los árboles binarios con valores en las hojas y en los nodos se definen por

   data Arbol a = H a
                | N a (Arbol a) (Arbol a) 
     deriving (Eq, Show)

Por ejemplo, el árbol

         10
        /  \
       /    \
      8      1
     / \    / \
    3   9  2   6

se puede representar por

   ejArbol :: Arbol Int
   ejArbol = N 10 (N 8 (H 3) (H 9))
                  (N 1 (H 2) (H 6))

La distancia entre un padre y un hijo en el árbol es el valor absoluto de la diferencia de sus valores. Por ejemplo, la distancia de 10 a 8 es 2 y de 1 a 6 es 5.

Definir la función

   maximaDistancia :: (Num a, Ord a) => Arbol a -> a

tal que (maximaDistancia a) es la máxima distancia entre un padre y un hijo del árbol a. Por ejemplo,

   maximaDistancia ejArbol                                     ==  9
   maximaDistancia (N 1 (N 8 (H 3) (H 9)) (N 1  (H 2) (H 6)))  ==  7
   maximaDistancia (N 8 (N 8 (H 3) (H 9)) (N 10 (H 2) (H 6)))  ==  8

Soluciones

Cadenas opuestas

La opuesta de una cadena de letras es la cadena obtenida cambiando las minúsculas por mayúsculas y las minúsculas por mayúsculas. Por ejemplo, la opuesta de “SeViLLa” es “sEvIllA”.

Definir la función

   esOpuesta :: String -> String -> Bool

tal que (esOpuesta s1 s2) se verifica si las cadenas de letras s1 y s2 son opuestas. Por ejemplo,

   esOpuesta "ab" "AB"      `== True
   esOpuesta "aB" "Ab"      `== True
   esOpuesta "aBcd" "AbCD"  `== True
   esOpuesta "aBcde" "AbCD" `== False
   esOpuesta "AB" "Ab"      `== False
   esOpuesta "" ""          `== True

Soluciones

import Data.Char (isLower, isUpper, ord, toLower, toUpper)
 
-- 1ª solución (por comprensión)
-- =============================
 
esOpuesta1 :: String -> String -> Bool
esOpuesta1 s1 s2 = [opuesto c | c <- s1] == s2
 
opuesto :: Char -> Char
opuesto c
  | isLower c = toUpper c
  | otherwise = toLower c
 
-- 2ª solución (con map)
-- =====================
 
esOpuesta2 :: String -> String -> Bool
esOpuesta2 s1 s2 = map opuesto s1 == s2
 
-- 3ª solución (por recursión)
-- ===========================
 
esOpuesta3 :: String -> String -> Bool
esOpuesta3 "" ""           = True
esOpuesta3 (c1:s1) (c2:s2) = esOpuesto c1 c2 && esOpuesta3 s1 s2
esOpuesta3 _ _             = False
 
esOpuesto :: Char -> Char -> Bool
esOpuesto c1 c2 = abs (ord c1 - ord c2) == 32

Distancias entre primos consecutivos

Los 15 primeros números primos son

   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47

Las distancias entre los elementos consecutivos son

   1, 2, 2, 4, 2,  4,  2,  4,  6,  2,  6,  4,  2,  4

La distribución de las distancias es

   (1,1), (2,6), (4,5), (6,2)

(es decir, el 1 aparece una vez, el 2 aparece 6 veces, etc.) La frecuencia de las distancias es

   (1,7.142857), (2,42.857143), (4,35.714287), (6,14.285714)

(es decir, el 1 aparece el 7.142857%, el 2 el 42.857143% etc.)

Definir las funciones

  cuentaDistancias        :: Int -> [(Int,Int)]
  frecuenciasDistancias   :: Int -> [(Int,Float)]
  graficas                :: [Int] -> IO ()
  distanciasMasFrecuentes :: Int -> [Int]

tales que

  • (cuentaDistancias n) es la distribución de distancias entre los n primeros primos consecutivos. Por ejemplo,
     λ> cuentaDistancias 15
     [(1,1),(2,6),(4,5),(6,2)]
     λ> cuentaDistancias 100
     [(1,1),(2,25),(4,26),(6,25),(8,7),(10,7),(12,4),(14,3),(18,1)]
  • (frecuenciasDistancias n) es la frecuencia de distancias entre los n primeros primos consecutivos. Por ejemplo,
     λ> frecuenciasDistancias 15
     [(1,7.142857),(2,42.857143),(4,35.714287),(6,14.285714)]
     λ> frecuenciasDistancias 30
     [(1,3.4482758),(2,34.482758),(4,34.482758),(6,24.137932),(8,3.4482758)]
  • (graficas ns) dibuja las gráficas de (frecuenciasDistancias k) para k en ns. Por ejemplo, (graficas [10,20,30]) dibuja
    Distancias_entre_primos_consecutivos1
    (graficas [1000,2000,3000]) dibuja
    Distancias_entre_primos_consecutivos2
    y (graficas [100000,200000,300000]) dibuja
    Distancias_entre_primos_consecutivos3
  • (distanciasMasFrecuentes n) es la lista de las distancias más frecuentes entre los elementos consecutivos de la lista de los n primeros primos. Por ejemplo,
     distanciasMasFrecuentes 15     ==  [2]
     distanciasMasFrecuentes 26     ==  [2,4]
     distanciasMasFrecuentes 32     ==  [4]
     distanciasMasFrecuentes 41     ==  [2,4,6]
     distanciasMasFrecuentes 77     ==  [6]
     distanciasMasFrecuentes 160    ==  [4,6]
     distanciasMasFrecuentes (10^6) == [6]

Comprobar con QuickCheck si para todo n > 160 se verifica que (distanciasMasFrecuentes n) es [6].

Soluciones

import Data.Numbers.Primes
import qualified Data.Map as M 
import Graphics.Gnuplot.Simple
import Test.QuickCheck
 
cuentaDistancias :: Int -> [(Int,Int)]
cuentaDistancias = M.toList . dicDistancias
 
dicDistancias :: Int -> M.Map Int Int
dicDistancias n =
  M.fromListWith (+) (zip (distancias n) (repeat 1))
 
distancias :: Int -> [Int]
distancias n =
  zipWith (-) (tail xs) xs
  where xs = take n primes
 
frecuenciasDistancias :: Int -> [(Int,Float)]
frecuenciasDistancias n =
  [(k,(100 * fromIntegral x) / n1) | (k,x) <- cuentaDistancias n]
  where n1 = fromIntegral (n-1)
 
graficas :: [Int] -> IO ()
graficas ns =
  plotLists [Key Nothing]
            (map frecuenciasDistancias ns)
 
distanciasMasFrecuentes :: Int -> [Int]
distanciasMasFrecuentes n =
  M.keys (M.filter (==m) d)
  where d = dicDistancias n
        m = maximum (M.elems d)
 
-- La propiedad es
prop_distanciasMasFrecuentes :: Int -> Bool
prop_distanciasMasFrecuentes n =
  distanciasMasFrecuentes (161 + abs n) == [6]

Cálculo de pi mediante la serie de Nilakantha

Una serie infinita para el cálculo de pi, publicada por Nilakantha en el siglo XV, es
Calculo_de_pi_mediante_la_serie_de_Nilakantha

Definir las funciones

   aproximacionPi :: Int -> Double
   tabla          :: FilePath -> [Int] -> IO ()

tales que

  • (aproximacionPi n) es la n-ésima aproximación de pi obtenido sumando los n primeros términos de la serie de Nilakantha. Por ejemplo,
     aproximacionPi 0        ==  3.0
     aproximacionPi 1        ==  3.1666666666666665
     aproximacionPi 2        ==  3.1333333333333333
     aproximacionPi 3        ==  3.145238095238095
     aproximacionPi 4        ==  3.1396825396825396
     aproximacionPi 5        ==  3.1427128427128426
     aproximacionPi 10       ==  3.1414067184965018
     aproximacionPi 100      ==  3.1415924109719824
     aproximacionPi 1000     ==  3.141592653340544
     aproximacionPi 10000    ==  3.141592653589538
     aproximacionPi 100000   ==  3.1415926535897865
     aproximacionPi 1000000  ==  3.141592653589787
     pi                      ==  3.141592653589793
  • (tabla f ns) escribe en el fichero f las n-ésimas aproximaciones de pi, donde n toma los valores de la lista ns, junto con sus errores. Por ejemplo, al evaluar la expresión
     tabla "AproximacionesPi.txt" [0,10..100]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|   10 | 3.141406718497 | 0.000185935093 |
|   20 | 3.141565734659 | 0.000026918931 |
|   30 | 3.141584272675 | 0.000008380915 |
|   40 | 3.141589028941 | 0.000003624649 |
|   50 | 3.141590769850 | 0.000001883740 |
|   60 | 3.141591552546 | 0.000001101044 |
|   70 | 3.141591955265 | 0.000000698325 |
|   80 | 3.141592183260 | 0.000000470330 |
|   90 | 3.141592321886 | 0.000000331704 |
|  100 | 3.141592410972 | 0.000000242618 |
+------+----------------+----------------+

y al evaluar la expresión

     tabla "AproximacionesPi.txt" [0,500..5000]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|  500 | 3.141592651602 | 0.000000001988 |
| 1000 | 3.141592653341 | 0.000000000249 |
| 1500 | 3.141592653516 | 0.000000000074 |
| 2000 | 3.141592653559 | 0.000000000031 |
| 2500 | 3.141592653574 | 0.000000000016 |
| 3000 | 3.141592653581 | 0.000000000009 |
| 3500 | 3.141592653584 | 0.000000000006 |
| 4000 | 3.141592653586 | 0.000000000004 |
| 4500 | 3.141592653587 | 0.000000000003 |
| 5000 | 3.141592653588 | 0.000000000002 |
+------+----------------+----------------+

Nota: Este ejercicio ha sido propuesto por Manuel Herrera.

Referencias

Soluciones

import Text.Printf
 
-- 1ª definición
-- =============
 
aproximacionPi :: Int -> Double
aproximacionPi n = serieNilakantha !! n
 
serieNilakantha :: [Double]
serieNilakantha = scanl1 (+) terminosNilakantha
 
terminosNilakantha :: [Double]
terminosNilakantha = zipWith (/) numeradores denominadores
  where numeradores   = 3 : cycle [4,-4]
        denominadores = 1 : [n*(n+1)*(n+2) | n <- [2,4..]]
 
tabla :: FilePath -> [Int] -> IO ()
tabla f ns = do
  writeFile f (tablaAux ns)
 
tablaAux :: [Int] -> String
tablaAux ns =
     linea
  ++ cabecera
  ++ linea
  ++ concat [printf "| %4d | %.12f | %.12f |\n" n a e
            | n <- ns
            , let a = aproximacionPi n
            , let e = abs (pi - a)]
  ++ linea
 
linea :: String
linea = "+------+----------------+----------------+\n"
 
cabecera :: String
cabecera = "| n    | Aproximación   | Error          |\n"

Cálculo de pi usando la fórmula de Vieta

La fórmula de Vieta para el cálculo de pi es la siguiente
Calculo_de_pi_usando_la_formula_de_Vieta

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

tales que

  • (aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,
     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795
  • (errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,
     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))
 
-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)
 
-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]
 
-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]
 
-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

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)

Árboles continuos

Los árboles binarios se pueden representar con el de tipo de dato algebraico

   data Arbol a = H
                | N a (Arbol a) (Arbol a)
     deriving Show

Por ejemplo, los árboles

       3                7     
      / \              / \    
     2   4            5   8   
    / \   \          / \   \  
   1   3   5        6   4   10

se representan por

   ej1, ej2 :: Arbol Int
   ej1 = N 3 (N 2 (N 1 H H) (N 3 H H)) (N 4 H (N 5 H H))
   ej2 = N 7 (N 5 (N 6 H H) (N 4 H H)) (N 8 H (N 10 H H))

Un árbol binario es continuo si el valor absoluto de la diferencia de los elementos adyacentes es 1. Por ejemplo, el árbol ej1 es continuo ya que el valor absoluto de sus pares de elementos adyacentes son

   |3-2| = |2-1| = |2-3| = |3-4| = |4-5| = 1

En cambio, el ej2 no lo es ya que |8-10| ≠ 1.

Definir la función

   esContinuo :: (Num a, Eq a) => Arbol a -> Bool

tal que (esContinuo x) se verifica si el árbol x es continuo. Por ejemplo,

   esContinuo ej1  ==  True
   esContinuo ej2  ==  False

Soluciones

data Arbol a = H
             | N a (Arbol a) (Arbol a)
  deriving Show
 
ej1, ej2 :: Arbol Int
ej1 = N 3 (N 2 (N 1 H H) (N 3 H H)) (N 4 H (N 5 H H))
ej2 = N 7 (N 5 (N 6 H H) (N 4 H H)) (N 8 H (N 10 H H))
 
-- 1ª solución
-- ===========
 
esContinuo :: (Num a, Eq a) => Arbol a -> Bool
esContinuo H         = True
esContinuo (N _ H H) = True
esContinuo (N x i@(N y _ _) H) =
  abs (x - y) == 1 && esContinuo i
esContinuo (N x H d@(N y _ _)) =
  abs (x - y) == 1 && esContinuo d
esContinuo (N x i@(N y _ _) d@(N z _ _)) =
  abs (x - y) == 1 && esContinuo i && abs (x - z) == 1 && esContinuo d 
 
-- 2ª solución
-- ===========
 
esContinuo2 :: (Num a, Eq a) => Arbol a -> Bool
esContinuo2 x =
  all esContinua (ramas x)
 
-- (ramas x) es la lista de las ramas del árbol x. Por ejemplo,
--    ramas ej1  ==  [[3,2,1],[3,2,3],[3,4,5]]
--    ramas ej2  ==  [[7,5,6],[7,5,4],[7,8,10]]
ramas :: Arbol a -> [[a]]
ramas H         = []
ramas (N x H H) = [[x]]
ramas (N x i d) = [x : xs | xs <- ramas i ++ ramas d]
 
-- (esContinua xs) se verifica si el valor absoluto de la diferencia de
-- los elementos adyacentes de xs es 1. Por ejemplo, 
--    esContinua [3,2,3]   ==  True
--    esContinua [7,8,10]  ==  False
esContinua :: (Num a, Eq a) => [a] -> Bool
esContinua xs =
  and [abs (x - y) == 1 | (x, y) <- zip xs (tail xs)]

Referencias