Polinomios de Bell

   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

   polBell :: Integer -> Polinomio Integer

   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

   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

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)