Menu Close

Mes: mayo 2022

Polinomios de Bell

José A. Alonso, ,

Los polinomios de Bell forman una sucesión de polinomios, definida como sigue:

  • B₀(x) = 1 (polinomio unidad)
  • Bₙ(x) = x·[Bₙ(x) + Bₙ'(x)] (donde Bₙ'(x) es la derivada de Bₙ(x))

Por ejemplo,

   B₀(x) = 1                     = 1
   B₁(x) = x·(1+0)               = x     
   B₂(x) = x·(x+1)               = x²+x         
   B₃(x) = x·(x²+x+2x+1)         = x³+3x²+x    
   B₄(x) = x·(x³+3x²+x+3x²+6x+1) = x⁴+6x³+7x²+x

Definir la función

   polBell :: Integer -> Polinomio Integer

tal que (polBell n) es el polinomio de Bell de grado n. Por ejemplo,

   polBell 4                    ==  x^4 + 6*x^3 + 7*x^2 + 1*x
   coeficiente 2 (polBell 4)    ==  7
   coeficiente 2 (polBell 30)   ==  536870911
   coeficiente 1 (polBell 1000) == 1
   length (show (coeficiente 9 (polBell 2000)))  ==  1903

Notas: Se usa la librería I1M.PolOperaciones que se encuentra aquí y se describe aquí. Además, en el último ejemplo se usa la función coeficiente tal que (coeficiente k p) es el coeficiente del término de grado k en el polinomio p definida por

   coeficiente :: Num a => Int -> Polinomio a -> a
   coeficiente k p | k == n                 = coefLider p
                   | k > grado (restoPol p) = 0
                   | otherwise              = coeficiente k (restoPol p)
                   where n = grado p

Soluciones

import Data.List          (genericIndex)
import I1M.PolOperaciones (Polinomio, coefLider, consPol, derivada,
                           grado, multPol, polCero, polUnidad, restoPol,
                           sumaPol) 
import Test.QuickCheck    (Positive (Positive), quickCheck)
 
-- Función auxiliar
-- ================
 
-- (coeficiente k p) es el coeficiente del término de grado k en el
-- polinomio p.
coeficiente :: Num a => Int -> Polinomio a -> a
coeficiente k p | k == n                 = coefLider p
                | k > grado (restoPol p) = 0
                | otherwise              = coeficiente k (restoPol p)
                where n = grado p
 
-- 1ª solución
-- ===========
 
polBell1 :: Integer -> Polinomio Integer
polBell1 0 = polUnidad
polBell1 n = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))
  where p = polBell1 (n-1)
 
-- 2ª solución
-- ===========
 
polBell2 :: Integer -> Polinomio Integer
polBell2 n = sucPolinomiosBell `genericIndex` n
 
sucPolinomiosBell :: [Polinomio Integer]
sucPolinomiosBell = iterate f polUnidad
  where f p = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_polBell :: Positive Integer -> Bool 
prop_polBell (Positive n) =
  polBell1 n == polBell2 n
 
-- La comprobación es
--    λ> quickCheck prop_polBell
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> length (show (coeficiente 9 (polBell1 2000)))
--    1903
--    (5.37 secs, 4,829,322,368 bytes)
--    λ> length (show (coeficiente 9 (polBell2 2000)))
--    1903
--    (4.03 secs, 4,825,094,064 bytes)

El código se encuentra en GitHub.

Leer más →

Ordenación de los racionales

José A. Alonso, ,

En este ejercicio, representamos las fracciones mediante pares de números de enteros.

Definir la función

   fraccionesOrd :: Integer -> [(Integer,Integer)]

tal que (fraccionesOrd n) es la lista con las fracciones propias positivas ordenadas, con denominador menor o igual que n. Por ejemplo,

   λ> fraccionesOrd 4
   [(1,4),(1,3),(1,2),(2,3),(3,4)]
   λ> fraccionesOrd 5
   [(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

Soluciones

import Data.List (sort, sortBy)
import Data.Ratio ((%), numerator, denominator)
import Test.QuickCheck
 
-- 1ª solución
-- ===========
 
fraccionesOrd1 :: Integer -> [(Integer,Integer)]
fraccionesOrd1 n = 
  [(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))
                              | y <- [2..n], 
                                x <- [1..y-1], 
                                gcd x y == 1]]
 
