Menu Close

Etiqueta: product

La sucesión de Sylvester

José A. Alonso---6 Comentarios

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]]
Seguir leyendo →

Cálculo de pi mediante la variante de Euler de la serie armónica

José A. Alonso---3 Comentarios

En el artículo El desarrollo más bello de Pi como suma infinita, Miguel Ángel Morales comenta el desarrollo de pi publicado por Leonhard Euler en su libro “Introductio in Analysis Infinitorum” (1748).

El desarrollo es el siguiente
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_1
y se obtiene a partir de la serie armónica
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_2
modificando sólo el signo de algunos términos según el siguiente criterio:

  • Dejamos un + cuando el denominador de la fracción sea un 2 o un primo de la forma 4m-1.
  • Cambiamos a – si el denominador de la fracción es un primo de la forma 4m+1.
  • Si el número es compuesto ponemos el signo que quede al multiplicar los signos correspondientes a cada factor.

Por ejemplo,

  • la de denominador 3 = 4×1-1 lleva un +,
  • la de denominador 5 = 4×1+1 lleva un -,
  • la de denominador 13 = 4×3+1 lleva un -,
  • la de denominador 6 = 2×3 lleva un + (porque los dos llevan un +),
  • la de denominador 10 = 2×5 lleva un – (porque el 2 lleva un + y el 5 lleva un -) y
  • la de denominador 50 = 5x5x2 lleva un + (un – por el primer 5, otro – por el segundo 5 y un + por el 2).

Definir las funciones

  aproximacionPi :: Int -> Double
  grafica        :: Int -> IO ()

tales que

  • (aproximacionPi n) es la aproximación de pi obtenida sumando los n primeros términos de la serie de Euler. Por ejemplo.
     aproximacionPi 1        ==  1.0
     aproximacionPi 10       ==  2.3289682539682537
     aproximacionPi 100      ==  2.934318000847734
     aproximacionPi 1000     ==  3.0603246224585128
     aproximacionPi 10000    ==  3.1105295744825403
     aproximacionPi 100000   ==  3.134308801935256
     aproximacionPi 1000000  ==  3.1395057903490806
  • (grafica n) dibuja la gráfica de las aproximaciones de pi usando k sumando donde k toma los valores de la lista [100,110..n]. Por ejemplo, al evaluar (grafica 4000) se obtiene
    Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_3.png

Soluciones

import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
 
-- 1ª definición
-- =============
 
aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [1 / fromIntegral (k * signo k) | k <- [1..n]] 
 
signoPrimo :: Int -> Int
signoPrimo 2 = 1
signoPrimo p | p `mod` 4 == 3 = 1
             | otherwise      = -1
 
signo :: Int -> Int
signo n | isPrime n = signoPrimo n
        | otherwise = product (map signoPrimo (primeFactors n))
 
-- 2ª definición
-- =============
 
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = serieEuler !! (n-1)
 
serieEuler :: [Double]
serieEuler =
  scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]
 
-- Definición de grafica
-- =====================
 
grafica :: Int -> IO ()
grafica n = 
    plotList [Key Nothing]
             [(k,aproximacionPi2 k) | k <- [100,110..n]]

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>
Seguir leyendo →

Productos de dos y tres números consecutivos

José A. Alonso---6 Comentarios

Definir la función

   productos :: Integer -> Integer -> [[Integer]]

tal que (productos n x) es las listas de n elementos consecutivos cuyo producto es x. Por ejemplo,

   productos 2 6     ==  [[2,3]]
   productos 3 6     ==  [[1,2,3]]
   productos 4 1680  ==  [[5,6,7,8]]
   productos 2 5     ==  []

Comprobar con QuickCheck que si n > 0 y x > 0, entonces

   productos n (product [x..x+n-1]) == [[x..x+n-1]]

Usando productos, definir la función

   productosDe2y3consecutivos :: [Integer]

cuyos elementos son los números naturales (no nulos) que pueden expresarse simultáneamente como producto de dos y tres números consecutivos. Por ejemplo,

   head productosDe2y3consecutivos  ==  6

Nota. Según demostró Mordell en 1962, productosDe2y3consecutivos sólo tiene dos elementos.

