zipWith – Exercitium

Emparejamiento de árboles

PorJosé A. Alonso 18-abril-202230-abril-2022

Los árboles se pueden representar mediante el siguiente tipo de datos

   data Arbol a = N a [Arbol a]
     deriving (Show, Eq)

1 2	data Arbol a = N a [Arbol a] deriving (Show, Eq)

Por ejemplo, los árboles

     1               3
    / \             /|\
   6   3           / | \
       |          5  4  7
       5          |     /\
                  6    2  1

1 3

/ \ /|\

6 3 / | \

| 5 4 7

5 | /\

6 2 1

se representan por

   ej1, ej2 :: Arbol Int
   ej1 = N 1 [N 6 [],N 3 [N 5 []]]
   ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

ej1, ej2 :: Arbol Int

ej1 = N 1 [N 6 [],N 3 [N 5 []]]

ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

Definir la función

   emparejaArboles :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c

1	emparejaArboles :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c

tal que (emparejaArboles f a1 a2) es el árbol obtenido aplicando la función f a los elementos de los árboles a1 y a2 que se encuentran en la misma posición. Por ejemplo,

   λ> emparejaArboles (+) (N 1 [N 2 [], N 3[]]) (N 1 [N 6 []])
   N 2 [N 8 []]
   λ> emparejaArboles (+) ej1 ej2
   N 4 [N 11 [],N 7 []]
   λ> emparejaArboles (+) ej1 ej1
   N 2 [N 12 [],N 6 [N 10 []]]

λ> emparejaArboles (+) (N 1 [N 2 [], N 3[]]) (N 1 [N 6 []])

N 2 [N 8 []]

λ> emparejaArboles (+) ej1 ej2

N 4 [N 11 [],N 7 []]

λ> emparejaArboles (+) ej1 ej1

N 2 [N 12 [],N 6 [N 10 []]]

Soluciones

module Emparejamiento_de_arboles where

import Data.Tree (Tree (..))
import Control.Monad.Zip (mzipWith)
import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, sized, quickCheck)

data Arbol a = N a [Arbol a]
  deriving (Show, Eq)

ej1, ej2 :: Arbol Int
ej1 = N 1 [N 6 [],N 3 [N 5 []]]
ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

-- 1ª solución
-- ===========

emparejaArboles1 :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c
emparejaArboles1 f (N x xs) (N y ys) =
  N (f x y) (zipWith (emparejaArboles1 f) xs ys)

-- 2ª solución
-- ===========

emparejaArboles2 :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c
emparejaArboles2 f x y =
  treeAarbol (mzipWith f (arbolAtree x) (arbolAtree y))

arbolAtree :: Arbol a -> Tree a
arbolAtree (N x xs) = Node x (map arbolAtree xs)

treeAarbol :: Tree a -> Arbol a
treeAarbol (Node x xs) = N x (map treeAarbol xs)

-- Comprobación de equivalencia
-- ============================

-- (arbolArbitrario n) es un árbol aleatorio de orden n. Por ejemplo,
--    λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))
--    N (-26) [N 8 [N 6 [N 11 []]],N 7 []]
--    λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))
--    N 1 []
--    λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))
--    N (-19) [N (-11) [],N 25 [],N 19 [N (-27) [],N (-19) [N 17 []]]]
arbolArbitrario :: Arbitrary a => Int -> Gen (Arbol a)
arbolArbitrario n = do
  x  <- arbitrary
  ms <- sublistOf [0 .. n `div` 2]
  as <- mapM arbolArbitrario ms
  return (N x as)

-- Arbol es una subclase de Arbitraria
instance Arbitrary a => Arbitrary (Arbol a) where
  arbitrary = sized arbolArbitrario

-- La propiedad es
prop_emparejaArboles :: Arbol Int -> Arbol Int -> Bool
prop_emparejaArboles x y =
  emparejaArboles1 (+) x y == emparejaArboles2 (+) x y &&
  emparejaArboles1 (*) x y == emparejaArboles2 (*) x y

