-- | Bernoulli and Euler polynomials
--
-- See <https://en.wikipedia.org/wiki/Bernoulli_polynomials>
-- 

{-# 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

--------------------------------------------------------------------------------
-- * Bernoulli polynomials

-- | Bernoulli polynomials
bernoulliB :: (Field c, KnownSymbol v) => Int -> Univariate c v
bernoulliB :: Int -> Univariate c v
bernoulliB = Univariate Rational v -> Univariate c v
forall c (v :: Symbol).
(Field c, KnownSymbol v) =>
Univariate Rational v -> Univariate c v
fromQUni (Univariate Rational v -> Univariate c v)
-> (Int -> Univariate Rational v) -> Int -> Univariate c v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Univariate Rational "x" -> Univariate Rational v
forall c (var1 :: Symbol) (var2 :: Symbol).
Univariate c var1 -> Univariate c var2
renameUniVar (Univariate Rational "x" -> Univariate Rational v)
-> (Int -> Univariate Rational "x") -> Int -> Univariate Rational v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Univariate Rational "x"
rationalBernoulliB

-- | Euler polynomials (not to be confused with the related Eulerian polynomials!)
eulerE :: (Field c, KnownSymbol v) => Int -> Univariate c v
eulerE :: Int -> Univariate c v
eulerE = Univariate Rational v -> Univariate c v
forall c (v :: Symbol).
(Field c, KnownSymbol v) =>
Univariate Rational v -> Univariate c v
fromQUni (Univariate Rational v -> Univariate c v)
-> (Int -> Univariate Rational v) -> Int -> Univariate c v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Univariate Rational "x" -> Univariate Rational v
forall c (var1 :: Symbol) (var2 :: Symbol).
Univariate c var1 -> Univariate c var2
renameUniVar (Univariate Rational "x" -> Univariate Rational v)
-> (Int -> Univariate Rational "x") -> Int -> Univariate Rational v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Univariate Rational "x"
rationalEulerE

--------------------------------------------------------------------------------

rationalBernoulliB :: Int -> Univariate Rational "x"
rationalBernoulliB :: Int -> Univariate Rational "x"
rationalBernoulliB Int
n = FreeMod Rational (U "x") -> Univariate Rational "x"
forall coeff (var :: Symbol).
FreeMod coeff (U var) -> Univariate coeff var
Uni (FreeMod Rational (U "x") -> Univariate Rational "x")
-> FreeMod Rational (U "x") -> Univariate Rational "x"
forall a b. (a -> b) -> a -> b
$ [(U "x", Rational)] -> FreeMod Rational (U "x")
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList
  [ (Int -> U "x"
forall (var :: Symbol). Int -> U var
U Int
k , Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial Int
n Int
k) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Int -> Rational
forall a. Integral a => a -> Rational
bernoulliNumber (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k) ) | Int
k<-[Int
0..Int
n] ]

rationalEulerE :: Int -> Univariate Rational "x"
rationalEulerE :: Int -> Univariate Rational "x"
rationalEulerE Int
n = [Univariate Rational "x"] -> Univariate Rational "x"
forall p. AlmostPolynomial p => [p] -> p
sumP
  [ CoeffP (Univariate Rational "x")
-> Univariate Rational "x" -> Univariate Rational "x"
forall p. AlmostPolynomial p => CoeffP p -> p -> p
scaleP Rational
CoeffP (Univariate Rational "x")
coeff (Univariate Rational "x" -> Univariate Rational "x")
-> Univariate Rational "x" -> Univariate Rational "x"
forall a b. (a -> b) -> a -> b
$ Int -> Univariate Rational "x"
xMinusHalfPowN (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)
  | Int
k <- [Int
0,Int
2..Int
n]
  , let coeff :: Rational
coeff = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial Int
n Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
eulerNumber Int
k) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2Rational -> Int -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^Int
k 
  ]
  where
    eulerNumber :: Int -> Integer
