Menu Close

Etiqueta: Lista infinita

Operaciones con series de potencias

Una serie de potencias es una serie de la forma

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

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

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

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

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

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

   type Serie a = [a]

Definir las siguientes funciones

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

tales que

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

Soluciones

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

Cálculo de pi mediante la serie de Nilakantha

Una serie infinita para el cálculo de pi, publicada por Nilakantha en el siglo XV, es

Definir las funciones

   aproximacionPi :: Int -> Double
   tabla          :: FilePath -> [Int] -> IO ()

tales que

  • (aproximacionPi n) es la n-ésima aproximación de pi obtenido sumando los n primeros términos de la serie de Nilakantha. Por ejemplo,
     aproximacionPi 0        ==  3.0
     aproximacionPi 1        ==  3.1666666666666665
     aproximacionPi 2        ==  3.1333333333333333
     aproximacionPi 3        ==  3.145238095238095
     aproximacionPi 4        ==  3.1396825396825396
     aproximacionPi 5        ==  3.1427128427128426
     aproximacionPi 10       ==  3.1414067184965018
     aproximacionPi 100      ==  3.1415924109719824
     aproximacionPi 1000     ==  3.141592653340544
     aproximacionPi 10000    ==  3.141592653589538
     aproximacionPi 100000   ==  3.1415926535897865
     aproximacionPi 1000000  ==  3.141592653589787
     pi                      ==  3.141592653589793
  • (tabla f ns) escribe en el fichero f las n-ésimas aproximaciones de pi, donde n toma los valores de la lista ns, junto con sus errores. Por ejemplo, al evaluar la expresión
     tabla "AproximacionesPi.txt" [0,10..100]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|   10 | 3.141406718497 | 0.000185935093 |
|   20 | 3.141565734659 | 0.000026918931 |
|   30 | 3.141584272675 | 0.000008380915 |
|   40 | 3.141589028941 | 0.000003624649 |
|   50 | 3.141590769850 | 0.000001883740 |
|   60 | 3.141591552546 | 0.000001101044 |
|   70 | 3.141591955265 | 0.000000698325 |
|   80 | 3.141592183260 | 0.000000470330 |
|   90 | 3.141592321886 | 0.000000331704 |
|  100 | 3.141592410972 | 0.000000242618 |
+------+----------------+----------------+

al evaluar la expresión

     tabla "AproximacionesPi.txt" [0,500..5000]

hace que el contenido del fichero “AproximacionesPi.txt” sea

+------+----------------+----------------+
| n    | Aproximación   | Error          |
+------+----------------+----------------+
|    0 | 3.000000000000 | 0.141592653590 |
|  500 | 3.141592651602 | 0.000000001988 |
| 1000 | 3.141592653341 | 0.000000000249 |
| 1500 | 3.141592653516 | 0.000000000074 |
| 2000 | 3.141592653559 | 0.000000000031 |
| 2500 | 3.141592653574 | 0.000000000016 |
| 3000 | 3.141592653581 | 0.000000000009 |
| 3500 | 3.141592653584 | 0.000000000006 |
| 4000 | 3.141592653586 | 0.000000000004 |
| 4500 | 3.141592653587 | 0.000000000003 |
| 5000 | 3.141592653588 | 0.000000000002 |
+------+----------------+----------------+

Soluciones

import Text.Printf
 
-- 1ª solución
-- ===========
 
aproximacionPi :: Int -> Double
aproximacionPi n = serieNilakantha !! n
 
serieNilakantha :: [Double]
serieNilakantha = scanl1 (+) terminosNilakantha
 
terminosNilakantha :: [Double]
terminosNilakantha = zipWith (/) numeradores denominadores
  where numeradores   = 3 : cycle [4,-4]
        denominadores = 1 : [n*(n+1)*(n+2) | n <- [2,4..]]
 
-- 2ª solución
-- ===========
 
aproximacionPi2 :: Int -> Double
aproximacionPi2 = aux 3 2 1
  where aux x _ _ 0 = x
        aux x y z m =
          aux (x+4/product[y..y+2]*z) (y+2) (negate z) (m-1)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> aproximacionPi (2*10^5)