-- La comprobación es
--    λ> quickCheck prop_emparejaArboles
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> a500 <- generate (arbolArbitrario 500 :: Gen (Arbol Int))
--    λ> emparejaArboles1 (+) a500 a500 == emparejaArboles1 (+) a500 a500
--    True
--    (1.92 secs, 1,115,813,352 bytes)
--    λ> emparejaArboles2 (+) a500 a500 == emparejaArboles2 (+) a500 a500
--    True
--    (3.28 secs, 2,212,257,928 bytes)
--
--    λ> b500 = arbolAtree a500
--    λ> mzipWith (+) b500 b500 == mzipWith (+) b500 b500
--    True
--    (0.21 secs, 563,503,112 bytes)

module Emparejamiento_de_arboles where

import Data.Tree (Tree (..))

import Control.Monad.Zip (mzipWith)

import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, sized, quickCheck)

data Arbol a = N a [Arbol a]

deriving (Show, Eq)

ej1, ej2 :: Arbol Int

ej1 = N 1 [N 6 [],N 3 [N 5 []]]

ej2 = N 3 [N 5 [N 6 []], N 4 [], N 7 [N 2 [], N 1 []]]

-- 1ª solución

-- ===========

emparejaArboles1 :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c

emparejaArboles1 f (N x xs) (N y ys) =

N (f x y) (zipWith (emparejaArboles1 f) xs ys)

-- 2ª solución

-- ===========

emparejaArboles2 :: (a -> b -> c) -> Arbol a -> Arbol b -> Arbol c

emparejaArboles2 f x y =

treeAarbol (mzipWith f (arbolAtree x) (arbolAtree y))

arbolAtree :: Arbol a -> Tree a

arbolAtree (N x xs) = Node x (map arbolAtree xs)

treeAarbol :: Tree a -> Arbol a

treeAarbol (Node x xs) = N x (map treeAarbol xs)

-- Comprobación de equivalencia

-- ============================

-- (arbolArbitrario n) es un árbol aleatorio de orden n. Por ejemplo,

-- λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))

-- N (-26) [N 8 [N 6 [N 11 []]],N 7 []]

-- λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))

-- N 1 []

-- λ> generate (arbolArbitrario 5 :: Gen (Arbol Int))

-- N (-19) [N (-11) [],N 25 [],N 19 [N (-27) [],N (-19) [N 17 []]]]

arbolArbitrario :: Arbitrary a => Int -> Gen (Arbol a)

arbolArbitrario n = do

x <- arbitrary

ms <- sublistOf [0 .. n `div` 2]

as <- mapM arbolArbitrario ms

return (N x as)

-- Arbol es una subclase de Arbitraria

instance Arbitrary a => Arbitrary (Arbol a) where

arbitrary = sized arbolArbitrario

-- La propiedad es

prop_emparejaArboles :: Arbol Int -> Arbol Int -> Bool

prop_emparejaArboles x y =

emparejaArboles1 (+) x y == emparejaArboles2 (+) x y &&

emparejaArboles1 (*) x y == emparejaArboles2 (*) x y

-- La comprobación es

-- λ> quickCheck prop_emparejaArboles

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> a500 <- generate (arbolArbitrario 500 :: Gen (Arbol Int))

-- λ> emparejaArboles1 (+) a500 a500 == emparejaArboles1 (+) a500 a500

-- True

-- (1.92 secs, 1,115,813,352 bytes)

-- λ> emparejaArboles2 (+) a500 a500 == emparejaArboles2 (+) a500 a500

-- True

-- (3.28 secs, 2,212,257,928 bytes)

-- λ> b500 = arbolAtree a500

-- λ> mzipWith (+) b500 b500 == mzipWith (+) b500 b500

-- True

