Menu Close

Etiqueta: zipWith

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..]

Números de Perrin

Los números de Perrin se definen por la elación de recurrencia

   P(n) = P(n - 2) + P(n - 3) si n > 2,

con los valores iniciales

   P(0) = 3, P(1) = 0 y P(2) = 2.

Definir la sucesión

   sucPerrin :: [Integer]

cuyos elementos son los números de Perrin. Por ejemplo,

   λ> take 15 sucPerrin
   [3,0,2,3,2,5,5,7,10,12,17,22,29,39,51]
   λ> length (show (sucPerrin !! (2*10^5)))
   24425

Comprobar con QuickCheck si se verifica la siguiente propiedad: para todo entero n > 1, el n-ésimo término de la sucesión de Perrin es divisible por n si y sólo si n es primo.

Soluciones

import Data.List (genericIndex, unfoldr)
import Data.Numbers.Primes (isPrime)
import Test.QuickCheck
 
-- 1ª solución
sucPerrin1 :: [Integer]
sucPerrin1 = 3 : 0 : 2 : zipWith (+) sucPerrin1 (tail sucPerrin1)
 
-- 2ª solución
sucPerrin2 :: [Integer]
sucPerrin2 = [x | (x,_,_) <- iterate op (3,0,2)]
  where op (a,b,c) = (b,c,a+b)
 
-- 3ª solución
sucPerrin3 :: [Integer]
sucPerrin3 =
  unfoldr (\(a, (b,c)) -> Just (a, (b,(c,a+b)))) (3,(0,2))
 
-- 4ª solución
sucPerrin4 :: [Integer]
sucPerrin4 = [vectorPerrin n ! n | n <- [0..]]
 
vectorPerrin :: Integer -> Array Integer Integer
vectorPerrin n = v where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 3
  f 1 = 0
  f 2 = 2
  f i = v ! (i-2) + v ! (i-3)
 
-- Comparación de eficiencia
--    λ> length (show (sucPerrin1 !! (3*10^5)))
--    36638
--    (1.62 secs, 2,366,238,984 bytes)
--    λ> length (show (sucPerrin2 !! (3*10^5)))
--    36638
--    (1.40 secs, 2,428,701,384 bytes)
--    λ> length (show (sucPerrin3 !! (3*10^5)))
--    36638
--    (1.14 secs, 2,409,504,864 bytes)
--    λ> length (show (sucPerrin4 !! (3*10^5)))
--    36638
--    (1.78 secs, 2,585,400,776 bytes)
 
 
-- Usaremos la 3ª
sucPerrin :: [Integer]
sucPerrin = sucPerrin3
 
-- La propiedad es  
conjeturaPerrin :: Integer -> Property
conjeturaPerrin n =
  n > 1 ==>
  (perrin n `mod` n == 0) == isPrime n
 
-- (perrin n) es el n-ésimo término de la sucesión de Perrin. Por
-- ejemplo,
--    perrin 4  ==  2
--    perrin 5  ==  5
--    perrin 6  ==  5
perrin :: Integer -> Integer
perrin n = sucPerrin `genericIndex` n
 
-- La comprobación es
--    λ> quickCheck conjeturaPerrin
--    +++ OK, passed 100 tests.
 
-- Nota: Aunque QuickCheck no haya encontrado contraejemplos, la
-- propiedad no es cierta. Sólo lo es una de las implicaciones: si n es
-- primo, entonces el  n-ésimo término de la sucesión de Perrin es
-- divisible por n. La otra es falsa y los primeros contraejemplos son
--    271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291

Problema de las puertas

Un hotel dispone de n habitaciones y n camareros. Los camareros tienen la costumbre de cambiar de estado las puertas (es decir, abrir las cerradas y cerrar las abiertas). El proceso es el siguiente:

  • Inicialmente todas las puertas están cerradas.
  • El primer camarero cambia de estado las puertas de todas las habitaciones.
  • El segundo cambia de estado de las puertas de las habitaciones pares.
  • El tercero cambia de estado todas las puertas que son múltiplos de 3.
  • El cuarto cambia de estado todas las puertas que son múltiplos de 4
  • Así hasta que ha pasado el último camarero.