--    3.141592653589787
--    (0.82 secs, 145,964,728 bytes)
--    λ> aproximacionPi2 (2*10^5)
--    3.141592653589787
--    (2.27 secs, 432,463,496 bytes)
--    λ> aproximacionPi (3*10^5)
--    3.141592653589787
--    (0.34 secs, 73,056,488 bytes)
--    λ> aproximacionPi2 (3*10^5)
--    3.141592653589787
--    (3.24 secs, 648,603,824 bytes)
 
-- Definicioń de tabla
-- ===================
 
tabla :: FilePath -> [Int] -> IO ()
tabla f ns = do
  writeFile f (tablaAux ns)
 
tablaAux :: [Int] -> String
tablaAux ns =
     linea
  ++ cabecera
  ++ linea
  ++ concat [printf "| %4d | %.12f | %.12f |\n" n a e
            | n <- ns
            , let a = aproximacionPi n
            , let e = abs (pi - a)]
  ++ linea
 
linea :: String
linea = "+------+----------------+----------------+\n"
 
cabecera :: String
cabecera = "| n    | Aproximación   | Error          |\n"

Sucesiones sin progresiones aritméticas de longitud 3

Tres números x, y, z está en progresión aritmética (PA) si existe un d tal que y = x+d y z = y+d. Por ejemplo, 1, 3, 5 están en PA ya que 3 = 1+2 y 5 = 3+2.

Se considera la sucesión donde cada uno de sus términos es el número natural tal que no está en PA con cualesquiera dos términos anteriores de la sucesión. Por ejemplo, si representamos por f(n) el n-ésimo término de la sucesión, entonces

f(0) = 0, que es el menor número natural;
f(1) = 1, que es el menor número natural que no está en la sucesión; 
f(2) = 3, ya que [0, 1, 2] están en PA y
                 [0, 1, 3] no están en PA;
f(3) = 4, ya que [0, 1, 4], [0, 3, 4] y [1, 3, 4] no están en PA;
f(4) = 9, ya que se descartan
          + el 5 porque [1, 3, 5] están en PA
          + el 6 porque [0, 3, 6] están en PA
          + el 7 porque [1, 4, 7] están en PA
          + el 8 porque [0, 4, 8] estan en PA
          y se acepta el 9 porque no están en PA niguna de [0,1,9],
          [0,3,9], [0,4,9], [1,3,9], [1,4,9], [3,4,9].

Definir la sucesión

   sucesionSin3enPA :: [Integer]

donde cada uno de sus términos es el menor número natural tal que no está en PA con cualesquiera dos términos anteriores de la sucesión. Por ejemplo,

   λ> take 20 sucesionSin3enPA
   [0,1,3,4,9,10,12,13,27,28,30,31,36,37,39,40,81,82,84,85]
   λ> sucesionSin3enPA !! 250
   3270

Soluciones

import Data.List ((\\), delete, tails)
 
-- 1ª solución
-- ===========
 
sucesionSin3enPA :: [Integer]
sucesionSin3enPA = aux []  
  where aux xs = x : aux (x:xs)
          where x = head (noEnPAConDos xs)
 
noEnPAConDos :: [Integer] -> [Integer]
noEnPAConDos xs = [z | z <- [0..]
                   , z `notElem` xs
                   , and [z - y /= y - x | x <- xs
                                         , y <- dropWhile (<= x) xs]]
 
-- 2ª solución
-- ===========
 
-- Si x, y, z están en progresión aritmética, entonces
--    y = (x + z) / 2
-- Por tanto,
--    z = 2 * y - x
--
-- Teniendo en cuenta la relación auterior se define (aux xs ys) tal que
-- xs es la lista de los términos actuales de la sucesión e ys es la
-- lista de los números que no están en PA con cualesquieras dos de xs.
 
sucesionSin3enPA2 :: [Integer]
sucesionSin3enPA2 = aux [] [0..] 
  where
    aux xs (y:ys) = y : aux (y:xs) (ys \\ [2 * y - x | x <- xs])
 
