-- | Some basic univariate power series expansions.
-- This module is not re-exported by "Math.Combinat".
--
-- Note: the \"@convolveWithXXX@\" functions are much faster than the equivalent
-- @(XXX \`convolve\`)@!
-- 
-- TODO: better names for these functions.
--

{-# LANGUAGE CPP, BangPatterns, GeneralizedNewtypeDeriving #-}
module Math.Combinat.Numbers.Series where

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

import Data.List

import Math.Combinat.Sign
import Math.Combinat.Numbers
import Math.Combinat.Partitions.Integer
import Math.Combinat.Helper

--------------------------------------------------------------------------------
-- * Trivial series

-- | The series [1,0,0,0,0,...], which is the neutral element for the convolution.
{-# SPECIALIZE unitSeries :: [Integer] #-}
unitSeries :: Num a => [a]
unitSeries :: forall a. Num a => [a]
unitSeries = a
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0

-- | Constant zero series
zeroSeries :: Num a => [a]
zeroSeries :: forall a. Num a => [a]
zeroSeries = forall a. a -> [a]
repeat a
0

-- | Power series representing a constant function
constSeries :: Num a => a -> [a]
constSeries :: forall a. Num a => a -> [a]
constSeries a
x = a
x forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0

-- | The power series representation of the identity function @x@
idSeries :: Num a => [a]
idSeries :: forall a. Num a => [a]
idSeries = a
0 forall a. a -> [a] -> [a]
: a
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0

-- | The power series representation of @x^n@
powerTerm :: Num a => Int -> [a]
powerTerm :: forall a. Num a => Int -> [a]
powerTerm Int
n = forall a. Int -> a -> [a]
replicate Int
n a
0 forall a. [a] -> [a] -> [a]
++ (a
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0)

--------------------------------------------------------------------------------
-- * Basic operations on power series

addSeries :: Num a => [a] -> [a] -> [a]
addSeries :: forall a. Num a => [a] -> [a] -> [a]
addSeries [a]
xs [a]
ys = forall a b c. a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith a
0 a
0 forall a. Num a => a -> a -> a
(+) [a]
xs [a]
ys

sumSeries :: Num a => [[a]] -> [a]
sumSeries :: forall a. Num a => [[a]] -> [a]
sumSeries [] = [a
0]
sumSeries [[a]]
xs = forall a. (a -> a -> a) -> [a] -> a
foldl1' forall a. Num a => [a] -> [a] -> [a]
addSeries [[a]]
xs

subSeries :: Num a => [a] -> [a] -> [a]
subSeries :: forall a. Num a => [a] -> [a] -> [a]
subSeries [a]
xs [a]
ys = forall a b c. a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith a
0 a
0 (-) [a]
xs [a]
ys

negateSeries :: Num a => [a] -> [a]
negateSeries :: forall a. Num a => [a] -> [a]
negateSeries = forall a b. (a -> b) -> [a] -> [b]
map forall a. Num a => a -> a
negate

scaleSeries :: Num a => a -> [a] -> [a]
scaleSeries :: forall a. Num a => a -> [a] -> [a]
scaleSeries a
s = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
*a
s)

-- | A different implementation, taken from:
--
-- M. Douglas McIlroy: Power Series, Power Serious 
mulSeries :: Num a => [a] -> [a] -> [a]
mulSeries :: forall a. Num a => [a] -> [a] -> [a]
mulSeries [a]
xs [a]
ys = forall a. Num a => [a] -> [a] -> [a]
go ([a]
xs forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) ([a]
ys forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) where
  go :: [a] -> [a] -> [a]
go (a
f:[a]
fs) ggs :: [a]
ggs@(a
g:[a]
gs) = a
fforall a. Num a => a -> a -> a
*a
g forall a. a -> [a] -> [a]
: (forall a. Num a => a -> [a] -> [a]
scaleSeries a
f [a]
gs) forall a. Num a => [a] -> [a] -> [a]
`addSeries` [a] -> [a] -> [a]
go [a]
fs [a]
ggs

-- | Multiplication of power series. This implementation is a synonym for 'convolve'
mulSeriesNaive :: Num a => [a] -> [a] -> [a]
mulSeriesNaive :: forall a. Num a => [a] -> [a] -> [a]
mulSeriesNaive = forall a. Num a => [a] -> [a] -> [a]
convolve

productOfSeries :: Num a => [[a]] -> [a]
productOfSeries :: forall a. Num a => [[a]] -> [a]
productOfSeries = forall a. Num a => [[a]] -> [a]
convolveMany

--------------------------------------------------------------------------------
-- * Convolution (product)

-- | Convolution of series (that is, multiplication of power series). 
-- The result is always an infinite list. Warning: This is slow!
convolve :: Num a => [a] -> [a] -> [a]
convolve :: forall a. Num a => [a] -> [a] -> [a]
convolve [a]
xs1 [a]
ys1 = [a]
res where
  res :: [a]
res = [ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' forall a. Num a => a -> a -> a
(+) a
0 (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [a]
xs (forall a. [a] -> [a]
reverse (forall a. Int -> [a] -> [a]
take Int
n [a]
ys)))
        | Int
n<-[Int
1..] 
        ]
  xs :: [a]
xs = [a]
xs1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
  ys :: [a]
ys = [a]
ys1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0

-- | Convolution (= product) of many series. Still slow!
convolveMany :: Num a => [[a]] -> [a]
convolveMany :: forall a. Num a => [[a]] -> [a]
convolveMany []  = a
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0
convolveMany [[a]]
xss = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall a. Num a => [a] -> [a] -> [a]
convolve [[a]]
xss

--------------------------------------------------------------------------------
-- * Reciprocals of general power series

-- | Division of series.
--
-- Taken from: M. Douglas McIlroy: Power Series, Power Serious 
divSeries :: (Eq a, Fractional a) => [a] -> [a] -> [a]
divSeries :: forall a. (Eq a, Fractional a) => [a] -> [a] -> [a]
divSeries [a]
xs [a]
ys = forall a. (Eq a, Fractional a) => [a] -> [a] -> [a]
go ([a]
xs forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) ([a]
ys forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) where
  go :: [a] -> [a] -> [a]
go (a
0:[a]
fs)     (a
0:[a]
gs) = [a] -> [a] -> [a]
go [a]
fs [a]
gs
  go (a
f:[a]
fs) ggs :: [a]
ggs@(a
g:[a]
gs) = let q :: a
q = a
fforall a. Fractional a => a -> a -> a
/a
g in a
q forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
go ([a]
fs forall a. Num a => [a] -> [a] -> [a]
`subSeries` forall a. Num a => a -> [a] -> [a]
scaleSeries a
q [a]
gs) [a]
ggs

-- | Given a power series, we iteratively compute its multiplicative inverse
reciprocalSeries :: (Eq a, Fractional a) => [a] -> [a]
reciprocalSeries :: forall a. (Eq a, Fractional a) => [a] -> [a]
reciprocalSeries [a]
series = case [a]
series of
  [] -> forall a. HasCallStack => [Char] -> a
error [Char]
"reciprocalSeries: empty input series (const 0 function does not have an inverse)"
  (a
a:[a]
as) -> case a
a of
    a
0 -> forall a. HasCallStack => [Char] -> a
error [Char]
"reciprocalSeries: input series has constant term 0"
    a
_ -> forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/a
a) forall a b. (a -> b) -> a -> b
$ forall a. (Eq a, Num a) => [a] -> [a]
integralReciprocalSeries forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/a
a) [a]
series

-- | Given a power series starting with @1@, we can compute its multiplicative inverse
-- without divisions.
--
{-# SPECIALIZE integralReciprocalSeries :: [Int]     -> [Int]     #-}
{-# SPECIALIZE integralReciprocalSeries :: [Integer] -> [Integer] #-}
integralReciprocalSeries :: (Eq a, Num a) => [a] -> [a]
integralReciprocalSeries :: forall a. (Eq a, Num a) => [a] -> [a]
integralReciprocalSeries [a]
series = case [a]
series of 
  [] -> forall a. HasCallStack => [Char] -> a
error [Char]
"integralReciprocalSeries: empty input series (const 0 function does not have an inverse)"
  (a
a:[a]
as) -> case a
a of
    a
1 -> a
1 forall a. a -> [a] -> [a]
: [a] -> [a]
worker [a
1]
    a
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"integralReciprocalSeries: input series must start with 1"
  where
    worker :: [a] -> [a]
worker [a]
bs = let b' :: a
b' = - forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (forall a. [a] -> [a]
tail [a]
series) [a]
bs) 
                in  a
b' forall a. a -> [a] -> [a]
: [a] -> [a]
worker (a
b'forall a. a -> [a] -> [a]
:[a]
bs)

--------------------------------------------------------------------------------
-- * Composition of formal power series

-- | @g \`composeSeries\` f@ is the power series expansion of @g(f(x))@.
-- This is a synonym for @flip substitute@.
--
-- This implementation is taken from
--
-- M. Douglas McIlroy: Power Series, Power Serious 
composeSeries :: (Eq a, Num a) => [a] -> [a] -> [a]
composeSeries :: forall a. (Eq a, Num a) => [a] -> [a] -> [a]
composeSeries [a]
xs [a]
ys = forall a. (Eq a, Num a) => [a] -> [a] -> [a]
go ([a]
xs forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) ([a]
ys forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) where
  go :: [a] -> [a] -> [a]
go (a
f:[a]
fs) (a
0:[a]
gs) = a
f forall a. a -> [a] -> [a]
: forall a. Num a => [a] -> [a] -> [a]
mulSeries [a]
gs ([a] -> [a] -> [a]
go [a]
fs (a
0forall a. a -> [a] -> [a]
:[a]
gs))
  go (a
f:[a]
fs) (a
_:[a]
gs) = forall a. HasCallStack => [Char] -> a
error [Char]
"PowerSeries/composeSeries: we expect the the constant term of the inner series to be zero"

-- | @substitute f g@ is the power series corresponding to @g(f(x))@. 
-- Equivalently, this is the composition of univariate functions (in the \"wrong\" order).
--
-- Note: for this to be meaningful in general (not depending on convergence properties),
-- we need that the constant term of @f@ is zero.
substitute :: (Eq a, Num a) => [a] -> [a] -> [a]
substitute :: forall a. (Eq a, Num a) => [a] -> [a] -> [a]
substitute [a]
f [a]
g = forall a. (Eq a, Num a) => [a] -> [a] -> [a]
composeSeries [a]
g [a]
f

-- | Naive implementation of 'composeSeries' (via 'substituteNaive')
composeSeriesNaive :: (Eq a, Num a) => [a] -> [a] -> [a]
composeSeriesNaive :: forall a. (Eq a, Num a) => [a] -> [a] -> [a]
composeSeriesNaive [a]
g [a]
f = forall a. (Eq a, Num a) => [a] -> [a] -> [a]
substituteNaive [a]
f [a]
g

-- | Naive implementation of 'substitute'
substituteNaive :: (Eq a, Num a) => [a] -> [a] -> [a]
substituteNaive :: forall a. (Eq a, Num a) => [a] -> [a] -> [a]
substituteNaive [a]
as_ [a]
bs_ = 
  case forall a. [a] -> a
head [a]
as of
    a
0 -> [ Int -> a
f Int
n | Int
n<-[Int
0..] ]
    a
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"PowerSeries/substituteNaive: we expect the the constant term of the inner series to be zero"
  where
    as :: [a]
as = [a]
as_ forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
    bs :: [a]
bs = [a]
bs_ forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
    a :: Int -> a
a Int
i = [a]
as forall a. [a] -> Int -> a
!! Int
i
    b :: Int -> a
b Int
j = [a]
bs forall a. [a] -> Int -> a
!! Int
j
    f :: Int -> a
f Int
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
            [ Int -> a
b Int
m forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)forall a b. (Num a, Integral b) => a -> b -> a
^Int
j | (Int
i,Int
j)<-[(Int, Int)]
es ] forall a. Num a => a -> a -> a
* forall a. Num a => Integer -> a
fromInteger (forall a. Integral a => [a] -> Integer
multinomial (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd [(Int, Int)]
es))
            | Partition
p <- Int -> [Partition]
partitions Int
n 
            , let es :: [(Int, Int)]
es = Partition -> [(Int, Int)]
toExponentialForm Partition
p
            , let m :: Int
m  = Partition -> Int
partitionWidth    Partition
p
            ]

--------------------------------------------------------------------------------
-- * Lagrange inversions

-- | We expect the input series to match @(0:a1:_)@. with a1 nonzero The following is true for the result (at least with exact arithmetic):
--
-- > substitute f (lagrangeInversion f) == (0 : 1 : repeat 0)
-- > substitute (lagrangeInversion f) f == (0 : 1 : repeat 0)
--
-- This implementation is taken from:
--
-- M. Douglas McIlroy: Power Series, Power Serious 
lagrangeInversion :: (Eq a, Fractional a) => [a] -> [a]
lagrangeInversion :: forall a. (Eq a, Fractional a) => [a] -> [a]
lagrangeInversion [a]
xs = forall a. (Eq a, Fractional a) => [a] -> [a]
go ([a]
xs forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0) where
  go :: [a] -> [a]
go (a
0:[a]
fs) = [a]
rs where rs :: [a]
rs = a
0 forall a. a -> [a] -> [a]
: forall a. (Eq a, Fractional a) => [a] -> [a] -> [a]
divSeries forall a. Num a => [a]
unitSeries (forall a. (Eq a, Num a) => [a] -> [a] -> [a]
composeSeries [a]
fs [a]
rs)
  go (a
_:[a]
fs) = forall a. HasCallStack => [Char] -> a
error [Char]
"lagrangeInversion: the series should start with (0 + a1*x + a2*x^2 + ...) where a1 is non-zero"

-- | Coefficients of the Lagrange inversion
lagrangeCoeff :: Partition -> Integer
lagrangeCoeff :: Partition -> Integer
lagrangeCoeff Partition
p = forall a. Integral a => a -> a -> a
div Integer
numer Integer
denom where
  numer :: Integer
numer = (-Integer
1)forall a b. (Num a, Integral b) => a -> b -> a
^Int
m forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (Integral a, Num b) => a -> b
fromIntegral [Int
nforall a. Num a => a -> a -> a
+Int
1..Int
nforall a. Num a => a -> a -> a
+Int
m])
  denom :: Integer
denom = forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nforall a. Num a => a -> a -> a
+Int
1) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (forall a b. (a -> b) -> [a] -> [b]
map (forall a. Integral a => a -> Integer
factorial forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd) [(Int, Int)]
es)
  m :: Int
m  = Partition -> Int
partitionWidth    Partition
p
  n :: Int
n  = Partition -> Int
partitionWeight   Partition
p
  es :: [(Int, Int)]
es = Partition -> [(Int, Int)]
toExponentialForm Partition
p

-- | We expect the input series to match @(0:1:_)@. The following is true for the result (at least with exact arithmetic):
--
-- > substitute f (integralLagrangeInversion f) == (0 : 1 : repeat 0)
-- > substitute (integralLagrangeInversion f) f == (0 : 1 : repeat 0)
--
integralLagrangeInversionNaive :: (Eq a, Num a) => [a] -> [a]
integralLagrangeInversionNaive :: forall a. (Eq a, Num a) => [a] -> [a]
integralLagrangeInversionNaive [a]
series_ = 
  case [a]
series of
    (a
0:a
1:[a]
rest) -> a
0 forall a. a -> [a] -> [a]
: a
1 forall a. a -> [a] -> [a]
: [ Int -> a
f Int
n | Int
n<-[Int
1..] ]
    [a]
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"integralLagrangeInversionNaive: the series should start with (0 + x + a2*x^2 + ...)"
  where
    series :: [a]
series = [a]
series_ forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
    as :: [a]
as  = forall a. [a] -> [a]
tail [a]
series 
    a :: Int -> a
a Int
i = [a]
as forall a. [a] -> Int -> a
!! Int
i
    f :: Int -> a
f Int
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ forall a. Num a => Integer -> a
fromInteger (Partition -> Integer
lagrangeCoeff Partition
p) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)forall a b. (Num a, Integral b) => a -> b -> a
^Int
j | (Int
i,Int
j) <- Partition -> [(Int, Int)]
toExponentialForm Partition
p ]
              | Partition
p <- Int -> [Partition]
partitions Int
n
              ] 

-- | Naive implementation of 'lagrangeInversion'
lagrangeInversionNaive :: (Eq a, Fractional a) => [a] -> [a]
lagrangeInversionNaive :: forall a. (Eq a, Fractional a) => [a] -> [a]
lagrangeInversionNaive [a]
series_ = 
  case [a]
series of
    (a
0:a
a1:[a]
rest) -> if a
a1 forall a. Eq a => a -> a -> Bool
==a
0 
      then forall {a}. a
err 
      else a
0 forall a. a -> [a] -> [a]
: (a
1forall a. Fractional a => a -> a -> a
/a
a1) forall a. a -> [a] -> [a]
: [ Int -> a
f Int
n forall a. Fractional a => a -> a -> a
/ a
a1forall a b. (Num a, Integral b) => a -> b -> a
^(Int
nforall a. Num a => a -> a -> a
+Int
1) | Int
n<-[Int
1..] ]
    [a]
_ -> forall {a}. a
err
  where
    err :: a
err    = forall a. HasCallStack => [Char] -> a
error [Char]
"lagrangeInversionNaive: the series should start with (0 + a1*x + a2*x^2 + ...) where a1 is non-zero"
    series :: [a]
series = [a]
series_ forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
    a1 :: a
a1  = [a]
series forall a. [a] -> Int -> a
!! Int
1
    as :: [a]
as  = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/a
a1) (forall a. [a] -> [a]
tail [a]
series)
    a :: Int -> a
a Int
i = [a]
as forall a. [a] -> Int -> a
!! Int
i
    f :: Int -> a
f Int
n = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ forall a. Num a => Integer -> a
fromInteger (Partition -> Integer
lagrangeCoeff Partition
p) forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)forall a b. (Num a, Integral b) => a -> b -> a
^Int
j | (Int
i,Int
j) <- Partition -> [(Int, Int)]
toExponentialForm Partition
p ]
              | Partition
p <- Int -> [Partition]
partitions Int
n
              ] 


--------------------------------------------------------------------------------
-- * Differentiation and integration

differentiateSeries :: Num a => [a] -> [a]
differentiateSeries :: forall a. Num a => [a] -> [a]
differentiateSeries (a
y:[a]
ys) = forall {t} {a}. (Integral t, Num a) => t -> [a] -> [a]
go (Int
1::Int) [a]
ys where
  go :: t -> [a] -> [a]
go !t
n (a
x:[a]
xs) = forall a b. (Integral a, Num b) => a -> b
fromIntegral t
n forall a. Num a => a -> a -> a
* a
x forall a. a -> [a] -> [a]
: t -> [a] -> [a]
go (t
nforall a. Num a => a -> a -> a
+t
1) [a]
xs
  go t
_  []     = []

integrateSeries :: Fractional a => [a] -> [a]
integrateSeries :: forall a. Fractional a => [a] -> [a]
integrateSeries [a]
ys = a
0 forall a. a -> [a] -> [a]
: forall {a} {t}. (Fractional a, Integral t) => t -> [a] -> [a]
go (Int
1::Int) [a]
ys where
  go :: t -> [a] -> [a]
go !t
n (a
x:[a]
xs) = a
x forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral t
n) forall a. a -> [a] -> [a]
: t -> [a] -> [a]
go (t
nforall a. Num a => a -> a -> a
+t
1) [a]
xs
  go t
_  []     = []

--------------------------------------------------------------------------------
-- * Power series expansions of elementary functions

-- | Power series expansion of @exp(x)@
expSeries :: Fractional a => [a]
expSeries :: forall a. Fractional a => [a]
expSeries = forall {t}. Fractional t => t -> t -> [t]
go a
0 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
e forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
1) (t
e forall a. Fractional a => a -> a -> a
/ (t
iforall a. Num a => a -> a -> a
+t
1))

-- | Power series expansion of @cos(x)@
cosSeries :: Fractional a => [a]
cosSeries :: forall a. Fractional a => [a]
cosSeries = forall {t}. Fractional t => t -> t -> [t]
go a
0 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
e forall a. a -> [a] -> [a]
: t
0 forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
2) (-t
e forall a. Fractional a => a -> a -> a
/ ((t
iforall a. Num a => a -> a -> a
+t
1)forall a. Num a => a -> a -> a
*(t
iforall a. Num a => a -> a -> a
+t
2)))

-- | Power series expansion of @sin(x)@
sinSeries :: Fractional a => [a]
sinSeries :: forall a. Fractional a => [a]
sinSeries = forall {t}. Fractional t => t -> t -> [t]
go a
1 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
0 forall a. a -> [a] -> [a]
: t
e forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
2) (-t
e forall a. Fractional a => a -> a -> a
/ ((t
iforall a. Num a => a -> a -> a
+t
1)forall a. Num a => a -> a -> a
*(t
iforall a. Num a => a -> a -> a
+t
2)))

-- | Alternative implementation using differential equations.
--
-- Taken from: M. Douglas McIlroy: Power Series, Power Serious
cosSeries2, sinSeries2 :: Fractional a => [a]
cosSeries2 :: forall a. Fractional a => [a]
cosSeries2 = forall a. Num a => [a]
unitSeries forall a. Num a => [a] -> [a] -> [a]
`subSeries` forall a. Fractional a => [a] -> [a]
integrateSeries forall a. Fractional a => [a]
sinSeries2
sinSeries2 :: forall a. Fractional a => [a]
sinSeries2 =                        forall a. Fractional a => [a] -> [a]
integrateSeries forall a. Fractional a => [a]
cosSeries2

-- | Power series expansion of @cosh(x)@
coshSeries :: Fractional a => [a]
coshSeries :: forall a. Fractional a => [a]
coshSeries = forall {t}. Fractional t => t -> t -> [t]
go a
0 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
e forall a. a -> [a] -> [a]
: t
0 forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
2) (t
e forall a. Fractional a => a -> a -> a
/ ((t
iforall a. Num a => a -> a -> a
+t
1)forall a. Num a => a -> a -> a
*(t
iforall a. Num a => a -> a -> a
+t
2)))

-- | Power series expansion of @sinh(x)@
sinhSeries :: Fractional a => [a]
sinhSeries :: forall a. Fractional a => [a]
sinhSeries = forall {t}. Fractional t => t -> t -> [t]
go a
1 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
0 forall a. a -> [a] -> [a]
: t
e forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
2) (t
e forall a. Fractional a => a -> a -> a
/ ((t
iforall a. Num a => a -> a -> a
+t
1)forall a. Num a => a -> a -> a
*(t
iforall a. Num a => a -> a -> a
+t
2)))

-- | Power series expansion of @log(1+x)@
log1Series :: Fractional a => [a]
log1Series :: forall a. Fractional a => [a]
log1Series = a
0 forall a. a -> [a] -> [a]
: forall {t}. Fractional t => t -> t -> [t]
go a
1 a
1 where
  go :: t -> t -> [t]
go t
i t
e = (t
eforall a. Fractional a => a -> a -> a
/t
i) forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
iforall a. Num a => a -> a -> a
+t
1) (-t
e)

-- | Power series expansion of @(1-Sqrt[1-4x])/(2x)@ (the coefficients are the Catalan numbers)
dyckSeries :: Num a => [a]
dyckSeries :: forall a. Num a => [a]
dyckSeries = [ forall a. Num a => Integer -> a
fromInteger (forall a. Integral a => a -> Integer
catalan Int
i) | Int
i<-[(Int
0::Int)..] ]

--------------------------------------------------------------------------------
-- * \"Coin\" series

-- | Power series expansion of 
-- 
-- > 1 / ( (1-x^k_1) * (1-x^k_2) * ... * (1-x^k_n) )
--
-- Example:
--
-- @(coinSeries [2,3,5])!!k@ is the number of ways 
-- to pay @k@ dollars with coins of two, three and five dollars.
--
-- TODO: better name?
coinSeries :: [Int] -> [Integer]
coinSeries :: [Int] -> [Integer]
coinSeries [] = Integer
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat Integer
0
coinSeries (Int
k:[Int]
ks) = [Integer]
xs where
  xs :: [Integer]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) ([Int] -> [Integer]
coinSeries [Int]
ks) (forall a. Int -> a -> [a]
replicate Int
k Integer
0 forall a. [a] -> [a] -> [a]
++ [Integer]
xs) 

-- | Generalization of the above to include coefficients: expansion of 
--  
-- > 1 / ( (1-a_1*x^k_1) * (1-a_2*x^k_2) * ... * (1-a_n*x^k_n) ) 
-- 
coinSeries' :: Num a => [(a,Int)] -> [a]
coinSeries' :: forall a. Num a => [(a, Int)] -> [a]
coinSeries' [] = a
1 forall a. a -> [a] -> [a]
: forall a. a -> [a]
repeat a
0
coinSeries' ((a
a,Int
k):[(a, Int)]
aks) = [a]
xs where
  xs :: [a]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) (forall a. Num a => [(a, Int)] -> [a]
coinSeries' [(a, Int)]
aks) (forall a. Int -> a -> [a]
replicate Int
k a
0 forall a. [a] -> [a] -> [a]
++ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
*a
a) [a]
xs) 

convolveWithCoinSeries :: [Int] -> [Integer] -> [Integer]
convolveWithCoinSeries :: [Int] -> [Integer] -> [Integer]
convolveWithCoinSeries [Int]
ks [Integer]
series1 = [Int] -> [Integer]
worker [Int]
ks where
  series :: [Integer]
series = [Integer]
series1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat Integer
0
  worker :: [Int] -> [Integer]
worker [] = [Integer]
series
  worker (Int
k:[Int]
ks) = [Integer]
xs where
    xs :: [Integer]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) ([Int] -> [Integer]
worker [Int]
ks) (forall a. Int -> a -> [a]
replicate Int
k Integer
0 forall a. [a] -> [a] -> [a]
++ [Integer]
xs)

convolveWithCoinSeries' :: Num a => [(a,Int)] -> [a] -> [a]
convolveWithCoinSeries' :: forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithCoinSeries' [(a, Int)]
ks [a]
series1 = [(a, Int)] -> [a]
worker [(a, Int)]
ks where
  series :: [a]
series = [a]
series1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0
  worker :: [(a, Int)] -> [a]
worker [] = [a]
series
  worker ((a
a,Int
k):[(a, Int)]
aks) = [a]
xs where
    xs :: [a]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) ([(a, Int)] -> [a]
worker [(a, Int)]
aks) (forall a. Int -> a -> [a]
replicate Int
k a
0 forall a. [a] -> [a] -> [a]
++ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
*a
a) [a]
xs)

--------------------------------------------------------------------------------
-- * Reciprocals of products of polynomials

-- | Convolution of many 'pseries', that is, the expansion of the reciprocal
-- of a product of polynomials
productPSeries :: [[Int]] -> [Integer]
productPSeries :: [[Int]] -> [Integer]
productPSeries = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (forall a b c. (a -> b -> c) -> b -> a -> c
flip [Int] -> [Integer] -> [Integer]
convolveWithPSeries) forall a. Num a => [a]
unitSeries

-- | The same, with coefficients.
productPSeries' :: Num a => [[(a,Int)]] -> [a]
productPSeries' :: forall a. Num a => [[(a, Int)]] -> [a]
productPSeries' = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithPSeries') forall a. Num a => [a]
unitSeries

convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]
convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]
convolveWithProductPSeries [[Int]]
kss [Integer]
ser = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (forall a b c. (a -> b -> c) -> b -> a -> c
flip [Int] -> [Integer] -> [Integer]
convolveWithPSeries) [Integer]
ser [[Int]]
kss

-- | This is the most general function in this module; all the others
-- are special cases of this one.  
convolveWithProductPSeries' :: Num a => [[(a,Int)]] -> [a] -> [a] 
convolveWithProductPSeries' :: forall a. Num a => [[(a, Int)]] -> [a] -> [a]
convolveWithProductPSeries' [[(a, Int)]]
akss [a]
ser = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithPSeries') [a]
ser [[(a, Int)]]
akss
  
--------------------------------------------------------------------------------
-- * Reciprocals of polynomials

-- Reciprocals of polynomials, without coefficients

-- | The power series expansion of 
--
-- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)
--
pseries :: [Int] -> [Integer]
pseries :: [Int] -> [Integer]
pseries [Int]
ks = [Int] -> [Integer] -> [Integer]
convolveWithPSeries [Int]
ks forall a. Num a => [a]
unitSeries

-- | Convolve with (the expansion of) 
--
-- > 1 / (1 - x^k_1 - x^k_2 - x^k_3 - ... - x^k_n)
--
convolveWithPSeries :: [Int] -> [Integer] -> [Integer]
convolveWithPSeries :: [Int] -> [Integer] -> [Integer]
convolveWithPSeries [Int]
ks [Integer]
series1 = [Integer]
ys where 
  series :: [Integer]
series = [Integer]
series1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat Integer
0 
  ys :: [Integer]
ys = [Int] -> [Integer] -> [Integer]
worker [Int]
ks [Integer]
ys 
  worker :: [Int] -> [Integer] -> [Integer]
worker [] [Integer]
_ = [Integer]
series 
  worker (Int
k:[Int]
ks) [Integer]
ys = [Integer]
xs where
    xs :: [Integer]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) (forall a. Int -> a -> [a]
replicate Int
k Integer
0 forall a. [a] -> [a] -> [a]
++ [Integer]
ys) ([Int] -> [Integer] -> [Integer]
worker [Int]
ks [Integer]
ys)

--------------------------------------------------------------------------------
--  Reciprocals of polynomials, with coefficients

-- | The expansion of 
--
-- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)
--
pseries' :: Num a => [(a,Int)] -> [a]
pseries' :: forall a. Num a => [(a, Int)] -> [a]
pseries' [(a, Int)]
aks = forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithPSeries' [(a, Int)]
aks forall a. Num a => [a]
unitSeries

-- | Convolve with (the expansion of) 
--
-- > 1 / (1 - a_1*x^k_1 - a_2*x^k_2 - a_3*x^k_3 - ... - a_n*x^k_n)
--
convolveWithPSeries' :: Num a => [(a,Int)] -> [a] -> [a]
convolveWithPSeries' :: forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithPSeries' [(a, Int)]
aks [a]
series1 = [a]
ys where 
  series :: [a]
series = [a]
series1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat a
0 
  ys :: [a]
ys = [(a, Int)] -> [a] -> [a]
worker [(a, Int)]
aks [a]
ys 
  worker :: [(a, Int)] -> [a] -> [a]
worker [] [a]
_ = [a]
series
  worker ((a
a,Int
k):[(a, Int)]
aks) [a]
ys = [a]
xs where
    xs :: [a]
xs = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) (forall a. Int -> a -> [a]
replicate Int
k a
0 forall a. [a] -> [a] -> [a]
++ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
*a
a) [a]
ys) ([(a, Int)] -> [a] -> [a]
worker [(a, Int)]
aks [a]
ys)

{-
data Sign = Plus | Minus deriving (Eq,Show)

signValue :: Num a => Sign -> a
signValue Plus  =  1
signValue Minus = -1
-}

signedPSeries :: [(Sign,Int)] -> [Integer] 
signedPSeries :: [(Sign, Int)] -> [Integer]
signedPSeries [(Sign, Int)]
aks = [(Sign, Int)] -> [Integer] -> [Integer]
convolveWithSignedPSeries [(Sign, Int)]
aks forall a. Num a => [a]
unitSeries

-- | Convolve with (the expansion of) 
--
-- > 1 / (1 +- x^k_1 +- x^k_2 +- x^k_3 +- ... +- x^k_n)
--
-- Should be faster than using `convolveWithPSeries'`.
-- Note: 'Plus' corresponds to the coefficient @-1@ in `pseries'` (since
-- there is a minus sign in the definition there)!
convolveWithSignedPSeries :: [(Sign,Int)] -> [Integer] -> [Integer]
convolveWithSignedPSeries :: [(Sign, Int)] -> [Integer] -> [Integer]
convolveWithSignedPSeries [(Sign, Int)]
aks [Integer]
series1 = [Integer]
ys where 
  series :: [Integer]
series = [Integer]
series1 forall a. [a] -> [a] -> [a]
++ forall a. a -> [a]
repeat Integer
0 
  ys :: [Integer]
ys = [(Sign, Int)] -> [Integer] -> [Integer]
worker [(Sign, Int)]
aks [Integer]
ys 
  worker :: [(Sign, Int)] -> [Integer] -> [Integer]
worker [] [Integer]
_ = [Integer]
series
  worker ((Sign
a,Int
k):[(Sign, Int)]
aks) [Integer]
ys = [Integer]
xs where
    xs :: [Integer]
xs = case Sign
a of
      Sign
Minus -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) [Integer]
one [Integer]
two 
      Sign
Plus  -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Integer]
one [Integer]
two
    one :: [Integer]
one = [(Sign, Int)] -> [Integer] -> [Integer]
worker [(Sign, Int)]
aks [Integer]
ys
    two :: [Integer]
two = forall a. Int -> a -> [a]
replicate Int
k Integer
0 forall a. [a] -> [a] -> [a]
++ [Integer]
ys
     
--------------------------------------------------------------------------------