Por ejemplo, para n = 5

   Pase    | Puerta 1 | Puerta 2 | Puerta 3 | Puerta 4 | Puerta 5
   Inicial | Cerrada  | Cerrada  | Cerrada  | Cerrada  | Cerrada
   Pase 1  | Abierta  | Abierta  | Abierta  | Abierta  | Abierta
   Pase 2  | Abierta  | Cerrada  | Abierta  | Cerrada  | Abierta
   Pase 3  | Abierta  | Cerrada  | Cerrada  | Cerrada  | Abierta
   Pase 4  | Abierta  | Cerrada  | Cerrada  | Abierta  | Abierta
   Pase 5  | Abierta  | Cerrada  | Cerrada  | Abierta  | Cerrada

Los estados de las puertas se representan por el siguiente tipo de datos

   data Estado = Abierta | Cerrada deriving Show

Definir la función

   final :: Int -> [Estado]

tal que (final n) es la lista de los estados de las n puertas después de que hayan pasado los n camareros. Por ejemplo,

   ghci> final 5
   [Abierta,Cerrada,Cerrada,Abierta,Cerrada]
   ghci> final 7
   [Abierta,Cerrada,Cerrada,Abierta,Cerrada,Cerrada,Cerrada]

Soluciones

 
-- 1ª solución
-- ===========
 
data Estado = Abierta | Cerrada 
  deriving (Eq, Show)
 
cambia Abierta = Cerrada
cambia Cerrada = Abierta
 
-- (inicial n) es el estado inicial para el problema de las n
-- habitaciones. Por ejemplo,
--    inicial 5  ==  [Cerrada,Cerrada,Cerrada,Cerrada,Cerrada]
inicial :: Int -> [Estado]
inicial n = replicate n Cerrada
 
-- (pase k es) es la lista de los estados de las puertas después de pasar el
-- camarero k que las encuentra en los estados es. Por ejemplo,
--    ghci> pase 1 (inicial 5)
--    [Abierta,Abierta,Abierta,Abierta,Abierta]
--    ghci> pase 2 it
--    [Abierta,Cerrada,Abierta,Cerrada,Abierta]
--    ghci> pase 3 it
--    [Abierta,Cerrada,Cerrada,Cerrada,Abierta]
--    ghci> pase 4 it
--    [Abierta,Cerrada,Cerrada,Abierta,Abierta]
--    ghci> pase 5 it
--    [Abierta,Cerrada,Cerrada,Abierta,Cerrada]
pase :: [Estado] -> Int -> [Estado] 
pase es k = zipWith cambiaK  es [1..] 
  where cambiaK e n | n `mod` k == 0 = cambia e
                    | otherwise      = e
 
final :: Int -> [Estado]
final n = aux [1..n] (inicial n) 
  where aux []     es = es  
        aux (k:ks) es = aux ks (pase es k)
 
-- 2ª solución
-- ===========
 
final2 :: Int -> [Estado]
final2 n = foldl pase (inicial n) [1..n] 
 
-- 3ª solución
-- =============
 
final3 :: Int -> [Estado]
final3 n = map f [1..n]
  where f x | even (length (divisores x)) = Cerrada
            | otherwise                   = Abierta
 
divisores :: Int -> [Int]
divisores n = [x | x <- [1..n], n `mod` x == 0]
 
-- 4ª solución
-- ===========
 
-- En primer lugar, vamos a determinar la lista de las posiciones
-- (comenzando a contar en 1) de las puertas que quedan abierta en el
-- problema de las n puertas. 
posicionesAbiertas :: Int -> [Int]
posicionesAbiertas n = 
  [x | (x,y) <- zip [1..] (final n), y == Abierta]
 
