Menu Close

Etiqueta: scanl1

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

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)
 
-- 3ª solución
-- ===========
 
aproximacionPi3 :: Int -> Double
aproximacionPi3 x =
  3 + sum [(((-1)**(n+1))*4)/(2*n*(2*n+1)*(2*n+2))
          | n <- [1..fromIntegral x]]
 
 
-- Comparación de eficiencia
-- =========================
 
--    λ> aproximacionPi (10^6)
--    3.141592653589787
--    (1.35 secs, 729,373,160 bytes)
--    λ> aproximacionPi2 (10^6)
--    3.141592653589787
--    (2.96 secs, 2,161,766,096 bytes)
--    λ> aproximacionPi3 (10^6)
--    3.1415926535897913
--    (2.02 secs, 1,121,372,536 bytes)
 
-- Definicioń de tabla
-- ===================
 
tabla :: FilePath -> [Int] -> IO ()
tabla f ns = 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"

Pensamiento

Bueno es saber que los vasos
nos sirven para beber;
lo malo es que no sabemos
para que sirve la sed.

Antonio Machado

Mayor prefijo con suma acotada

Definir la función

   mayorPrefijoAcotado :: [Int] -> Int -> [Int]

tal que (mayorPrefijoAcotado xs y) es el mayor prefijo de la lista de los números enteros positivos xs cuya suma es menor o igual que y. Por ejemplo,

   mayorPrefijoAcotado [45,30,55,20,80,20] 75   ==  [45,30]
   mayorPrefijoAcotado [45,30,55,20,80,20] 140  ==  [45,30,55]
   mayorPrefijoAcotado [45,30,55,20,80,20] 180  ==  [45,30,55,20]
   length (mayorPrefijoAcotado (repeat 1) (8*10^6)) == 8000000

Soluciones

import Data.List (inits)
 
-- 1ª solución
-- ===========
 
mayorPrefijoAcotado :: [Int] -> Int -> [Int]
mayorPrefijoAcotado [] _     = []
mayorPrefijoAcotado (x:xs) y
  | x > y     = []
  | otherwise = x : mayorPrefijoAcotado xs (y-x)
 
-- 2ª solución
-- ===========
 
mayorPrefijoAcotado2 :: [Int] -> Int -> [Int]
mayorPrefijoAcotado2 xs y =
  take (longitudMayorPrefijoAcotado2 xs y) xs
 
longitudMayorPrefijoAcotado2 :: [Int] -> Int -> Int
longitudMayorPrefijoAcotado2 xs y =
  length (takeWhile (<=y) (map sum (inits xs))) - 1
 
-- 3ª solución
-- ===========
 
mayorPrefijoAcotado3 :: [Int] -> Int -> [Int]
mayorPrefijoAcotado3 xs y =
  take (longitudMayorPrefijoAcotado3 xs y) xs
 
longitudMayorPrefijoAcotado3 :: [Int] -> Int -> Int
longitudMayorPrefijoAcotado3 xs y =
  length (takeWhile (<= y) (scanl1 (+) xs))
 
-- Equivalencia
-- ============
 
-- La propiedad es
prop_equiv :: [Int] -> Int -> Bool
prop_equiv xs y =
  mayorPrefijoAcotado xs' y' == mayorPrefijoAcotado2 xs' y' &&
  mayorPrefijoAcotado xs' y' == mayorPrefijoAcotado3 xs' y' 
  where xs' = map abs xs
        y'  = abs y
 
-- La comprobación es
--    λ> quickCheck prop_equiv
--    +++ OK, passed 100 tests.
--    (0.01 secs, 2,463,688 bytes)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> length (mayorPrefijoAcotado (repeat 1) (2*10^4))
--    20000
--    (0.04 secs, 5,086,544 bytes)
--    λ> length (mayorPrefijoAcotado2 (repeat 1) (2*10^4))
--    20000
--    (11.22 secs, 27,083,980,168 bytes)
--    λ> length (mayorPrefijoAcotado3 (repeat 1) (2*10^4))
--    20000
--    (0.02 secs, 4,768,992 bytes)
--    
--    λ> length (mayorPrefijoAcotado (repeat 1) (8*10^6))
--    8000000
--    (3.19 secs, 1,984,129,832 bytes)
--    λ> length (mayorPrefijoAcotado3 (repeat 1) (8*10^6))
--    8000000
--    (1.02 secs, 1,856,130,936 bytes)