-- 2ª solución
-- ===========
 
fraccionesOrd2 :: Integer -> [(Integer,Integer)]
fraccionesOrd2 n = 
  map snd (sort [(fromIntegral x/fromIntegral y,(x,y))
                 | y <- [2..n], 
                   x <- [1..y-1], 
                   gcd x y == 1])
 
-- 3ª solución
-- ===========
 
fraccionesOrd3 :: Integer -> [(Integer,Integer)]
fraccionesOrd3 n = 
  sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]
  where comp (a,b) (c,d) = compare (a*d) (b*c)
 
-- 4ª solución
-- ===========
 
fraccionesOrd4 :: Integer -> [(Integer,Integer)]
fraccionesOrd4 n = 
  [(numerator x, denominator x) | x <- racionalesOrd4 n]
 
-- (racionalesOrd4 n) es la lista con los racionales ordenados, con
-- denominador menor o igual que n. Por ejemplo,
--    λ> racionalesOrd4 4
--    [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]
racionalesOrd4 :: Integer -> [Rational]
racionalesOrd4 n =
  sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_fraccionesOrd :: Positive Integer -> Bool 
prop_fraccionesOrd (Positive n) =
  all (== fraccionesOrd1 n)
      [fraccionesOrd2 n,
       fraccionesOrd3 n,
       fraccionesOrd4 n]
 
-- La comprobación es
--    λ> quickCheck prop_fraccionesOrd
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> length (fraccionesOrd1 2000)
--    1216587
--    (3.65 secs, 2,879,842,368 bytes)
--    λ> length (fraccionesOrd2 2000)
--    1216587
--    (3.36 secs, 2,870,109,640 bytes)
--    λ> length (fraccionesOrd3 2000)
--    1216587
--    (8.83 secs, 5,700,519,584 bytes)
--    λ> length (fraccionesOrd4 2000)
--    1216587
--    (4.12 secs, 5,181,904,336 bytes)

El código se encuentra en GitHub.

Leer más →

Polinomios cuadráticos generadores de primos

José A. Alonso, ,

En 1772, Euler publicó que el polinomio n² + n + 41 genera 40 números primos para todos los valores de n entre 0 y 39. Sin embargo, cuando n = 40, 40²+40+41 = 40(40+1)+41 es divisible por 41.

Usando ordenadores, se descubrió que el polinomio n² – 79n + 1601 genera 80 números primos para todos los valores de n entre 0 y 79.

Definir la función

   generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

tal que (generadoresMaximales n) es el par (m,xs) donde

  • xs es la lista de pares (x,y) tales que n²+xn+y es uno de polinomios que genera un número máximo de números primos consecutivos a partir de cero entre todos los polinomios de la forma n²+an+b, con |a| ≤ n y |b| ≤ n y
  • m es dicho número máximo.

Por ejemplo,

   generadoresMaximales    4  ==  ( 3,[(-2,3),(-1,3),(3,3)])
   generadoresMaximales    6  ==  ( 5,[(-1,5),(5,5)])
   generadoresMaximales   41  ==  (41,[(-1,41)])
   generadoresMaximales   50  ==  (43,[(-5,47)])
   generadoresMaximales  100  ==  (48,[(-15,97)])
   generadoresMaximales  200  ==  (53,[(-25,197)])
   generadoresMaximales 1650  ==  (80,[(-79,1601)])

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (primes, isPrime)
import I1M.PolOperaciones (valor, consPol, polCero)
import Test.QuickCheck (Positive (Positive), quickCheck)
 
-- 1ª solución
-- ===========
 
generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales1 n =
  (m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])
  where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]
 
-- (nPrimos a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos (-1) 41     ==  41
--    nPrimos (-79) 1601  ==  80
nPrimos :: Integer -> Integer -> Int
nPrimos a b =
  length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])
 