-- Al calcularlas,
--    ghci> posicionesAbiertas 200
--    [1,4,9,16,25,36,49,64,81,100,121,144,169,196]
-- Se observa las que quedan abiertas son las que sus posiciones son
-- cuadrados perfectos. Usando esta observación se construye la
-- siguiente definición
 
final4 :: Int -> [Estado]
final4 n = aux [1..n] [k*k | k <- [1..]] 
  where aux (x:xs) (y:ys) | x == y  =  Abierta : aux xs ys
        aux (x:xs) ys               =  Cerrada : aux xs ys
        aux []     _                =  []
 
-- ---------------------------------------------------------------------
-- § Comparación de eficiencia                                        --
-- ---------------------------------------------------------------------
 
--    ghci> last (final 1000)
--    Cerrada
--    (0.23 secs, 218727400 bytes)
--    ghci> last (final 2000)
--    Cerrada
--    (1.78 secs, 868883080 bytes)
--    ghci> last (final2 1000)
--    Cerrada
--    (0.08 secs, 218729392 bytes)
--    ghci> last (final2 2000)
--    Cerrada
--    (1.77 secs, 868948600 bytes)
--    ghci> last (final3 1000)
--    Cerrada
--    (0.01 secs, 1029256 bytes)
--    ghci> last (final3 2000)
--    Cerrada
--    (0.01 secs, 2121984 bytes)
--    ghci> last (final4 1000)
--    Cerrada
--    (0.01 secs, 1029328 bytes)
--    ghci> last (final4 2000)
--    Cerrada
--    (0.01 secs, 1578504 bytes)
--    ghci> last (final3 10000)
--    Abierta
--    (0.01 secs, 4670104 bytes)
--    ghci> last (final3 100000)
--    Cerrada
--    (0.09 secs, 38717032 bytes)
--    ghci> last (final3 1000000)
--    Abierta
--    (1.27 secs, 377100832 bytes)
--    ghci> last (final4 1000000)
--    Abierta
--    (1.41 secs, 273292448 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>

Huecos de Aquiles

Un número de Aquiles es un número natural n que es potente (es decir, si p es un divisor primo de n, entonces p² también lo es) y no es una potencia perfecta (es decir, no existen números naturales m y k tales que n es igual a m^k). Por ejemplo,

  • 108 es un número de Aquiles proque es un número potente (ya que su factorización es 2^2 · 3^3, sus divisores primos son 2 and 3 y sus cuadrados (2^2 = 4 y 3^2 = 9) son divisores de 108. Además, 108 no es una potencia perfecta.
  • 360 no es un número de Aquiles ya que 5 es un divisor primo de 360, pero 5^2 = 15 no lo es.
  • 784 no es un número de Aquiles porque, aunque es potente, es una potencia perfecta ya que 784 = 28^2.

Los primeros números de Aquiles son

   72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, ...

Definir las funciones

   esAquiles              :: Integer -> Bool
   huecosDeAquiles        :: [Integer]
   graficaHuecosDeAquiles :: Int -> IO ()

tales que

  • (esAquiles x) se verifica si x es un número de Aquiles. Por ejemplo,
     esAquiles 108         ==  True
     esAquiles 360         ==  False
     esAquiles 784         ==  False
     esAquiles 5425069447  ==  True
     esAquiles 5425069448  ==  True
  • huecosDeAquiles es la sucesión de la diferencias entre los números de Aquiles consecutivos. Por ejemplo,
     λ> take 15 huecosDeAquiles
     [36,92,88,104,40,68,148,27,125,64,104,4,153,27,171]
  • (graficaHuecosDeAquiles n) dibuja la gráfica de los n primeros huecos de Aquiles. Por ejemplo, (graficaHuecosDeAquiles 160) dibuja

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple
 
-- Definición de esAquiles
-- =======================
 
esAquiles :: Integer -> Bool
esAquiles x = esPotente x && noEsPotenciaPerfecta x
 
-- (esPotente x) se verifica si x es potente. Por ejemplo,
--    esPotente 108  ==  True
--    esPotente 360  ==  False
--    esPotente 784  ==  True
esPotente :: Integer -> Bool
esPotente x = all (>1) (exponentes x)
 
-- (exponentes x) es la lista de los exponentes en la factorización de
-- x. Por ejemplo,
--    exponentes 108  ==  [2,3]
--    exponentes 360  ==  [3,2,1]
--    exponentes 784  ==  [4,2]
exponentes :: Integer -> [Int]
exponentes x = map length (group (primeFactors x))
 
-- (noEsPotenciaPerfecta x) se verifica si x no es una potencia
-- perfecta. Por ejemplo,
--    noEsPotenciaPerfecta 108  ==  True
--    noEsPotenciaPerfecta 360  ==  True
--    noEsPotenciaPerfecta 784  ==  False
noEsPotenciaPerfecta :: Integer -> Bool
noEsPotenciaPerfecta x = foldl1 gcd (exponentes x) == 1 
 
-- Definición de huecosDeAquiles
-- =============================
 
huecosDeAquiles :: [Integer]
huecosDeAquiles = zipWith (-) (tail aquiles) aquiles
 
-- aquiles es la sucesión de los números de Aquiles. Por ejemplo, 
--    λ> take 15 aquiles
--    [72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152]
aquiles :: [Integer]
aquiles = filter esAquiles [2..]
 
-- Definición de graficaHuecosDeAquiles
-- ====================================
 
graficaHuecosDeAquiles :: Int -> IO ()
graficaHuecosDeAquiles n =
  plotList [ Key Nothing
           , PNG "Huecos_de_Aquiles.png"
           ]
           (take n huecosDeAquiles)

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>

Término ausente en una progresión aritmética

Una progresión aritmética es una sucesión de números tales que la diferencia de dos términos sucesivos cualesquiera de la sucesión es constante.

Definir la función

   ausente :: Integral a => [a] -> a

tal que (ausente xs) es el único término ausente de la progresión aritmética xs. Por ejemplo,

   ausente [3,7,9,11]               ==  5
   ausente [3,5,9,11]               ==  7
   ausente [3,5,7,11]               ==  9
   ausente ([1..9]++[11..])         ==  10
   ausente ([1..10^6] ++ [2+10^6])  ==  1000001

Nota. Se supone que la lista tiene al menos 3 elementos, que puede ser infinita y que sólo hay un término de la progresión aritmética que no está en la lista.

Soluciones

import Data.List (group, genericLength)
 
-- 1ª solución
ausente :: Integral a => [a] -> a
ausente (x1:xs@(x2:x3:_))
  | d1 == d2     = ausente xs
  | d1 == 2 * d2 = x1 + d2
  | d2 == 2 * d1 = x2 + d1
  where d1 = x2 - x1
        d2 = x3 - x2          
 
-- 2ª solución
ausente2 :: Integral a => [a] -> a
ausente2 s@(x1:x2:x3:xs) 
  | x1 + x3 /= 2 * x2 = x1 + (x3 - x2)
  | otherwise         = head [a | (a,b) <- zip [x1,x2..] s
                                , a /= b]
 
-- 3ª solución
ausente3 :: Integral a => [a] -> a
ausente3  xs@(x1:x2:_) 
  | null us   = x1 + v
  | otherwise = x2 + u * genericLength (u:us) 
  where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)
 
-- Comparación de eficiencia
--    ghci> let n = 10^6 in ausente1 ([1..n] ++ [n+2])
--    1000001
--    (3.53 secs, 634729880 bytes)
--    
--    ghci> let n = 10^6 in ausente2 ([1..n] ++ [n+2])
--    1000001
--    (0.86 secs, 346910784 bytes)
--    
--    ghci> let n = 10^6 in ausente3 ([1..n] ++ [n+2])
--    1000001
--    (1.22 secs, 501521888 bytes)
--    
--    ghci> let n = 10^7 in ausente2 ([1..n] ++ [n+2])
--    10000001
--    (8.68 secs, 3444142568 bytes)
--    
--    ghci> let n = 10^7 in ausente3 ([1..n] ++ [n+2])
--    10000001
--    (12.59 secs, 4975932088 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>

Cálculo de pi mediante los métodos de Gregory-Leibniz y de Beeler

La fórmula de Gregory-Leibniz para calcular pi es
Calculo_de_pi_mediante_los_metodos_de_Gregory-Leibniz_y_de_Beeler_1
y la de Beeler es
Calculo_de_pi_mediante_los_metodos_de_Gregory-Leibniz_y_de_Beeler_2

Definir las funciones

   aproximaPiGL     :: Int -> Double
   aproximaPiBeeler :: Int -> Double
   graficas         :: [Int] -> IO ()

tales que

  • (aproximaPiGL n) es la aproximación de pi con los primeros n términos de la fórmula de Gregory-Leibniz. Por ejemplo,
     aproximaPiGL 1       ==  4.0
     aproximaPiGL 2       ==  2.666666666666667
     aproximaPiGL 3       ==  3.466666666666667
     aproximaPiGL 10      ==  3.0418396189294032
     aproximaPiGL 100     ==  3.1315929035585537
     aproximaPiGL 1000    ==  3.140592653839794
     aproximaPiGL 10000   ==  3.1414926535900345
     aproximaPiGL 100000  ==  3.1415826535897198
  • (aproximaPiBeeler n) es la aproximación de pi con los primeros n términos de la fórmula de Beeler. Por ejemplo,
     aproximaPiBeeler 1   ==  2.0
     aproximaPiBeeler 2   ==  2.6666666666666665
     aproximaPiBeeler 3   ==  2.933333333333333
     aproximaPiBeeler 10  ==  3.140578169680337
     aproximaPiBeeler 60  ==  3.141592653589793
     pi                   ==  3.141592653589793
  • (graficas xs) dibuja la gráfica de las k-ésimas aproximaciones de pi, donde k toma los valores de la lista xs, con las fórmulas de Gregory-Leibniz y de Beeler. Por ejemplo, (graficas [1..25]) dibuja
    Calculo_de_pi_mediante_los_metodos_de_Gregory-Leibniz_y_de_Beeler_3
    donde la línea morada corresponde a la aproximación de Gregory-Leibniz y la verde a la de Beeler.

Soluciones

import Graphics.Gnuplot.Simple
 
-- Definiciones de aproximaPiGL
-- ============================
 
-- 1ª definición de aproximaPiGL
aproximaPiGL :: Int -> Double
aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]
  where sumaA (x:y:xs) = x:(-y):sumaA xs
 
-- 2ª definición de aproximaPiGL
aproximaPiGL2 :: Int -> Double
aproximaPiGL2 n =
  4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))
 