-- 3º solución
-- ===========
 
sucesionSin3enPA3 :: [Integer]
sucesionSin3enPA3 = 0 : aux [1]
  where aux (x:xs) = x : aux (xs ++ [3*x,3*x+1])
 
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> sucesionSin3enPA !! 70
--    741
--    (7.97 secs, 5,444,053,240 bytes)
--    λ> sucesionSin3enPA2 !! 70
--    741
--    (0.03 secs, 35,781,952 bytes)
--    λ> sucesionSin3enPA3 !! 70
--    741
--    (0.01 secs, 192,864 bytes)
--    
--    λ> sucesionSin3enPA2 !! 350
--    7410
--    (9.46 secs, 10,119,851,016 bytes)
--    λ> sucesionSin3enPA3 !! 350
--    7410
--    (0.01 secs, 1,931,296 bytes)

Otras soluciones

  • Se pueden escribir otras soluciones en los comentarios.
  • El código se debe escribir entre una línea con <pre lang=”haskell”> y otra con </pre>

Pensamiento

Quien se vive se pierde, Abel decía.
¡Oh, distancia, distancia!, que la estrella
que nadie toca, guía.
¿Quién navegó sin ella?

Antonio Machado

Múltiplos con ceros y unos

Se observa que todos los primeros números naturales tienen al menos un múltiplo no nulo que está formado solamente por ceros y unos. Por ejemplo, 1×10=10, 2×5=10, 3×37=111, 4×25=100, 5×2=10, 6×185=1110; 7×143=1001; 8X125=1000; 9×12345679=111111111.

Definir la función

   multiplosCon1y0 :: Integer -> [Integer]

tal que (multiplosCon1y0 n) es la lista de los múltiplos de n cuyos dígitos son 1 ó 0. Por ejemplo,

   take 4 (multiplosCon1y0 3)      ==  [111,1011,1101,1110]
   take 3 (multiplosCon1y0 23)     ==  [110101,1011011,1101010]
   head (multiplosCon1y0 1234658)  ==  110101101101000000110

Comprobar con QuickCheck que todo entero positivo tiene algún múltiplo cuyos dígitos son 1 ó 0.

Soluciones

import Test.QuickCheck
 
-- 1ª definición
-- =============
 
multiplosCon1y0 :: Integer -> [Integer]
multiplosCon1y0 n = [x | x <- multiplos n
                       , todos1y0 x]
 
-- (multiplos n) es la lista de los múltiplos de n. Por ejemplo, 
--    take 12 (multiplos 5)  ==  [5,10,15,20,25,30,35,40,45,50,55,60]
multiplos :: Integer -> [Integer]
multiplos n = [n,2*n..]
 
-- (todos1y0 n) se verifica si todos los dígitos de n son el 1 o el
-- 0. Por ejmplo,
--    todos1y0 1101110  ==  True
--    todos1y0 1102110  ==  False
todos1y0 :: Integer -> Bool
todos1y0 n = all (`elem` "01") (show n)
 
-- 2ª definición
-- =============
 
multiplosCon1y0b :: Integer -> [Integer] 
multiplosCon1y0b n = 
    [x | x <- numerosCon1y0
       , x `rem` n == 0] 
 
-- numerosCon1y0 es la lista de los números cuyos dígitos son 1 ó 0. Por
-- ejemplo,  
--    ghci> take 15 numerosCon1y0
--    [1,10,11,100,101,110,111,1000,1001,1010,1011,1100,1101,1110,1111]
numerosCon1y0 :: [Integer]
numerosCon1y0 = 1 : concat [[10*x,10*x+1] | x <- numerosCon1y0]
 
-- Comparación de eficiencia
-- =========================
 
--    λ> head (multiplosCon1y0 9)
--    111111111
--    (7.70 secs, 10,853,320,456 bytes)
--    λ> head (multiplosCon1y0b 9)
--    111111111
--    (0.01 secs, 167,992 bytes)
 