Soluciones

import Test.QuickCheck
 
-- 1ª definición
productos1 :: Integer -> Integer -> [[Integer]]
productos1 n x =
  [[y..y+n-1] | y <- [1..x]
              , product [y..y+n-1] == x]
 
-- 2ª definición
productos2 :: Integer -> Integer -> [[Integer]]
productos2 n x =
  [[z..z+n-1] | z <- [1..y]
              , product [z..z+n-1] == x]
  where y = head (filter (\y -> y^n >= x) [2..])
 
productos :: Integer -> Integer -> [[Integer]]
productos = productos2
 
prop_productos n x =
  n > 0 && x > 0 ==> productos n (product [x..x+n-1]) == [[x..x+n-1]]
 
-- La comprobación es
--    λ> quickCheck prop_productos
--    +++ OK, passed 100 tests.
--    (0.10 secs, 26409644 bytes)
 
productosDe2y3consecutivos :: [Integer]
productosDe2y3consecutivos =
  [x | x <- [1..] 
     , (not . null) (productos 2 x)
     , (not . null) (productos 3 x)]
 
-- El cálculo es
--    λ> take 2 productosDe2y3consecutivos
--    [6,210]
--    λ> productos 2 210
--    [[14,15]]
--    λ> productos 3 210
--    [[5,6,7]]

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

“El verdadero viaje de descubrimiento no consiste en buscar nuevos paisajes sino en tener nuevos ojos.”

Marcel Proust.

Seguir leyendo →

Números de Bell

José A. Alonso---1 Comentario

Una partición de un conjunto A es un conjunto de subconjuntos no vacíos de A, disjuntos dos a dos y cuya unión es A. Por ejemplo, el conjunto {1, 2, 3} tiene exactamente 5 particiones:

   {{1}, {2}, {3}}
   {{1,2}, {3}}
   {{1,3}, {2}}
   {{1}, {2,3}}
   {{1,2,3}}

El n-ésimo número de Bell, B(n), es el número de particiones de un conjunto de n elementos. Por lo visto anteriormentem B(3) = 5.

Definir las funciones

   particiones :: [a] -> [[[a]]]
   bell :: Integer -> Integer

tales que

  • (particiones xs) es el conjunto de las particiones de xs. Por ejemplo,
     λ> particiones [1,2]
     [[[1,2]],[[1],[2]]]
     λ> particiones [1,2,3]
     [[[1,2,3]],[[1],[2,3]],[[1,2],[3]],[[2],[1,3]],[[1],[2],[3]]]
     λ> particiones "abcd"
     [["abcd"],["a","bcd"],["ab","cd"],["b","acd"],["abc","d"],["bc","ad"],
      ["ac","bd"],["c","abd"],["a","b","cd"],["a","bc","d"],["a","c","bd"],
      ["ab","c","d"],["b","ac","d"],["b","c","ad"],["a","b","c","d"]]
  • (bell n) es el n-ésimo número de Bell. Por ejemplo,
     λ> bell 3
     5
     λ> map bell [0..10]
     [1,1,2,5,15,52,203,877,4140,21147,115975]

Comprobar con QuickCheck que (bell n) es equivalente a la función B(n) definida por

  • B(0) = 1
  • B(n) = \displaystyle \sum_{k=0}^{n-1} \binom{n-1}{k} B(k)

Soluciones

import Data.List (genericLength)
import Test.QuickCheck
 
-- Definición de particiones
-- =========================
 
particiones :: [a] -> [[[a]]]
particiones [] = [[]]
particiones (x:xs) =
  concat [([x] : yss) : inserta x yss | yss <- ysss]
  where ysss = particiones xs
 
-- (inserta x yss) es la lista obtenida insertando x en cada uno de los
-- elementos de yss. Por ejemplo, 
--    λ> inserta 1 [[2,3],[4],[5,6,7]]
--    [[[1,2,3],[4],[5,6,7]],[[2,3],[1,4],[5,6,7]],[[2,3],[4],[1,5,6,7]]]
inserta :: a -> [[a]] -> [[[a]]]
inserta _ []       = []
inserta x (ys:yss) = ((x:ys):yss) : [ys : zs | zs <- inserta x yss] 
 