-- 3ª definición de aproximaPiGL
aproximaPiGL3 :: Int -> Double
aproximaPiGL3 n =
  4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]
 
-- 4ª definición de aproximaPiGL
aproximaPiGL4 :: Int -> Double
aproximaPiGL4 n = serieGL !! (n-1)
 
serieGL :: [Double]
serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)
  where numeradores   = cycle [4,-4]
        denominadores = [1,3..]
 
-- Definición de aproximaPiBeeler
aproximaPiBeeler :: Int -> Double
aproximaPiBeeler n = 2 * aux (fromIntegral n) 1
  where
    aux :: Double -> Double -> Double 
    aux n k | n == k    = 1
            | otherwise = 1 + (k/(2*k+1)) * aux n (1+k)
 
-- Definición de graficas
graficas :: [Int] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
             [[(k,aproximaPiGL k)     | k <- xs],
              [(k,aproximaPiBeeler k) | k <- xs]]

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>

Cálculo de pi mediante el método de Newton

El método de Newton para el cálculo de pi se basa en la relación
Calculo_de_pi_mediante_el_metodo_de_Newton_1
y en el desarrollo del arco seno
Calculo_de_pi_mediante_el_metodo_de_Newton_2
de donde se obtiene la fórmula
Calculo_de_pi_mediante_el_metodo_de_Newton_3