-- Comprobación de la propiedad
-- ============================
 
-- La propiedad es
prop_existe_multiplosCon1y0 :: Integer -> Property
prop_existe_multiplosCon1y0 n = 
    n > 0 ==> (not . null) (multiplosCon1y0b n)
 
-- La comprobación es
--    λ> quickCheck prop_existe_multiplosCon1y0
--    +++ OK, passed 100 tests.

Pensamiento

Huye del triste amor, amor pacato,
sin peligro, sin venda ni aventura,
que espera del amor prenda segura,
porque en amor locura es lo sensato.

Antonio Machado

Números triangulares

La sucesión de los números triangulares se obtiene sumando los números naturales.

   *     *      *        *         *   
        * *    * *      * *       * *  
              * * *    * * *     * * * 
                      * * * *   * * * *
                               * * * * * 
   1     3      6        10        15

Así, los 5 primeros números triangulares son

    1 = 1
    3 = 1+2
    6 = 1+2+3
   10 = 1+2+3+4
   15 = 1+2+3+4+5

Definir la función

   triangulares :: [Integer]

tal que triangulares es la lista de los números triangulares. Por ejemplo,

   take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
   maximum (take (5*10^6) triangulares4)  ==  12500002500000

Comprobar con QuickCheck que entre dos números triangulares consecutivos siempre hay un número primo.

Soluciones

import Test.QuickCheck (Property, (==>), quickCheck)
import Data.Numbers.Primes (primes)
 
-- 1ª solución
-- ===========
 
triangulares :: [Integer]
triangulares = [sum [1..n] | n <- [1..]]
 
-- 2ª solución
-- ===========
 
triangulares2 :: [Integer]
triangulares2 = map triangular [1..]
 
-- (triangular n) es el n-ésimo número triangular. Por ejemplo, 
--    triangular 5  ==  15
triangular :: Integer -> Integer
triangular 1 = 1
triangular n = n + triangular (n-1)
 
-- 3ª solución
-- ===========
 
triangulares3 :: [Integer]
triangulares3 = 1 : [x+y | (x,y) <- zip [2..] triangulares]
 
-- 4ª solución
-- ===========
 
triangulares4 :: [Integer]
triangulares4 = scanl1 (+) [1..]
 
-- 5ª solución
-- ===========
 
triangulares5 :: [Integer]
triangulares5 = [(n*(n+1)) `div` 2 | n <- [1..]]
 
-- Comparación de eficiencia
-- =========================
 
--    λ> maximum (take (10^4) triangulares)
--    50005000
--    (2.10 secs, 8,057,774,104 bytes)
--    λ> maximum (take (10^4) triangulares2)
--    50005000
--    (18.89 secs, 12,142,690,784 bytes)
--    λ> maximum (take (10^4) triangulares3)
--    50005000
--    (0.01 secs, 4,600,976 bytes)
--    λ> maximum (take (10^4) triangulares4)
--    50005000
--    (0.01 secs, 3,643,192 bytes)
--    λ> maximum (take (10^4) triangulares5)
--    50005000
--    (0.02 secs, 5,161,464 bytes)
--    
--    λ> maximum (take (3*10^4) triangulares3)
--    450015000
--    (26.06 secs, 72,546,027,136 bytes)
--    λ> maximum (take (3*10^4) triangulares4)
--    450015000
--    (0.02 secs, 10,711,600 bytes)
--    λ> maximum (take (3*10^4) triangulares5)
--    450015000
--    (0.03 secs, 15,272,320 bytes)
--    
--    λ> maximum (take (5*10^6) triangulares4)
--    12500002500000
--    (1.67 secs, 1,772,410,336 bytes)
--    λ> maximum (take (5*10^6) triangulares5)
--    12500002500000
--    (4.09 secs, 2,532,407,720 bytes)
 
-- La propiedad es
prop_triangulares :: Int -> Property
prop_triangulares n =
  n >= 0 ==> siguientePrimo x < y
  where (x:y:_) = drop n triangulares4
 