-- Definición de Bell
-- ==================
 
bell :: Integer -> Integer
bell n = genericLength (particiones [1..n])
 
-- Propiedad
-- =========
 
prop_Bell :: Integer -> Property
prop_Bell n =
  n >= 0 ==> bell n == b n
 
b :: Integer -> Integer
b 0 = 1
b n = sum [comb (n-1) k * b k | k <- [0..n-1]]
 
comb :: Integer -> Integer -> Integer
comb n k = product [n-k+1..n] `div` product [1..k]
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_Bell
--    +++ 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

“Cambiemos nuestra actitud tradicional en la construcción de programas. En lugar de imaginar que nuestra tarea principal es indicarle a una computadora lo que debe hacer, concentrémonos más bien en explicarle a los seres humanos lo que queremos que haga una computadora.”

Donald Knuth.

Seguir leyendo →

Acotación del primorial

José A. Alonso---3 Comentarios

El primorial de un número natural n es el producto de todos los números primos menores o iguales a n. Por ejemplo, el primorial de 5 es 30 porque el producto de los primos menores o iguales que 5 es

   2 * 3 * 5 = 30

La propiedad de Erdös de acotación de los primoriales afirma que

Para todo número natural n, su primorial es menor o igual que 4ⁿ.

Definir las funciones

   primorial :: Integer -> Integer
   primoriales :: [Integer]

tales que

  • (primorial n) es el primorial de n. Por ejemplo,
     primorial 3  ==  6
     primorial 5  ==  30
     primorial 8  ==  210
  • primoriales es la sucesión de los primoriales. Por ejemplo,
   λ> take 15 primoriales
   [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

Comprobar con QuickCheck la propiedad de Erdös de acotación de los primoriales.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck
 
-- 1ª definición de primorial
-- ==========================
 
primorial :: Integer -> Integer
primorial n = product (takeWhile (<= n) primes)
 
-- 2ª definición de primorial
-- ==========================
 
primorial2 :: Integer -> Integer
primorial2 0 = 1
primorial2 n | gcd n x == 1 = n*x
             | otherwise    = x
  where x = primorial2 (n-1)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> length (show (primorial (5*10^5)))
--    216852
--    (1.65 secs, 2,472,977,584 bytes)
--    λ> length (show (primorial2 (5*10^5)))
--    216852
--    (3.56 secs, 2,719,162,272 bytes)
 
-- 1ª definición de primoriales
-- ============================
 
--    λ> take 15 primoriales
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales :: [Integer]
primoriales = map primorial [0..]
 
-- 2ª definición de primoriales
-- ============================
 
--    λ> take 15 primoriales2
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales2 :: [Integer]
primoriales2 = map primorial2 [0..]
 
-- 3ª definición de primoriales
-- ============================
 
--    λ> take 15 primoriales3
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales3 :: [Integer]
primoriales3 = scanl1 f [1..]
  where f x n | gcd n x == 1 = n*x
              | otherwise    = x
 
-- Comparación de eficiencia
-- =========================
 
--    λ> minimum (take 5000 primoriales)
--    1
--    (1.56 secs, 4,857,760,464 bytes)
--    λ> minimum (take 5000 primoriales2)
--    1
--    (9.39 secs, 10,942,848,240 bytes)
--    λ> minimum (take 5000 primoriales3)
--    1
--    (0.01 secs, 5,575,024 bytes)
--    
--    λ> minimum (take 6000 primoriales)
--    1
--    (2.22 secs, 7,013,937,248 bytes)
--    λ> minimum (take 6000 primoriales3)
--    1
--    (0.01 secs, 6,737,328 bytes)
 
-- Propiedad
-- =========
 
prop_primorial :: Integer -> Property
prop_primorial n =
  n >= 0 ==> primorial n <= 4^n
 
-- La comprobación es
--    λ> quickCheck prop_primorial
--    +++ 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

“Las matemáticas son la reina de las ciencias y la teoría de los números es la reina de las matemáticas.”

Carl Friedrich Gauss.

Seguir leyendo →