La primeras aproximaciones son

   a(0) = 6*(1/2)                               = 3.0
   a(1) = 6*(1/2+1/(2*3*2^3))                   = 3.125
   a(2) = 6*(1/2+1/(2*3*2^3)+(1*3)/(2*4*5*2^5)) = 3.1390625

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 Newton. Por ejemplo,
     aproximacionPi 0   ==  3.0
     aproximacionPi 1   ==  3.125
     aproximacionPi 2   ==  3.1390625
     aproximacionPi 10  ==  3.1415926468755613
     aproximacionPi 21  ==  3.141592653589793
     pi                 ==  3.141592653589793
  • (grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi donde k toma los valores de la lista xs. Por ejemplo, (grafica [1..30]) dibuja
    Calculo_de_pi_mediante_el_metodo_de_Newton_4

Soluciones

import Graphics.Gnuplot.Simple
 
-- 1ª definición
-- =============
 
aproximacionPi :: Int -> Double
aproximacionPi n = 6 * arcsinX
  where arcsinX = 0.5 + sum (take n factoresN)
 
factoresN :: [Double]
factoresN = zipWith (*) (potenciasK 3) fraccionesPI
 
potenciasK :: Double -> [Double]
potenciasK k = (0.5**k)/k : potenciasK (k+2)
 
fraccionesPI :: [Double]
fraccionesPI =
  scanl (*) (1/2) (tail (zipWith (/) [1,3..] [2,4..]))
 
-- 2ª definición
-- =============
 
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = 6 * (serie !! n)
 
serie :: [Double]
serie = scanl1 (+) (zipWith (/)
                            (map fromIntegral numeradores)
                            (map fromIntegral denominadores))
  where numeradores    = 1 : scanl1 (*) [1,3..]
        denominadores  = zipWith (*) denominadores1 denominadores2
        denominadores1 = 2 : scanl1 (*) [2,4..]
        denominadores2 = 1 : [n * 2^n | n <- [3,5..]]
 
grafica :: [Int] -> IO ()
grafica xs = 
    plotList [Key Nothing]
             [(k,aproximacionPi k) | k <- xs]

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>

Pensamiento

“Mi trabajo siempre trató de unir lo verdadero con lo bello; pero cuando tuve que elegir uno u otro, generalmente elegí lo bello.”

Hermann Weyl.

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)
 