Pensamiento

Sed hombres de mal gusto. Yo os aconsejo el mal gusto para combatir los excesos de la moda.

Antonio Machado

Números construidos con los dígitos de un conjunto dado

Definir las siguientes funciones

   numerosCon      :: [Integer] -> [Integer]
   numeroDeDigitos :: [Integer] -> Integer -> Int

tales que

  • (numerosCon ds) es la lista de los números que se pueden construir con los dígitos de ds (cuyos elementos son distintos elementos del 1 al 9) . Por ejemplo,
     λ> take 22 (numerosCon [1,4,6,9])
     [1,4,6,9,11,14,16,19,41,44,46,49,61,64,66,69,91,94,96,99,111,114]
     λ> take 15 (numerosCon [4,6,9])
     [4,6,9,44,46,49,64,66,69,94,96,99,444,446,449]
     λ> take 15 (numerosCon [6,9])
     [6,9,66,69,96,99,666,669,696,699,966,969,996,999,6666]
  • (numeroDeDigitos ds k) es el número de dígitos que tiene el k-ésimo elemento (empezando a contar en 0) de la sucesión (numerosCon ds). Por ejemplo,
     numeroDeDigitos [1,4,6,9]   3  ==  1
     numeroDeDigitos [1,4,6,9]   6  ==  2
     numeroDeDigitos [1,4,6,9]  22  ==  3
     numeroDeDigitos [4,6,9]    15  ==  3
     numeroDeDigitos [6,9]      15  ==  4
     numeroDeDigitos [1,4,6,9] (10^(10^5))  ==  166097
     numeroDeDigitos   [4,6,9] (10^(10^5))  ==  209590
     numeroDeDigitos     [6,9] (10^(10^5))  ==  332192

Soluciones

import Data.List (genericIndex, genericLength)
 
-- Definición de numerosCon
-- ========================
 
numerosCon :: [Integer] -> [Integer]
numerosCon xs = [n | n <- [0..]
                   , n `tieneSusDigitosEn` xs]
 
-- (tieneSusDigitosEn xs n) se verifica si los dígitos de n están contenidos en
-- xs. Por ejemplo,
--    4149 `tieneSusDigitosEn` [1,4,6,9]  ==  True
--    4143 `tieneSusDigitosEn` [1,4,6,9]  ==  False
tieneSusDigitosEn :: Integer -> [Integer] -> Bool
tieneSusDigitosEn n xs =
  digitos n `esSubconjunto` xs
 
-- (digitos n) es la lista de los dígitos de n. Por ejemplo,
--    digitos 325  ==  [3,2,5]
digitos :: Integer -> [Integer]
digitos n = [read [c] | c <- show n]
 
-- (esSubconjunto xs ys) se verifica si xs es subconjunto de ys. Por
-- ejemplo, 
--    esSubconjunto [3,2,5] [4,2,5,7,3]  ==  True
--    esSubconjunto [3,2,5] [4,2,5,7]  ==  False
esSubconjunto :: Eq a => [a] -> [a] -> Bool
esSubconjunto xs ys = all (`elem` ys) xs
 
-- 1ª definición de numeroDeDigitos
-- ================================
 
numeroDeDigitos :: [Integer] -> Integer -> Int
numeroDeDigitos xs n =
  length (show (numerosCon xs `genericIndex` n))
 
-- 2ª definición de numeroDeDigitos
-- ================================
 
-- Observando que si el conjunto de dígitos tiene n elementos, entonces
-- hay n^k números con k dígitos.
 
numeroDeDigitos2 :: [Integer] -> Integer -> Int
numeroDeDigitos2 xs n =
  1 + length (takeWhile (<= n) (sumasDePotencias (genericLength xs)))
 
-- (potencia n) son las potencias de n. Por ejemplo,
--    take 5 (potencias 3)  ==  [3,9,27,81,243]
potencias :: Integer -> [Integer]
potencias k = iterate (*k) k 
 