-- (0.21 secs, 563,503,112 bytes)

El código se encuentra en GitHub.

La elaboración de las soluciones se describe en el siguiente vídeo

Nuevas soluciones

En los comentarios se pueden escribir nuevas soluciones.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

Operaciones con series de potencias

PorJosé A. Alonso 29-mayo-202025-marzo-2022

Una serie de potencias es una serie de la forma

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

1	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! + ...

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

1	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

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]

1 2	λ> 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]

1 2	λ> 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]

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

1 2	λ> 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]

1 2	λ> 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]

1 2	λ> 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]

1 2	λ> 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]

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

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

PorJosé A. Alonso 19-mayo-202026-mayo-2020

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

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

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

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

Definir la sucesión

   sucPerrin :: [Integer]

1	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

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

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

PorJosé A. Alonso 28-abril-20205-mayo-2020

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

1	data Estado = Abierta \| Cerrada deriving Show

Definir la función

   final :: Int -> [Estado]

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

ghci> final 5

[Abierta,Cerrada,Cerrada,Abierta,Cerrada]

ghci> final 7

[Abierta,Cerrada,Cerrada,Abierta,Cerrada,Cerrada,Cerrada]

Soluciones

[schedule on=’2020-06-05′ at=»06:00″]

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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

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

PorJosé A. Alonso 23-abril-202026-marzo-2022

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

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

Definir las funciones

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

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

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]

1 2	λ> 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)

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>

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

PorJosé A. Alonso 7-abril-202026-marzo-2022

La fórmula de Gregory-Leibniz para calcular pi es

y la de Beeler es

Definir las funciones

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

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

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

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

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

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

PorJosé A. Alonso 18-marzo-202025-marzo-2020

El método de Newton para el cálculo de pi se basa en la relación

y en el desarrollo del arco seno

de donde se obtiene la fórmula

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

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

1 2	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

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

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]

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

PorJosé A. Alonso 12-marzo-202025-marzo-2022

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

1 2	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

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]

1	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 |
+------+----------------+----------------+

+------+----------------+----------------+

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

1	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 |
+------+----------------+----------------+

+------+----------------+----------------+

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

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

PorJosé A. Alonso 6-febrero-202013-febrero-2020

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

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

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

1 2	λ> 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.

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.

Suma de números de Fibonacci con índice impar

PorJosé A. Alonso 23-diciembre-201930-diciembre-2019

La sucesión de Fibonacci, F(n), es la siguiente sucesión infinita de números naturales:

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

La sucesión comienza con los números 0 y 1. A partir de estos, cada término es la suma de los dos anteriores.

Definir la función

   sumaFibsIndiceImpar :: Int -> Integer

1	sumaFibsIndiceImpar :: Int -> Integer

tal que (sumaFibsIndiceImpar n) es la suma de los n primeros términos de la sucesión de Fibonacci no índice impar; es decir,

   sumaFibsIndiceImpar n = F(1) + F(3) + ... + F(2*n-1)

1	sumaFibsIndiceImpar n = F(1) + F(3) + ... + F(2*n-1)

Por ejemplo,

   sumaFibsIndiceImpar 1  ==  1
   sumaFibsIndiceImpar 2  ==  3
   sumaFibsIndiceImpar 3  ==  8
   sumaFibsIndiceImpar 4  ==  21
   sumaFibsIndiceImpar 5  ==  55
   sumaFibsIndiceImpar (10^4) `rem` (10^9)  ==  213093125

sumaFibsIndiceImpar 1 == 1

sumaFibsIndiceImpar 2 == 3

sumaFibsIndiceImpar 3 == 8

sumaFibsIndiceImpar 4 == 21

sumaFibsIndiceImpar 5 == 55

sumaFibsIndiceImpar (10^4) `rem` (10^9) == 213093125