-- 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"

Triángulo de Bell

El triágulo de Bell es el triángulo numérico, cuya primera fila es [1] y en cada fila, el primer elemento es el último de la fila anterior y el elemento en la posición j se obtiene sumando el elemento anterior de su misma fila y de la fila anterior. Sus primeras filas son

   1 
   1   2
   2   3   5
   5   7  10  15
   15 20  27  37  52
   52 67  87 114 151 203

Definir la función

   trianguloDeBell :: [[Integer]]

tal que trianguloDeBell es la lista con las filas de dicho triángulo. Por ejemplo

   λ> take 5 trianguloDeBell
   [[1],[1,2],[2,3,5],[5,7,10,15],[15,20,27,37,52]]

Comprobar con QuickCheck que los números que aparecen en la primera columna del triángulo coinciden con los números de Bell; es decir, el primer elemento de la n-ésima fila es el n-ésimo número de Bell.

Soluciones

import Data.List (genericIndex, genericLength)
import Test.QuickCheck
 
trianguloDeBell :: [[Integer]]
trianguloDeBell = iterate siguiente [1]
 
-- (siguiente xs) es la fila siguiente de xs en el triángulo de
-- Bell. Por ejemplo,
--    siguiente [1]     ==  [1,2]
--    siguiente [1,2]   ==  [2,3,5]
--    siguiente [2,3,5] ==  [5,7,10,15]
siguiente :: [Integer] -> [Integer]
siguiente xs = last xs : zipWith (+) xs (siguiente xs)
 