-- (sumasDePotencias n) es la lista de las sumas acumuladas de las
-- potencias de n. Por ejemplo,
--    take 5 (sumasDePotencias 3)  ==  [3,12,39,120,363]
sumasDePotencias :: Integer -> [Integer]
sumasDePotencias k = scanl1 (+) (potencias k)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> numeroDeDigitos [1,4,6,9] 2000
--    6
--    (3.98 secs, 1,424,390,064 bytes)
--    λ> numeroDeDigitos2 [1,4,6,9] 2000
--    6
--    (0.01 secs, 127,464 bytes)

Operaciones con series de potencias

Una serie de potencias es una serie de la forma

   a₀ + a₁x + a₂x² + a₃x³ + ...

Las series de potencias se pueden representar mediante listas infinitas. Por ejemplo, la serie de la función exponencial es

   e^x = 1 + x + x²/2! + x³/3! + ...

y se puede representar por [1, 1, 1/2, 1/6, 1/24, 1/120, …]

Las operaciones con series se pueden ver como una generalización de las de los polinomios.

En lo que sigue, usaremos el tipo (Serie a) para representar las series de potencias con coeficientes en a y su definición es

   type Serie a = [a]

Definir las siguientes funciones

   opuesta      :: Num a => Serie a -> Serie a
   suma         :: Num a => Serie a -> Serie a -> Serie a
   resta        :: Num a => Serie a -> Serie a -> Serie a
   producto     :: Num a => Serie a -> Serie a -> Serie a
   cociente     :: Fractional a => Serie a -> Serie a -> Serie a
   derivada     :: (Num a, Enum a) => Serie a -> Serie a
   integral     :: (Fractional a, Enum a) => Serie a -> Serie a
   expx         :: Serie Rational

tales que

  • (opuesta xs) es la opuesta de la serie xs. Por ejemplo,
     λ> take 7 (opuesta [-6,-4..])
     [6,4,2,0,-2,-4,-6]
  • (suma xs ys) es la suma de las series xs e ys. Por ejemplo,
     λ> take 7 (suma [1,3..] [2,4..])
     [3,7,11,15,19,23,27]
  • (resta xs ys) es la resta de las series xs es ys. Por ejemplo,
     λ> take 7 (resta [3,5..] [2,4..])
     [1,1,1,1,1,1,1]
     λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])
     [2,4,6,8,10,12,14]
  • (producto xs ys) es el producto de las series xs e ys. Por ejemplo,
     λ> take 7 (producto [3,5..] [2,4..])
     [6,22,52,100,170,266,392]
  • (cociente xs ys) es el cociente de las series xs e ys. Por ejemplo,
     λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..])
     [2.0,4.0,6.0,8.0,10.0,12.0,14.0]
  • (derivada xs) es la derivada de la serie xs. Por ejemplo,
     λ> take 7 (derivada [2,4..])
     [4,12,24,40,60,84,112]
  • (integral xs) es la integral de la serie xs. Por ejemplo,
     λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0))
     [0.0,4.0,6.0,8.0,10.0,12.0,14.0]
  • expx es la serie de la función exponencial. Por ejemplo,
     λ> take 8 expx
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (derivada expx)
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (integral expx)
     [0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

Soluciones

type Serie a = [a] 
 
opuesta :: Num a => Serie a -> Serie a
opuesta = map negate
 
suma :: Num a => Serie a -> Serie a -> Serie a
suma = zipWith (+)
 
resta :: Num a => Serie a -> Serie a -> Serie a
resta xs ys = suma xs (opuesta ys)
 
producto :: Num a => Serie a -> Serie a -> Serie a
producto (x:xs) zs@(y:ys) = 
    x*y : suma (producto xs zs) (map (x*) ys)
 
cociente :: Fractional a => Serie a -> Serie a -> Serie a
cociente (x:xs) (y:ys) = zs 
    where zs = x/y : map (/y) (resta xs (producto zs ys))  
 
derivada :: (Num a, Enum a) => Serie a -> Serie a
derivada (_:xs) = zipWith (*) xs [1..]
 
integral :: (Fractional a, Enum a) => Serie a -> Serie a
integral xs = 0 : zipWith (/) xs [1..]
 
expx :: Serie Rational
expx = map (1/) (map fromIntegral factoriales)
 
-- factoriales es la lista de los factoriales. Por ejemplo, 
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

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"