eulerNumber Int
k = if Int -> Bool
forall a. Integral a => a -> Bool
even Int
k 
      then Int -> Integer
signedEulerNumber (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
k Int
2) 
      else Integer
0

xMinusHalfPowN :: Int -> Univariate Rational "x"
xMinusHalfPowN :: Int -> Univariate Rational "x"
xMinusHalfPowN = ((Int -> Univariate Rational "x")
 -> Int -> Univariate Rational "x")
-> Int -> Univariate Rational "x"
forall a. ((Int -> a) -> Int -> a) -> Int -> a
intCache (Int -> Univariate Rational "x") -> Int -> Univariate Rational "x"
forall (var :: Symbol) t.
(KnownSymbol var, Eq t, Num t) =>
(t -> Univariate Rational var) -> t -> Univariate Rational var
compute where
  x_minus_half :: Univariate Rational var
x_minus_half = FreeMod Rational (U var) -> Univariate Rational var
forall coeff (var :: Symbol).
FreeMod coeff (U var) -> Univariate coeff var
Uni (FreeMod Rational (U var) -> Univariate Rational var)
-> FreeMod Rational (U var) -> Univariate Rational var
forall a b. (a -> b) -> a -> b
$ [(U var, Rational)] -> FreeMod Rational (U var)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList [ (Int -> U var
forall (var :: Symbol). Int -> U var
U Int
1 , Rational
1) , (Int -> U var
forall (var :: Symbol). Int -> U var
U Int
0 , -Rational
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2) ]
  compute :: (t -> Univariate Rational var) -> t -> Univariate Rational var
compute t -> Univariate Rational var
recur t
n = case t
n of 
    t
0 -> Univariate Rational var
1
    t
n -> Univariate Rational var
forall (var :: Symbol). Univariate Rational var
x_minus_half Univariate Rational var
-> Univariate Rational var -> Univariate Rational var
forall a. Num a => a -> a -> a
* t -> Univariate Rational var
recur (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1)

--------------------------------------------------------------------------------
-- * Bernoulli numbers

-- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for
-- k>2 and /odd/. This function uses the formula involving Stirling numbers
-- of the second kind. Numerators: A027641, denominators: A027642.
bernoulliNumber :: Integral a => a -> Rational
bernoulliNumber :: a -> Rational
bernoulliNumber a
n 
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<  a
0    = [Char] -> Rational
forall a. HasCallStack => [Char] -> a
error [Char]
"bernoulli: n should be nonnegative"
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = Rational
1
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1    = -Rational
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2
  | Bool
otherwise = [Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a -> Rational
f a
k | a
k<-[a
1..a
n] ] 
  where
    f :: a -> Rational
f a
k = Integer -> Rational
forall a. Real a => a -> Rational
toRational (a -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
k) (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ a -> Integer
forall a. Integral a => a -> Integer
factorial a
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* a -> a -> Integer
forall a. Integral a => a -> a -> Integer
stirling2nd a
n a
k) 
        Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ a -> Rational
forall a. Real a => a -> Rational
toRational (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
1)

--------------------------------------------------------------------------------
-- * Euler numbers

{-
x, oneminusx, xoneminusx, halfoneminusx :: Univariate Rational "x"
x = variableP ()
oneminusx = 1 - x
xoneminusx = x * (1 - x)
halfoneminusx = scaleP (1/2) oneminusx
nx :: Int -> Univariate Rational "x"
nx n = monomP' (U 1) (fromIntegral n)

scaledEulerianPolynomial :: Int -> Univariate Rational "x"
scaledEulerianPolynomial = intCache compute where
  compute recur n = case n of 
    0 -> 1
    n -> (nx n + halfoneminusx) * recur (n-1) + xoneminusx * differentiateUni (recur (n-1))
-}