-- (siguientePrimo n) es el menor primo mayor o igual que n. Por
-- ejemplo, 
--    siguientePrimo 14  ==  17
--    siguientePrimo 17  ==  17
siguientePrimo :: Integer -> Integer
siguientePrimo n = head (dropWhile (< n) primes)
 
-- La comprobación es
--    λ> quickCheck prop_triangulares
--    +++ OK, passed 100 tests.

Pensamiento

Autores, la escena acaba
con un dogma de teatro:
En el principio era la máscara.

Antonio Machado

Sumas de dos primos

Definir la sucesión

   sumasDeDosPrimos :: [Integer]

cuyos elementos son los números que se pueden escribir como suma de dos números primos. Por ejemplo,

    λ> take 23 sumasDeDosPrimos
    [4,5,6,7,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,28,30,31]
    λ> sumasDeDosPrimos !! (5*10^4)
    83674

Soluciones

import Data.Numbers.Primes (isPrime, primes)
 
-- 1ª definición
-- =============
 
sumasDeDosPrimos1 :: [Integer]
sumasDeDosPrimos1 =
  [n | n <- [1..], not (null (sumaDeDosPrimos1 n))]
 
-- (sumasDeDosPrimos1 n) es la lista de pares de primos cuya suma es
-- n. Por ejemplo,
--    sumaDeDosPrimos  9  ==  [(2,7),(7,2)]
--    sumaDeDosPrimos 16  ==  [(3,13),(5,11),(11,5),(13,3)]
--    sumaDeDosPrimos 17  ==  []
sumaDeDosPrimos1 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos1 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (< n) primes
 
-- 2ª definición
-- =============
 
sumasDeDosPrimos2 :: [Integer]
sumasDeDosPrimos2 =
  [n | n <- [1..], not (null (sumaDeDosPrimos2 n))]
 
-- (sumasDeDosPrimos2 n) es la lista de pares (x,y) de primos cuya suma
-- es n y tales que x <= y. Por ejemplo,
--    sumaDeDosPrimos2  9  ==  [(2,7)]
--    sumaDeDosPrimos2 16  ==  [(3,13),(5,11)]
--    sumaDeDosPrimos2 17  ==  []
sumaDeDosPrimos2 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos2 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (<= (n `div` 2)) primes
 
-- 3ª definición
-- =============
 
sumasDeDosPrimos3 :: [Integer]
sumasDeDosPrimos3 = filter esSumaDeDosPrimos3 [4..]
 
