Menu Close

Etiqueta: scanl1

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

La sucesión de Sylvester

La sucesión de Sylvester es la sucesión que comienza en 2 y sus restantes términos se obtienen multiplicando los anteriores y sumándole 1.

Definir las funciones

   sylvester        :: Integer -> Integer
   graficaSylvester :: Integer -> Integer -> IO ()

tales que

  • (sylvester n) es el n-ésimo término de la sucesión de Sylvester. Por ejemplo,
     λ> [sylvester n | n <- [0..7]]
     [2,3,7,43,1807,3263443,10650056950807,113423713055421844361000443]
     λ> length (show (sylvester 25))
     6830085
  • (graficaSylvester d n) dibuja la gráfica de los d últimos dígitos de los n primeros términos de la sucesión de Sylvester. Por ejemplo,
    • (graficaSylvester 3 30) dibuja
      La_sucesion_de_Sylvester_(3,30)
    • (graficaSylvester 4 30) dibuja
      La_sucesion_de_Sylvester_(4,30)
    • (graficaSylvester 5 30) dibuja
      La_sucesion_de_Sylvester_(5,30)

Nota: Se puede usar programación dinámica para aumentar la eficiencia.

Soluciones

import Data.List               (genericIndex)
import Data.Array              ((!), array)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG))
 
-- 1ª solución (por recursión)
-- ===========================
 
sylvester1 :: Integer -> Integer
sylvester1 0 = 2
sylvester1 n = 1 + product [sylvester1 k | k <- [0..n-1]]
 
-- 2ª solución (con programación dinámica)
-- =======================================
 
sylvester2 :: Integer -> Integer
sylvester2 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + product [v!k | k <- [0..m-1]]
 
-- 3ª solución
-- ===========
 
-- Observando que
--    S(n) = 1 + S(0)*S(1)*...*S(n-2)*S(n-1)
--         = 1 + (1 + S(0)*S(1)*...*S(n-2))*S(n-1) - S(n-1)
--         = 1 + S(n-1)*S(n-1) - S(n-1)
--         = 1 + S(n-1)^2 - S(n-1)
-- se obtiene la siguiente definición.
sylvester3 :: Integer -> Integer
sylvester3 0 = 2
sylvester3 n = 1 + x^2 - x
  where x = sylvester3 (n-1)
 
-- 4ª solución
-- ===========
 
sylvester4 :: Integer -> Integer
sylvester4 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + x^2 - x
    where x = v ! (m-1)
 
-- 5ª solución
-- ===========
 
sylvester5 :: Integer -> Integer
sylvester5 n = sucSylvester5 `genericIndex` n
 
sucSylvester5 :: [Integer]
sucSylvester5 = iterate (\x -> (x-1)*x+1) 2 
 
-- La comparación es
--    λ> length (show (sylvester1 23))
--    1707522
--    (6.03 secs, 4,090,415,704 bytes)
--    λ> length (show (sylvester2 23))
--    1707522
--    (0.33 secs, 109,477,296 bytes)
--    λ> length (show (sylvester3 23))
--    1707522
--    (0.35 secs, 109,395,136 bytes)
--    λ> length (show (sylvester4 23))
--    1707522
--    (0.33 secs, 109,402,440 bytes)
--    λ> length (show (sylvester5 23))
--    1707522
--    (0.30 secs, 108,676,256 bytes)
 
graficaSylvester :: Integer -> Integer -> IO ()
graficaSylvester d n =
  plotList [ Key Nothing
           , PNG ("La_sucesion_de_Sylvester_" ++ show (d,n) ++ ".png")
           ]
           [sylvester5 k `mod` (10^d) | k <- [0..n]]

Suma de segmentos iniciales

Los segmentos iniciales de [3,1,2,5] son [3], [3,1], [3,1,2] y [3,1,2,5]. Sus sumas son 3, 4, 6 y 9, respectivamente. La suma de dichas sumas es 24.

Definir la función

   sumaSegmentosIniciales :: [Integer] -> Integer

tal que (sumaSegmentosIniciales xs) es la suma de las sumas de los segmentos iniciales de xs. Por ejemplo,

   sumaSegmentosIniciales [3,1,2,5]     ==  24
   sumaSegmentosIniciales [1..3*10^6]  ==  4500004500001000000

Comprobar con QuickCheck que la suma de las sumas de los segmentos iniciales de la lista formada por n veces el número uno es el n-ésimo número triangular; es decir que

   sumaSegmentosIniciales (genericReplicate n 1)

es igual a

   n * (n + 1) `div` 2

Soluciones

import Data.List (genericLength, genericReplicate)
import Test.QuickCheck
 
-- 1ª solución
-- ===========
 
sumaSegmentosIniciales :: [Integer] -> Integer
sumaSegmentosIniciales xs =
  sum [sum (take k xs) | k <- [1.. length xs]]
 
-- 2ª solución
-- ===========
 
sumaSegmentosIniciales2 :: [Integer] -> Integer
sumaSegmentosIniciales2 xs =
  sum (zipWith (*) [n,n-1..1] xs)
  where n = genericLength xs
 
-- 3ª solución
-- ===========
 
sumaSegmentosIniciales3 :: [Integer] -> Integer
sumaSegmentosIniciales3 xs =
  sum (scanl1 (+) xs)
 
-- Comprobación de la equivalencia
-- ===============================
 
-- La propiedad es
prop_sumaSegmentosInicialesEquiv :: [Integer] -> Bool
prop_sumaSegmentosInicialesEquiv xs =
  all (== sumaSegmentosIniciales xs) [f xs | f <- [ sumaSegmentosIniciales2
                                                  , sumaSegmentosIniciales3]]
 
-- La comprobación es
--   λ> quickCheck prop_sumaSegmentosInicialesEquiv
--   +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
--   λ> sumaSegmentosIniciales [1..10^4]
--   166716670000
--   (2.42 secs, 7,377,926,824 bytes)
--   λ> sumaSegmentosIniciales2 [1..10^4]
--   166716670000
--   (0.01 secs, 4,855,176 bytes)
--   
--   λ> sumaSegmentosIniciales2 [1..3*10^6]
--   4500004500001000000
--   (2.68 secs, 1,424,404,168 bytes)
--   λ> sumaSegmentosIniciales3 [1..3*10^6]
--   4500004500001000000
--   (1.54 secs, 943,500,384 bytes)
 
-- Comprobación de la propiedad
-- ============================
 
-- La propiedad es
prop_sumaSegmentosIniciales :: Positive Integer -> Bool
prop_sumaSegmentosIniciales (Positive n) =
  sumaSegmentosIniciales3 (genericReplicate n 1) ==
  n * (n + 1) `div` 2
 
-- La compronación es
--   λ> quickCheck prop_sumaSegmentosIniciales
--   +++ 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>

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

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

“¿Por qué son hermosos los números? Es como preguntar por qué es bella la Novena Sinfonía de Beethoven. Si no ves por qué, alguien no puede decírtelo. Yo sé que los números son hermosos. Si no son hermosos, nada lo es.”

Paul Erdös.