En los ejemplos anteriores se observa que

   sumaFibsIndiceImpar 1  ==  F(2)
   sumaFibsIndiceImpar 2  ==  F(4)
   sumaFibsIndiceImpar 3  ==  F(6)
   sumaFibsIndiceImpar 4  ==  F(8)
   sumaFibsIndiceImpar 5  ==  F(10)

sumaFibsIndiceImpar 1 == F(2)

sumaFibsIndiceImpar 2 == F(4)

sumaFibsIndiceImpar 3 == F(6)

sumaFibsIndiceImpar 4 == F(8)

sumaFibsIndiceImpar 5 == F(10)

Comprobar con QuickCheck que (sumaFibsIndiceImpar n) es F(2n); es decir, el 2n-ésimo número de Fibonacci

Soluciones

import Test.QuickCheck

-- 1ª solución
-- ===========

sumaFibsIndiceImpar :: Int -> Integer
sumaFibsIndiceImpar n = sum [fib (2*k-1) | k <- [1..n]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por
-- ejemplo,
--    fib 6  ==  8
fib :: Int -> Integer
fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,
--    λ> take 20 fibs
--    [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

-- 2ª solución
-- ============

sumaFibsIndiceImpar2 :: Int -> Integer
sumaFibsIndiceImpar2 n =
  sum [a | (a,b) <- zip fibs [0..2*n], odd b]

-- Comparación de eficiencia
-- =========================

--    λ> sumaFibsIndiceImpar (10^4) `rem` (10^9)
--    213093125
--    (0.98 secs, 13,889,312 bytes)
--    λ> sumaFibsIndiceImpar2 (10^4) `rem` (10^9)
--    213093125
--    (0.05 secs, 18,047,720 bytes)

-- Comprobación
-- ============

-- La propiedad es
prop_sumaFibsIndiceImpar :: Int -> Property
prop_sumaFibsIndiceImpar n =
  n >= 0 ==> sumaFibsIndiceImpar n == fib (2*n)

-- La comprobación es
--    λ> quickCheck prop_sumaFibsIndiceImpar
--    +++ OK, passed 100 tests.

import Test.QuickCheck

-- 1ª solución

-- ===========

sumaFibsIndiceImpar :: Int -> Integer

sumaFibsIndiceImpar n = sum [fib (2*k-1) | k <- [1..n]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por

-- ejemplo,

-- fib 6 == 8

fib :: Int -> Integer

fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,

-- λ> take 20 fibs

-- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]

fibs :: [Integer]

fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

-- 2ª solución

-- ============

sumaFibsIndiceImpar2 :: Int -> Integer

sumaFibsIndiceImpar2 n =

sum [a | (a,b) <- zip fibs [0..2*n], odd b]

-- Comparación de eficiencia

-- =========================

-- λ> sumaFibsIndiceImpar (10^4) `rem` (10^9)

-- 213093125

-- (0.98 secs, 13,889,312 bytes)

-- λ> sumaFibsIndiceImpar2 (10^4) `rem` (10^9)

-- 213093125

-- (0.05 secs, 18,047,720 bytes)

-- Comprobación

-- ============

-- La propiedad es

prop_sumaFibsIndiceImpar :: Int -> Property

prop_sumaFibsIndiceImpar n =

n >= 0 ==> sumaFibsIndiceImpar n == fib (2*n)

-- La comprobación es

-- λ> quickCheck prop_sumaFibsIndiceImpar

-- +++ OK, passed 100 tests.

Referencia

Sum of sequence of odd index Fibonacci numbers en ProofWiki.

Pensamiento

El corazón del poeta, tan rico en sonoridades, es casi un insulto a la afonía cordial de la masa.

Antonio Machado

Operaciones con series de potencias

PorJosé A. Alonso 26-abril-201925-marzo-2022

Una serie de potencias es una serie de la forma

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

1	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! + ...

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

1	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

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]

1 2	λ> 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]

1 2	λ> 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]

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

1 2	λ> 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]