x, oneminusx, two_xoneminusx :: Univariate Integer "x"
x :: Univariate Integer "x"
x = VarP (Univariate Integer "x") -> Univariate Integer "x"
forall p. AlmostPolynomial p => VarP p -> p
variableP ()
oneminusx :: Univariate Integer "x"
oneminusx = Univariate Integer "x"
1 Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
- Univariate Integer "x"
x
two_xoneminusx :: Univariate Integer "x"
two_xoneminusx = Univariate Integer "x"
2 Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
* Univariate Integer "x"
x Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
* (Univariate Integer "x"
1 Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
- Univariate Integer "x"
x)
two_nx :: Int -> Univariate Integer "x"
two_nx :: Int -> Univariate Integer "x"
two_nx Int
n = MonomP (Univariate Integer "x")
-> CoeffP (Univariate Integer "x") -> Univariate Integer "x"
forall p. AlmostPolynomial p => MonomP p -> CoeffP p -> p
monomP' (Int -> U "x"
forall (var :: Symbol). Int -> U var
U Int
1) (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)

-- | Eulerian polynomials (row polynomials of A060187)
eulerianPolynomial :: Int -> Univariate Integer "x"
eulerianPolynomial :: Int -> Univariate Integer "x"
eulerianPolynomial = ((Int -> Univariate Integer "x") -> Int -> Univariate Integer "x")
-> Int -> Univariate Integer "x"
forall a. ((Int -> a) -> Int -> a) -> Int -> a
intCache (Int -> Univariate Integer "x") -> Int -> Univariate Integer "x"
compute where
  compute :: (Int -> Univariate Integer "x") -> Int -> Univariate Integer "x"
compute Int -> Univariate Integer "x"
recur Int
n = case Int
n of 
    Int
0 -> Univariate Integer "x"
1
    Int
n -> (Int -> Univariate Integer "x"
two_nx Int
n Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
+ Univariate Integer "x"
oneminusx) Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
* Int -> Univariate Integer "x"
recur (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
+ Univariate Integer "x"
two_xoneminusx Univariate Integer "x"
-> Univariate Integer "x" -> Univariate Integer "x"
forall a. Num a => a -> a -> a
* Univariate Integer "x" -> Univariate Integer "x"
forall c (var :: Symbol).
(Ring c, KnownSymbol var) =>
Univariate c var -> Univariate c var
differentiateUni (Int -> Univariate Integer "x"
recur (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))

-- | Signed Euler numbers (unsigned version: A000364)
--
-- See <https://en.wikipedia.org/wiki/Euler_number>
--
-- NOTE: we skip the zeros (every other index)
signedEulerNumber :: Int -> Integer
signedEulerNumber :: Int -> Integer
signedEulerNumber = ((Int -> Integer) -> Int -> Integer) -> Int -> Integer
forall a. ((Int -> a) -> Int -> a) -> Int -> a
intCache (Int -> Integer) -> Int -> Integer
forall p. p -> Int -> Integer
compute where
  compute :: p -> Int -> Integer
compute p
_ Int
n = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
div ((CoeffP (Univariate Integer "x") -> Integer)
-> (VarP (Univariate Integer "x") -> Integer)
-> Univariate Integer "x"
-> Integer
forall p d.
(Polynomial p, Num d) =>
(CoeffP p -> d) -> (VarP p -> d) -> p -> d
evalP CoeffP (Univariate Integer "x") -> Integer
forall a. a -> a
id (\VarP (Univariate Integer "x")
_ -> -Integer
1) (Univariate Integer "x" -> Integer)
-> Univariate Integer "x" -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Univariate Integer "x"
eulerianPolynomial (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)) (Integer
4Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
n)

-- | unsigned Euler numbers (A000364)
--
-- NOTE: we skip the zeros (every other index)
unsignedEulerNumber :: Int -> Integer
unsignedEulerNumber :: Int -> Integer
unsignedEulerNumber Int
n = Int -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd Int
n (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Integer
signedEulerNumber Int
n

--------------------------------------------------------------------------------