-- 2ª solución
-- ===========
 
-- Notas:
-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que
--    0²+a*0+b = b es primo.
-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para
--    n = 1.
 
generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales2 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    a <- [-n..n],
                                    isPrime(1+a+b)]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys
 
-- 3ª solución
-- ===========
 
generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales3 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    p <- takeWhile (<=1+2*n) primes,
                                    let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys
 
-- 4ª solución (con la librería de polinomios)
-- ===========================================
 
generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales4 n = (m,map snd zs)
  where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,
                                     p <- takeWhile (<=1+2*n) primes,
                                     let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys
 
-- (nPrimos2 a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos2 (-1) 41     ==  41
--    nPrimos2 (-79) 1601  ==  80
nPrimos2 :: Integer -> Integer -> Int
nPrimos2 a b =
  length (takeWhile isPrime [valor p n | n <- [0..]])
  where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_generadoresMaximales :: Positive Integer -> Bool
prop_generadoresMaximales (Positive n) =
  all (equivalentes (generadoresMaximales1 n'))
      [generadoresMaximales2 n',
       generadoresMaximales3 n',
       generadoresMaximales4 n']
  where n' = n+1
 
equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool
equivalentes (n,xs) (m,ys) =
  n == m && sort xs == sort ys
 
-- La comprobación es
--    λ> quickCheck prop_generadoresMaximales
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> generadoresMaximales1 300
--    (56,[(-31,281)])
--    (2.10 secs, 2,744,382,760 bytes)
--    λ> generadoresMaximales2 300
--    (56,[(-31,281)])
--    (0.17 secs, 382,103,656 bytes)
--    λ> generadoresMaximales3 300
--    (56,[(-31,281)])
--    (0.19 secs, 346,725,872 bytes)
--    λ> generadoresMaximales4 300
--    (56,[(-31,281)])
--    (0.20 secs, 388,509,808 bytes)

El código se encuentra en GitHub.

Leer más →

El triángulo de Floyd

José A. Alonso, ,

El triángulo de Floyd, llamado así en honor a Robert Floyd, es un triángulo rectángulo formado con números naturales. Para crear un triángulo de Floyd, se comienza con un 1 en la esquina superior izquierda, y se continúa escribiendo la secuencia de los números naturales de manera que cada línea contenga un número más que la anterior. Las 5 primeras líneas del triángulo de Floyd son

    1
    2   3
    4   5   6
    7   8   9  10
   11  12  13  14  15

Definir la función

   trianguloFloyd :: [[Integer]]

tal que trianguloFloyd es el triángulo de Floyd. Por ejemplo,

   λ> take 4 trianguloFloyd
   [[1],
    [2,3],
    [4,5,6],
    [7,8,9,10]]
  (trianguloFloyd !! (10^5)) !! 0  ==  5000050001
  (trianguloFloyd !! (10^6)) !! 0  ==  500000500001
  (trianguloFloyd !! (10^7)) !! 0  ==  50000005000001

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)
 
-- 1ª solución
-- ===========
 
trianguloFloyd1 :: [[Integer]]
trianguloFloyd1 = floyd 1 [1..]
  where floyd n xs = i : floyd (n+1) r
          where (i,r) = splitAt n xs
 
-- 2ª solución
-- ===========
 
trianguloFloyd2 :: [[Integer]]
trianguloFloyd2 = iterate siguienteF [1]
 
-- (siguienteF xs) es la lista de los elementos de la línea xs en el
-- triángulo de Floyd. Por ejemplo,
--    siguienteF [2,3]    ==  [4,5,6]
--    siguienteF [4,5,6]  ==  [7,8,9,10]
siguienteF :: [Integer] -> [Integer]
siguienteF xs = [a..a+n]
    where a = 1 + last xs
          n = genericLength xs
 
-- 3ª solución
-- ===========
 
trianguloFloyd3 :: [[Integer]]
trianguloFloyd3 =
  [[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]
 
-- 4ª solución
-- ===========
 
trianguloFloyd4 :: [[Integer]]
trianguloFloyd4 =
  scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_trianguloFloyd :: Positive Int -> Bool
prop_trianguloFloyd (Positive n) =
  all (== (trianguloFloyd1 !! n))
      [trianguloFloyd2 !! n,
       trianguloFloyd3 !! n,
       trianguloFloyd4 !! n]
 
-- La comprobación es
-- λ> quickCheck prop_trianguloFloyd
-- +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> (trianguloFloyd1 !! 5000) !! 5000
--    12507501
--    (1.47 secs, 2,505,005,752 bytes)
--    λ> (trianguloFloyd2 !! 5000) !! 5000
--    12507501
--    (0.79 secs, 2,416,259,176 bytes)
--    λ> (trianguloFloyd3 !! 5000) !! 5000
--    12507501
--    (0.00 secs, 1,809,152 bytes)
--    λ> (trianguloFloyd4 !! 5000) !! 5000
--    12507501
--    (0.01 secs, 3,517,896 bytes)
--
--    λ> (trianguloFloyd3 !! (10^7)) !! 0
--    50000005000001
--    (2.45 secs, 1,656,534,080 bytes)
--    λ> (trianguloFloyd4 !! (10^7)) !! 0
--    50000005000001
--    (10.86 secs, 5,302,760,752 bytes)

El código se encuentra en GitHub.

Leer más →

Numeración con base múltiple

José A. Alonso, ,

Sea (b(i) | i ≥ 1) una sucesión infinita de números enteros mayores que 1. Entonces todo entero x mayor que cero se puede escribir de forma única como

   x = x(0) + x(1)b(1) +x(2)b(1)b(2) + ... + x(n)b(1)b(2)...b(n)

donde cada x(i) satisface la condición 0 ≤ x(i) < b(i+1). Se dice que [x(n),x(n-1),…,x(2),x(1),x(0)] es la representación de x en la base (b(i)). Por ejemplo, la representación de 377 en la base (2, 6, 8, …) es [7,5,0,1] ya que

   377 = 1 + 0*2 + 5*2*4 + 7*2*4*6

y, además, 0 ≤ 1 < 2, 0 ≤ 0 < 4, 0 ≤ 5 < 6 y 0 ≤ 7 < 8.

Definir las funciones

   decimalAmultiple :: [Integer] -> Integer -> [Integer]
   multipleAdecimal :: [Integer] -> [Integer] -> Integer

tales que

  • (decimalAmultiple bs x) es la representación del número x en la base bs. Por ejemplo,
     decimalAmultiple [2,4..] 377                      ==  [7,5,0,1]
     decimalAmultiple [2,5..] 377                      ==  [4,5,3,1]
     decimalAmultiple [2^n | n < - [1..]] 2015          ==  [1,15,3,3,1]
     decimalAmultiple (repeat 10) 2015                 ==  [2,0,1,5]
  • (multipleAdecimal bs cs) es el número decimal cuya representación en la base bs es cs. Por ejemplo,
     multipleAdecimal [2,4..] [7,5,0,1]                ==  377
     multipleAdecimal [2,5..] [4,5,3,1]                ==  377
     multipleAdecimal [2^n | n < - [1..]] [1,15,3,3,1]  ==  2015
     multipleAdecimal (repeat 10) [2,0,1,5]            ==  2015

Comprobar con QuickCheck que se verifican las siguientes propiedades

  • Para cualquier base bs y cualquier entero positivo n,
     multipleAdecimal bs (decimalAmultiple bs x) == x
  • Para cualquier base bs y cualquier entero positivo n, el coefiente i-ésimo de la representación múltiple de n en la base bs es un entero no negativo menos que el i-ésimo elemento de bs.

Soluciones

import Test.QuickCheck
import Data.List (unfoldr)
 
-- 1ª solución de decimalAmultiple
-- ===============================
 
decimalAmultiple1 :: [Integer] -> Integer -> [Integer]
decimalAmultiple1 bs n = reverse (aux bs n)
  where aux _ 0      = []
        aux (d:ds) m = r : aux ds q
          where (q,r) = quotRem m d
 
-- 2ª solución de decimalAmultiple
-- ===============================
 
decimalAmultiple2 :: [Integer] -> Integer -> [Integer]
decimalAmultiple2 bs n = aux bs n []
  where aux _ 0  xs     = xs
        aux (d:ds) m xs = aux ds q (r:xs)
          where (q,r) = quotRem m d
 
-- 3ª solución de decimalAmultiple
-- ===============================
 
decimalAmultiple3 :: [Integer] -> Integer -> [Integer]
decimalAmultiple3 xs n = reverse (unfoldr f (xs,n))
  where f (_     ,0) = Nothing
        f ((y:ys),m) = Just (r,(ys,q))
                       where (q,r) = quotRem m y
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_decimalAmultiple :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_decimalAmultiple (InfiniteList xs _) (Positive n) =
  all (== decimalAmultiple1 xs' n)
      [decimalAmultiple2 xs' n,
       decimalAmultiple3 xs' n]
  where xs' = map getPositive xs
 
-- Comparación de eficiencia de decimalAmultiple
-- =============================================
 
-- La comparación es
--    λ> length (decimalAmultiple1 [2,7..] (10^(10^5)))
--    21731
--    (0.45 secs, 486,085,256 bytes)
--    λ> length (decimalAmultiple2 [2,7..] (10^(10^5)))
--    21731
--    (0.32 secs, 485,563,664 bytes)
--    λ> length (decimalAmultiple3 [2,7..] (10^(10^5)))
--    21731
--    (0.44 secs, 487,649,768 bytes)
 
-- 1ª solución de multipleAdecimal
-- ===============================
 
multipleAdecimal1  :: [Integer] -> [Integer] -> Integer
multipleAdecimal1 xs ns = aux xs (reverse ns)
  where aux (y:ys) (m:ms) = m + y * (aux ys ms)
        aux _ _           = 0
 
-- 2ª solución de multipleAdecimal
-- ===============================
 
multipleAdecimal2 :: [Integer] -> [Integer] -> Integer
multipleAdecimal2 bs xs =
  sum (zipWith (*) (reverse xs) (1 : scanl1 (*) bs))
 
-- Comprobación de equivalencia de multipleAdecimal
-- ================================================
 
-- La propiedad es
prop_multipleAdecimal :: InfiniteList (Positive Integer) -> [Positive Integer] -> Bool
prop_multipleAdecimal (InfiniteList xs _) ys =
  multipleAdecimal1 xs' ys' == multipleAdecimal2 xs' ys'
  where xs' = map getPositive xs
        ys' = map getPositive ys
 
-- La comprobación es
--    λ> quickCheck prop_multipleAdecimal
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia de multipleAdecimal
-- =============================================
 
-- La comparación es
--    λ> length (show (multipleAdecimal1 [2,3..] [1..10^4]))
--    35660
--    (0.14 secs, 179,522,152 bytes)
--    λ> length (show (multipleAdecimal2 [2,3..] [1..10^4]))
--    35660
--    (0.22 secs, 243,368,664 bytes)
 
-- Comprobación de las propiedades
-- ===============================
 
-- La primera propiedad es
prop_inversas :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_inversas (InfiniteList xs _) (Positive n) =
  multipleAdecimal1 xs' (decimalAmultiple1 xs' n) == n
  where xs' = map getPositive xs
 
-- Su comprobación es
--    λ> quickCheck prop_inversas
--    +++ OK, passed 100 tests.
 
-- la 2ª propiedad es
prop_coeficientes :: InfiniteList (Positive Integer) -> Positive Integer -> Bool
prop_coeficientes (InfiniteList xs _) (Positive n) =
  and [0 <= c && c < b | (c,b) <- zip cs xs']
  where xs' = map getPositive xs
        cs = reverse (decimalAmultiple1 xs' n)
 
-- Su comprobación es
--    λ> quickCheck prop_coeficientes
--    +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Leer más →