1 2	λ> 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]

1 2	λ> 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]

1 2	λ> 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]

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

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

Pensamiento

Ni mármol duro y eterno,
ni música ni pintura,
sino palabra en el tiempo.

Antonio Machado

Problema de las puertas

PorJosé A. Alonso 24-abril-201925-abril-2021

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

1	data Estado = Abierta \| Cerrada deriving Show

Definir la función

   final :: Int -> [Estado]

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

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)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

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

Pensamiento

… cuánto exilio en la presencia cabe.

Antonio Machado

Distribución de diferencias de dígitos consecutivos de pi

PorJosé A. Alonso 4-abril-201911-abril-2019

Usando la librería Data.Number.CReal, que se instala con

   cabal install number

1	cabal install number

se pueden calcular el número pi con la precisión que se desee. Por ejemplo,

   λ> import Data.Number.CReal
   λ> showCReal 60 pi
   "3.141592653589793238462643383279502884197169399375105820974945"

λ> import Data.Number.CReal

λ> showCReal 60 pi

"3.141592653589793238462643383279502884197169399375105820974945"

importa la librería y calcula el número pi con 60 decimales.

La distribución de las diferencias de los dígitos consecutivos para los 18 primeros n dígitos de pi se calcula como sigue: los primeros 18 dígitos de pi son

   3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

1	3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

Las diferencias de sus elementos consecutivos es

   2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

1	2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

y la distribución de sus frecuencias en el intervalo [-9,9] es

   0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

1	0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

es decir, el desde el -9 a -5 no aparecen, el -4 aparece 3 veces, el -2 aparece 2 veces y así sucesivamente.

Definir las funciones

   distribucionDDCpi :: Int -> [Int]
   graficas :: [Int] -> FilePath -> IO ()

1 2	distribucionDDCpi :: Int -> [Int] graficas :: [Int] -> FilePath -> IO ()

tales que

(distribucionDDCpi n) es la distribución de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi. Por ejemplo,

     λ> distribucionDDCpi 18
     [0,0,0,0,0,3,2,2,2,0,2,3,1,0,0,1,1,0,0]
     λ> distribucionDDCpi 100
     [1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]
     λ> distribucionDDCpi 200
     [3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]
     λ> distribucionDDCpi 1000
     [16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]
     λ> distribucionDDCpi 5000
     [67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

λ> distribucionDDCpi 18

[0,0,0,0,0,3,2,2,2,0,2,3,1,0,0,1,1,0,0]

λ> distribucionDDCpi 100

[1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]

λ> distribucionDDCpi 200

[3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]

λ> distribucionDDCpi 1000

[16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]

λ> distribucionDDCpi 5000

[67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

(graficas ns f) dibuja en el fichero f las gráficas de las distribuciones de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi, para n en ns. Por ejemplo, al evaluar (graficas [100,250..4000] «distribucionDDCpi.png» se escribe en el fichero «distribucionDDCpi.png» la siguiente gráfica

Soluciones

import Data.Number.CReal
import Graphics.Gnuplot.Simple
import Data.Array

--    λ> digitosPi 18
--    [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]
digitosPi :: Int -> [Int]
digitosPi n = init [read [c] | c <- (x:xs)]
  where (x:_:xs) = showCReal n pi

--    λ> diferenciasConsecutivos (digitosPi 18)
--    [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]
diferenciasConsecutivos :: Num a => [a] -> [a]
diferenciasConsecutivos xs =
  zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]
distribucionDDCpi =
  distribucion . diferenciasConsecutivos . digitosPi
  where distribucion xs =
          elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()
graficas ns f = 
  plotLists [Key Nothing, PNG f]
            [puntos n | n <- ns]
  where puntos :: Int -> [(Int,Int)]
        puntos n = zip [-9..9] (distribucionDDCpi n)

import Data.Number.CReal

import Graphics.Gnuplot.Simple

import Data.Array

-- λ> digitosPi 18

-- [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]

digitosPi :: Int -> [Int]

digitosPi n = init [read [c] | c <- (x:xs)]

where (x:_:xs) = showCReal n pi

-- λ> diferenciasConsecutivos (digitosPi 18)

-- [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]

diferenciasConsecutivos :: Num a => [a] -> [a]

diferenciasConsecutivos xs =

zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]

distribucionDDCpi =

distribucion . diferenciasConsecutivos . digitosPi

where distribucion xs =

elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()

graficas ns f =

plotLists [Key Nothing, PNG f]

[puntos n | n <- ns]

where puntos :: Int -> [(Int,Int)]

puntos n = zip [-9..9] (distribucionDDCpi n)

Pensamiento

Doy consejo, a fuer de viejo:
nunca sigas mi consejo.

Antonio Machado

Huecos de Aquiles

PorJosé A. Alonso 2-abril-201926-marzo-2022

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

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

Definir las funciones

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

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

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]