-- (esSumaDeDosPrimos3 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos3  9  ==  True
--    esSumaDeDosPrimos3 16  ==  True
--    esSumaDeDosPrimos3 17  ==  False
esSumaDeDosPrimos3 :: Integer -> Bool
esSumaDeDosPrimos3 n
  | odd n     = isPrime (n-2)
  | otherwise = any isPrime [n-x | x <- takeWhile (<= (n `div` 2)) primes]
 
-- 4ª definición
-- =============
 
-- Usando la conjetura de Goldbach que dice que "Todo número par mayor
-- que 2 puede escribirse como suma de dos números primos" .
 
sumasDeDosPrimos4 :: [Integer]
sumasDeDosPrimos4 = filter esSumaDeDosPrimos4 [4..]
 
-- (esSumaDeDosPrimos4 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos4  9  ==  True
--    esSumaDeDosPrimos4 16  ==  True
--    esSumaDeDosPrimos4 17  ==  False
esSumaDeDosPrimos4 :: Integer -> Bool
esSumaDeDosPrimos4 n = even n || isPrime (n-2)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> sumasDeDosPrimos1 !! 3000
--    4731
--    (6.66 secs, 3,278,830,304 bytes)
--    λ> sumasDeDosPrimos2 !! 3000
--    4731
--    (3.69 secs, 1,873,984,088 bytes)
--    λ> sumasDeDosPrimos3 !! 3000
--    4731
--    (0.35 secs, 175,974,016 bytes)
--    λ> sumasDeDosPrimos4 !! 3000
--    4731
--    (0.07 secs, 18,396,432 bytes)
--    
--    λ> sumasDeDosPrimos3 !! 30000
--    49785
--    (6.65 secs, 3,785,736,416 bytes)
--    λ> sumasDeDosPrimos4 !! 30000
--    49785
--    (1.06 secs, 590,767,736 bytes)

Referencia

Pensamiento

“Es el deber de todos los profesores, y de los profesores de matemáticas en particular, exponer a sus alumnos a problemas mucho más que a hechos.” ~ Paul Halmos.

Sucesión de sumas de dos números abundantes

Un número n es abundante si la suma de los divisores propios de n es mayor que n. El primer número abundante es el 12 (cuyos divisores propios son 1, 2, 3, 4 y 6 cuya suma es 16). Por tanto, el menor número que es la suma de dos números abundantes es el 24.

Definir la sucesión

   sumasDeDosAbundantes :: [Integer]

cuyos elementos son los números que se pueden escribir como suma de dos números abundantes. Por ejemplo,

   take 10 sumasDeDosAbundantes  ==  [24,30,32,36,38,40,42,44,48,50]

Soluciones

sumasDeDosAbundantes :: [Integer]
sumasDeDosAbundantes = [n | n <- [1..], esSumaDeDosAbundantes n]
 
-- (esSumaDeDosAbundantes n) se verifica si n es suma de dos números
-- abundantes. Por ejemplo,
--    esSumaDeDosAbundantes 24           ==  True
--    any esSumaDeDosAbundantes [1..22]  ==  False
esSumaDeDosAbundantes :: Integer -> Bool
esSumaDeDosAbundantes n = (not . null) [x | x <- xs, n-x `elem` xs] 
  where xs = takeWhile (<n) abundantes
 
-- abundantes es la lista de los números abundantes. Por ejemplo,
--    take 10 abundantes  ==  [12,18,20,24,30,36,40,42,48,54]
abundantes :: [Integer]
abundantes = [n | n <- [2..], abundante n]
 
-- (abundante n) se verifica si n es abundante. Por ejemplo,
--    abundante 12  ==  True
--    abundante 11  ==  False
abundante :: Integer -> Bool
abundante n = sum (divisores n) > n
 
-- (divisores n) es la lista de los divisores propios de n. Por ejemplo,
--    divisores 12  ==  [1,2,3,4,6]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n `div` 2], n `mod` x == 0]

Pensamiento

¿Dices que nada se crea?
Alfarero, a tus cacharros.
Haz tu copa y no te importe
si no puedes hacer barro.

Antonio Machado

Numeración de las ternas de números naturales

Las ternas de números naturales se pueden ordenar como sigue

   (0,0,0), 
   (0,0,1),(0,1,0),(1,0,0),
   (0,0,2),(0,1,1),(0,2,0),(1,0,1),(1,1,0),(2,0,0),
   (0,0,3),(0,1,2),(0,2,1),(0,3,0),(1,0,2),(1,1,1),(1,2,0),(2,0,1),...
   ...

Definir la función

   posicion :: (Int,Int,Int) -> Int

tal que (posicion (x,y,z)) es la posición de la terna de números
naturales (x,y,z) en la ordenación anterior. Por ejemplo,

   posicion (0,1,0)  ==  2
   posicion (0,0,2)  ==  4
   posicion (0,1,1)  ==  5

Comprobar con QuickCheck que

  • la posición de (x,0,0) es x(x²+6x+11)/6
  • la posición de (0,y,0) es y(y²+3y+ 8)/6
  • la posición de (0,0,z) es z(z²+3z+ 2)/6
  • la posición de (x,x,x) es x(9x²+14x+7)/2

Soluciones

import Test.QuickCheck
import Data.List (elemIndex)
import Data.Maybe (fromJust)
 
-- 1ª solución
-- ===========
 
posicion :: (Int,Int,Int) -> Int
posicion (x,y,z) = length (takeWhile (/= (x,y,z)) ternas)
 