-- Propiedad
-- =========
 
-- La propiedad es
prop_TrianguloDeBell :: Integer -> Property
prop_TrianguloDeBell n =
  n > 0 ==> head (trianguloDeBell `genericIndex` n) == bell n
 
-- (bell n) es el n-ésimo número de Bell definido en el ejercicio
-- anterior.  
bell :: Integer -> Integer
bell n = genericLength (particiones [1..n])
 
particiones :: [a] -> [[[a]]]
particiones [] = [[]]
particiones (x:xs) =
  concat [([x] : yss) : inserta x yss | yss <- ysss]
  where ysss = particiones xs
 
inserta :: a -> [[a]] -> [[[a]]]
inserta _ []       = []
inserta x (ys:yss) = ((x:ys):yss) : [ys : zs | zs <- inserta x yss] 
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_TrianguloDeBell
--    +++ OK, passed 100 tests.

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>

Pensamiento

“La ciencia es lo que entendemos lo suficientemente bien como para explicarle a una computadora. El arte es todo lo demás.”

Donald Knuth.

Codificación de Gödel

El Consejo Ejecutivo de la UNESCO ha proclamado el 14 de enero como Día mundial de la Lógica.

La fecha del 14 de enero se ha eligido para rendir homenaje a dos grandes lógicos del siglo XX: Kurt Gödel (que falleció el 14 de enero de 1978) y Alfred Tarski (que nació el 14 de enero de 1901).

Gödel demostró en 1931 los teoremas de incompletitud y para su demostración introdujo la actualmente conocida como codificación de Gödel que asigna a cada fórmula de un lenguaje formal un número único.

Dada una lista de números naturales xs, la codificación de Gödel de xs se obtiene multiplicando las potencias de los primos sucesivos, siendo los exponentes los sucesores de los elementos de xs. Por ejemplo, si xs es [6,0,4], la codificación de xs es

   2^(6+1) * 3^(0+1) * 5^(4+1) = 1200000

Definir las funciones

   codificaG   :: [Integer] -> Integer
   decodificaG :: Integer -> [Integer]

tales que

  • (codificaG xs) es la codificación de Gödel de xs. Por ejemplo,
     codificaG [6,0,4]            ==  1200000
     codificaG [3,1,1]            ==  3600
     codificaG [3,1,0,0,0,0,0,1]  ==  4423058640
     codificaG [1..6]             ==  126111168580452537982500
  • (decodificaG n) es la lista xs cuya codificación es n. Por ejemplo,
     decodificaG 1200000                   ==  [6,0,4]
     decodificaG 3600                      ==  [3,1,1]
     decodificaG 4423058640                ==  [3,1,0,0,0,0,0,1]
     decodificaG 126111168580452537982500  ==  [1,2,3,4,5,6]

Comprobar con QuickCheck que decodificaG es la inversa por la izquierda de codificaG; es decir, para toda lista xs de números enteros, se verifica que

   decodificaG (codificaG ys) == ys

donde ys es la lista de los valores absolutos de los elementos de xs.

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primes, primeFactors)
import Test.QuickCheck
 
codificaG   :: [Integer] -> Integer
codificaG xs = product [p^(x+1) | (p,x) <- zip primes xs] 
 
decodificaG :: Integer -> [Integer]
decodificaG n =
  [genericLength xs - 1 | xs <- group (primeFactors n)]
 
-- La propiedad es
propCodifica1 :: [Integer] -> Bool
propCodifica1 xs =
  decodificaG (codificaG ys) == ys
  where ys = map abs xs
 
-- La comprobación es
--    ghci> quickCheck propCodifica1
--    +++ OK, passed 100 tests.

Pensamiento

La verdad es la verdad, dígala Agamenón o su porquero.
– Agamenón: Conforme.
– El porquero: No me convence.

Antonio Machado