1 2	λ> 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)

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)

Pensamiento

Tengo a mis amigos
en mi soledad;
cuando estoy con ellos
¡qué lejos están!

Antonio Machado

Cálculo de pi mediante la serie de Nilakantha

PorJosé A. Alonso 13-febrero-201925-marzo-2022

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

1 2	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

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]

1	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 |
+------+----------------+----------------+

+------+----------------+----------------+

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

1	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 |
+------+----------------+----------------+

+------+----------------+----------------+

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

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

Impares en filas del triángulo de Pascal

PorJosé A. Alonso 30-enero-20196-febrero-2019

El triángulo de Pascal es un triángulo de números

         1
        1 1
       1 2 1
     1  3 3  1
    1 4  6  4 1
   1 5 10 10 5 1
  ...............

1 1

1 2 1

1 3 3 1

1 4 6 4 1

1 5 10 10 5 1

...............

construido de la siguiente forma

la primera fila está formada por el número 1;
las filas siguientes se construyen sumando los números adyacentes de la fila superior y añadiendo un 1 al principio y al final de la fila.

Definir las funciones

   imparesPascal          :: [[Integer]]
   nImparesPascal         :: [Int]
   grafica_nImparesPascal :: Int -> IO ()

imparesPascal :: [[Integer]]

nImparesPascal :: [Int]

grafica_nImparesPascal :: Int -> IO ()

tales que

imparesPascal es la lista de los elementos impares en cada una de las filas del triángulo de Pascal. Por ejemplo,

     λ> take 8 imparesPascal
     [[1],[1,1],[1,1],[1,3,3,1],[1,1],[1,5,5,1],[1,15,15,1],[1,7,21,35,35,21,7,1]]

1 2	λ> take 8 imparesPascal [[1],[1,1],[1,1],[1,3,3,1],[1,1],[1,5,5,1],[1,15,15,1],[1,7,21,35,35,21,7,1]]

nImparesPascal es la lista del número de elementos impares en cada una de las filas del triángulo de Pascal. Por ejemplo,

     λ> take 32 nImparesPascal
     [1,2,2,4,2,4,4,8,2,4,4,8,4,8,8,16,2,4,4,8,4,8,8,16,4,8,8,16,8,16,16,32]
     λ> maximum (take (10^6) nImparesPascal3)
     524288

λ> take 32 nImparesPascal

[1,2,2,4,2,4,4,8,2,4,4,8,4,8,8,16,2,4,4,8,4,8,8,16,4,8,8,16,8,16,16,32]

λ> maximum (take (10^6) nImparesPascal3)

524288

(grafica_nImparesPascal n) dibuja la gráfica de los n primeros términos de nImparesPascal. Por ejemplo, (grafica_nImparesPascal 50) dibuja

y (grafica_nImparesPascal 100) dibuja

Comprobar con QuickCheck que todos los elementos de nImparesPascal son potencias de dos.