-- ternas es la lista ordenada de las ternas de números naturales. Por ejemplo,
--    ghci> take 9 ternas
--    [(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,0,2),(0,1,1),(0,2,0),(1,0,1),(1,1,0)]
ternas :: [(Int,Int,Int)]
ternas = [(x,y,n-x-y) | n <- [0..], x <- [0..n], y <- [0..n-x]] 
 
-- 2ª solución
-- ===========
 
posicion2 :: (Int,Int,Int) -> Int
posicion2 t = aux 0 ternas
  where aux n (t':ts) | t' == t   = n
                      | otherwise = aux (n+1) ts
 
-- 3ª solución
-- ===========
 
posicion3 :: (Int,Int,Int) -> Int
posicion3 t = 
  head [n | (n,t') <- zip [0..] ternas, t' == t]
 
-- 4ª solución
-- ===========
 
posicion4 :: (Int,Int,Int) -> Int
posicion4 t = fromJust (elemIndex t ternas)
 
-- 5ª solución
-- ===========
 
posicion5 :: (Int,Int,Int) -> Int
posicion5 = fromJust . (`elemIndex` ternas)
 
-- Equivalencia
-- ============
 
-- La propiedad es
prop_posicion_equiv :: NonNegative Int
                       -> NonNegative Int
                       -> NonNegative Int
                       -> Bool
prop_posicion_equiv (NonNegative x) (NonNegative y) (NonNegative z) =
  all (== posicion (x,y,z)) [f (x,y,z) | f <- [ posicion2
                                              , posicion3
                                              , posicion4
                                              , posicion5
                                              ]]
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_posicion_equiv
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
--   λ> posicion (200,0,0)
--   1373700
--   (2.48 secs, 368,657,528 bytes)
--   λ> posicion2 (200,0,0)
--   1373700
--   (3.96 secs, 397,973,912 bytes)
--   λ> posicion3 (200,0,0)
--   1373700
--   (1.47 secs, 285,831,056 bytes)
--   λ> posicion4 (200,0,0)
--   1373700
--   (0.13 secs, 102,320 bytes)
--   λ> posicion5 (200,0,0)
--   1373700
--   (0.14 secs, 106,376 bytes)
 
-- Propiedades
-- ===========
 
-- La 1ª propiedad es
prop_posicion1 :: NonNegative Int -> Bool
prop_posicion1 (NonNegative x) =
  posicion (x,0,0) == x * (x^2 + 6*x + 11) `div` 6
 
-- Su comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_posicion1
--    +++ OK, passed 100 tests.
 
-- La 2ª propiedad es
prop_posicion2 :: NonNegative Int -> Bool
prop_posicion2 (NonNegative y) =
  posicion (0,y,0) == y * (y^2 + 3*y + 8) `div` 6
 
-- Su comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_posicion2
--    +++ OK, passed 100 tests.
 
-- La 3ª propiedad es
prop_posicion3 :: NonNegative Int -> Bool
prop_posicion3 (NonNegative z) =
  posicion (0,0,z) == z * (z^2 + 3*z + 2) `div` 6
 
-- Su comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_posicion3
--    +++ OK, passed 100 tests.
 
-- La 4ª propiedad es
prop_posicion4 :: NonNegative Int -> Bool
prop_posicion4 (NonNegative x) =
  posicion (x,x,x) == x * (9 * x^2 + 14 * x + 7) `div` 2
 
-- Su comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=20}) prop_posicion4
--    +++ OK, passed 100 tests.

Pensamiento

¿Cabe una comunión cordial entre hombres, que nos permita cantar en coro, animados de un mismo sentir?

Antonio Machado

El 2019 es un número de la suerte

Un número de la suerte es un número natural que se genera por una criba, similar a la criba de Eratóstenes, como se indica a continuación:

Se comienza con la lista de los números enteros a partir de 1:

   1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25...

Se eliminan los números de dos en dos

   1,  3,  5,  7,  9,   11,   13,   15,   17,   19,   21,   23,   25...

Como el segundo número que ha quedado es 3, se eliminan los números restantes de tres en tres:

   1,  3,      7,  9,         13,   15,         19,   21,         25...

