{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}
module Math.Algebra.Polynomial.Univariate.Bernoulli
( bernoulliB, eulerE
, rationalBernoulliB, rationalEulerE
, bernoulliNumber, signedEulerNumber, unsignedEulerNumber
, eulerianPolynomial
)
where
import Data.List
import Data.Ratio
import Data.Semigroup
import Data.Monoid
import GHC.TypeLits
import qualified Math.Algebra.Polynomial.FreeModule as ZMod
import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )
import Math.Algebra.Polynomial.Univariate
import Math.Algebra.Polynomial.Class
import Math.Algebra.Polynomial.Pretty
import Math.Algebra.Polynomial.Misc
bernoulliB :: (Field c, KnownSymbol v) => Int -> Univariate c v
bernoulliB = fromQUni . renameUniVar . rationalBernoulliB
eulerE :: (Field c, KnownSymbol v) => Int -> Univariate c v
eulerE = fromQUni . renameUniVar . rationalEulerE
rationalBernoulliB :: Int -> Univariate Rational "x"
rationalBernoulliB n = Uni $ ZMod.fromList
[ (U k , fromInteger (binomial n k) * bernoulliNumber (n-k) ) | k<-[0..n] ]
rationalEulerE :: Int -> Univariate Rational "x"
rationalEulerE n = sumP
[ scaleP coeff $ xMinusHalfPowN (n-k)
| k <- [0,2..n]
, let coeff = fromInteger (binomial n k * eulerNumber k) / 2^k
]
where
eulerNumber k = if even k
then signedEulerNumber (div k 2)
else 0
xMinusHalfPowN :: Int -> Univariate Rational "x"
xMinusHalfPowN = intCache compute where
x_minus_half = Uni $ ZMod.fromList [ (U 1 , 1) , (U 0 , -1/2) ]
compute recur n = case n of
0 -> 1
n -> x_minus_half * recur (n-1)
bernoulliNumber :: Integral a => a -> Rational
bernoulliNumber n
| n < 0 = error "bernoulli: n should be nonnegative"
| n == 0 = 1
| n == 1 = -1/2
| otherwise = sum [ f k | k<-[1..n] ]
where
f k = toRational (negateIfOdd (n+k) $ factorial k * stirling2nd n k)
/ toRational (k+1)
x, oneminusx, two_xoneminusx :: Univariate Integer "x"
x = variableP ()
oneminusx = 1 - x
two_xoneminusx = 2 * x * (1 - x)
two_nx :: Int -> Univariate Integer "x"
two_nx n = monomP' (U 1) (2 * fromIntegral n)
eulerianPolynomial :: Int -> Univariate Integer "x"
eulerianPolynomial = intCache compute where
compute recur n = case n of
0 -> 1
n -> (two_nx n + oneminusx) * recur (n-1) + two_xoneminusx * differentiateUni (recur (n-1))
signedEulerNumber :: Int -> Integer
signedEulerNumber = intCache compute where
compute _ n = div (evalP id (\_ -> -1) $ eulerianPolynomial (2*n)) (4^n)
unsignedEulerNumber :: Int -> Integer
unsignedEulerNumber n = negateIfOdd n $ signedEulerNumber n