Como el tercer número que ha quedado es 7, se eliminan los números restantes de siete en siete:

   1,  3,      7,  9,         13,   15,               21,         25...

Este procedimiento se repite indefinidamente y los supervivientes son los números de la suerte:

   1,3,7,9,13,15,21,25,31,33,37,43,49,51,63,67,69,73,75,79

Definir las funciones

   numerosDeLaSuerte  :: [Int]
   esNumeroDeLaSuerte :: Int -> Bool

tales que

  • numerosDeLaSuerte es la sucesión de los números de la suerte. Por ejemplo,
     λ> take 20 numerosDeLaSuerte
     [1,3,7,9,13,15,21,25,31,33,37,43,49,51,63,67,69,73,75,79]
     λ> numerosDeLaSuerte !! 277
     2019
     λ> numerosDeLaSuerte !! 2000
     19309
  • (esNumeroDeLaSuerte n) que se verifica si n es un número de la suerte. Por ejemplo,
   esNumeroDeLaSuerte 15    ==  True
   esNumeroDeLaSuerte 16    ==  False
   esNumeroDeLaSuerte 2019  ==  True

Soluciones

-- 1ª definición de numerosDeLaSuerte 
numerosDeLaSuerte :: [Int]
numerosDeLaSuerte = criba 3 [1,3..]
  where
    criba i (n:s:xs) =
      n : criba (i + 1) (s : [x | (k, x) <- zip [i..] xs
                                , rem k s /= 0])
 
-- 2ª definición de numerosDeLaSuerte 
numerosDeLaSuerte2 :: [Int]
numerosDeLaSuerte2 =  1 : criba 2 [1, 3..]
  where criba k xs = z : criba (k + 1) (aux xs)
          where z = xs !! (k - 1 )
                aux ws = us ++ aux vs
                  where (us, _:vs) = splitAt (z - 1) ws 
 
-- Comparación de eficiencia
-- =========================
 
--    λ> numerosDeLaSuerte2 !! 200
--    1387
--    (9.25 secs, 2,863,983,232 bytes)
--    λ> numerosDeLaSuerte !! 200
--    1387
--    (0.06 secs, 10,263,880 bytes)
 
-- Definición de esNumeroDeLaSuerte
esNumeroDeLaSuerte :: Int -> Bool
esNumeroDeLaSuerte n =
  n == head (dropWhile (<n) numerosDeLaSuerte)

Pensamiento

Ya es sólo brocal el pozo;
púlpito será mañana;
pasado mañana, trono.

Antonio Machado

La sucesión del reloj astronómico de Praga

La cadena infinita “1234321234321234321…”, formada por la repetición de los dígitos 123432, tiene una propiedad (en la que se basa el funcionamiento del reloj astronómico de Praga: la cadena se puede partir en una sucesión de números, de forma que la suma de los dígitos de dichos números es la serie de los números naturales, como se observa a continuación:

    1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, ...
    1, 2, 3, 4,  5,   6,  7,    8,   9,   10,    11, ...

Definir la lista

   reloj :: [Integer]

cuyos elementos son los términos de la sucesión anterior. Por ejemplo,

   ghci> take 11 reloj
   [1,2,3,4,32,123,43,2123,432,1234,32123]

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)
 
reloj :: [Int]
reloj = aux [1..] (cycle "123432")
    where aux (n:ns) xs = read ys : aux ns zs
              where (ys,zs) = prefijoSuma n xs
 
-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo
-- suma es n y el resto de xs. Por ejemplo,
--    prefijoSuma 6 "12343"  ==  ("123","43")
prefijoSuma :: Int -> String -> (String,String)
prefijoSuma n xs = 
    head [(us,vs) | (us,vs) <- zip (inits xs) (tails xs)
                  , sumaD us == n]
 
-- (sumaD xs) es la suma de los dígitos de xs. Por ejemplo,
--    sumaD "123"  ==  6
sumaD :: String -> Int
sumaD = sum . map digitToInt