-- | 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 :: [a]
unitSeries = a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
forall a. a -> [a]
repeat a
0

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

-- | Power series representing a constant function
constSeries :: Num a => a -> [a]
constSeries :: a -> [a]
constSeries a
x = a
x a -> [a] -> [a]
forall a. a -> [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 :: [a]
idSeries = a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
forall a. a -> [a]
repeat a
0

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

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

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

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

subSeries :: Num a => [a] -> [a] -> [a]
subSeries :: [a] -> [a] -> [a]
subSeries [a]
xs [a]
ys = a -> a -> (a -> a -> a) -> [a] -> [a] -> [a]
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 :: [a] -> [a]
negateSeries = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
negate

scaleSeries :: Num a => a -> [a] -> [a]
scaleSeries :: a -> [a] -> [a]
scaleSeries a
s = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
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 :: [a] -> [a] -> [a]
mulSeries [a]
xs [a]
ys = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
go ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0) ([a]
ys [a] -> [a] -> [a]
forall a. [a] -> [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
fa -> a -> a
forall a. Num a => a -> a -> a
*a
g a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleSeries a
f [a]
gs) [a] -> [a] -> [a]
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 :: [a] -> [a] -> [a]
mulSeriesNaive = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
convolve

productOfSeries :: Num a => [[a]] -> [a]
productOfSeries :: [[a]] -> [a]
productOfSeries = [[a]] -> [a]
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 :: [a] -> [a] -> [a]
convolve [a]
xs1 [a]
ys1 = [a]
res where
  res :: [a]
res = [ (a -> a -> a) -> a -> [a] -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0 ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
xs ([a] -> [a]
forall a. [a] -> [a]
reverse (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n [a]
ys)))
        | Int
n<-[Int
1..] 
        ]
  xs :: [a]
xs = [a]
xs1 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0
  ys :: [a]
ys = [a]
ys1 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0

-- | Convolution (= product) of many series. Still slow!
convolveMany :: Num a => [[a]] -> [a]
convolveMany :: [[a]] -> [a]
convolveMany []  = a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a]
forall a. a -> [a]
repeat a
0
convolveMany [[a]]
xss = ([a] -> [a] -> [a]) -> [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 [a] -> [a] -> [a]
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 :: [a] -> [a] -> [a]
divSeries [a]
xs [a]
ys = [a] -> [a] -> [a]
forall a. (Eq a, Fractional a) => [a] -> [a] -> [a]
go ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0) ([a]
ys [a] -> [a] -> [a]
forall a. [a] -> [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
fa -> a -> a
forall a. Fractional a => a -> a -> a
/a
g in a
q a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
go ([a]
fs [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
`subSeries` a -> [a] -> [a]
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 :: [a] -> [a]
reciprocalSeries [a]
series = case [a]
series of
  [] -> [Char] -> [a]
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 -> [Char] -> [a]
forall a. HasCallStack => [Char] -> a
error [Char]
"reciprocalSeries: input series has constant term 0"
    a
_ -> (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/a
a) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. (Eq a, Num a) => [a] -> [a]
integralReciprocalSeries ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
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 :: [a] -> [a]
integralReciprocalSeries [a]
series = case [a]
series of 
  [] -> [Char] -> [a]
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 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
worker [a
1]
    a
_ -> [Char] -> [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' = - [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
series) [a]
bs) 
                in  a
b' a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
worker (a
b'a -> [a] -> [a]
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 :: [a] -> [a] -> [a]
composeSeries [a]
xs [a]
ys = [a] -> [a] -> [a]
forall a. (Eq a, Num a) => [a] -> [a] -> [a]
go ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0) ([a]
ys [a] -> [a] -> [a]
forall a. [a] -> [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 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
mulSeries [a]
gs ([a] -> [a] -> [a]
go [a]
fs (a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
gs))
  go (a
f:[a]
fs) (a
_:[a]
gs) = [Char] -> [a]
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 :: [a] -> [a] -> [a]
substitute [a]
f [a]
g = [a] -> [a] -> [a]
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 :: [a] -> [a] -> [a]
composeSeriesNaive [a]
g [a]
f = [a] -> [a] -> [a]
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 :: [a] -> [a] -> [a]
substituteNaive [a]
as_ [a]
bs_ = 
  case [a] -> a
forall a. [a] -> a
head [a]
as of
    a
0 -> [ Int -> a
f Int
n | Int
n<-[Int
0..] ]
    a
_ -> [Char] -> [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_ [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0
    bs :: [a]
bs = [a]
bs_ [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0
    a :: Int -> a
a Int
i = [a]
as [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
i
    b :: Int -> a
b Int
j = [a]
bs [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
j
    f :: Int -> a
f Int
n = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
            [ Int -> a
b Int
m a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^Int
j | (Int
i,Int
j)<-[(Int, Int)]
es ] a -> a -> a
forall a. Num a => a -> a -> a
* Integer -> a
forall a. Num a => Integer -> a
fromInteger ([Int] -> Integer
forall a. Integral a => [a] -> Integer
multinomial (((Int, Int) -> Int) -> [(Int, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> Int
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 :: [a] -> [a]
lagrangeInversion [a]
xs = [a] -> [a]
forall a. (Eq a, Fractional a) => [a] -> [a]
go ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [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 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. (Eq a, Fractional a) => [a] -> [a] -> [a]
divSeries [a]
forall a. Num a => [a]
unitSeries ([a] -> [a] -> [a]
forall a. (Eq a, Num a) => [a] -> [a] -> [a]
composeSeries [a]
fs [a]
rs)
  go (a
_:[a]
fs) = [Char] -> [a]
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 = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
div Integer
numer Integer
denom where
  numer :: Integer
numer = (-Integer
1)Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral [Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
m])
  denom :: Integer
denom = Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (((Int, Int) -> Integer) -> [(Int, Int)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Integer
forall a. Integral a => a -> Integer
factorial (Int -> Integer) -> ((Int, Int) -> Int) -> (Int, Int) -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Int) -> Int
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 :: [a] -> [a]
integralLagrangeInversionNaive [a]
series_ = 
  case [a]
series of
    (a
0:a
1:[a]
rest) -> a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [ Int -> a
f Int
n | Int
n<-[Int
1..] ]
    [a]
_ -> [Char] -> [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_ [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0
    as :: [a]
as  = [a] -> [a]
forall a. [a] -> [a]
tail [a]
series 
    a :: Int -> a
a Int
i = [a]
as [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
i
    f :: Int -> a
f Int
n = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Integer -> a
forall a. Num a => Integer -> a
fromInteger (Partition -> Integer
lagrangeCoeff Partition
p) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)a -> Int -> a
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 :: [a] -> [a]
lagrangeInversionNaive [a]
series_ = 
  case [a]
series of
    (a
0:a
a1:[a]
rest) -> if a
a1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0 
      then [a]
forall a. a
err 
      else a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
a1) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [ Int -> a
f Int
n a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
a1a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) | Int
n<-[Int
1..] ]
    [a]
_ -> [a]
forall a. a
err
  where
    err :: a
err    = [Char] -> a
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_ [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
0
    a1 :: a
a1  = [a]
series [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
1
    as :: [a]
as  = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/a
a1) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
series)
    a :: Int -> a
a Int
i = [a]
as [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
i
    f :: Int -> a
f Int
n = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Integer -> a
forall a. Num a => Integer -> a
fromInteger (Partition -> Integer
lagrangeCoeff Partition
p) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ (Int -> a
a Int
i)a -> Int -> a
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 :: [a] -> [a]
differentiateSeries (a
y:[a]
ys) = Int -> [a] -> [a]
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) = t -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral t
n a -> a -> a
forall a. Num a => a -> a -> a
* a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: t -> [a] -> [a]
go (t
nt -> t -> t
forall a. Num a => a -> a -> a
+t
1) [a]
xs
  go t
_  []     = []

integrateSeries :: Fractional a => [a] -> [a]
integrateSeries :: [a] -> [a]
integrateSeries [a]
ys = a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [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 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (t -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral t
n) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: t -> [a] -> [a]
go (t
nt -> t -> t
forall 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 :: [a]
expSeries = a -> a -> [a]
forall t. Fractional t => t -> t -> [t]
go a
0 a
1 where
  go :: t -> t -> [t]
go t
i t
e = t
e t -> [t] -> [t]
forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
it -> t -> t
forall a. Num a => a -> a -> a
+t
1) (t
e t -> t -> t
forall a. Fractional a => a -> a -> a
/ (t
it -> t -> t
forall a. Num a => a -> a -> a
+t
1))

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

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

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

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

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

-- | Power series expansion of @log(1+x)@
log1Series :: Fractional a => [a]
log1Series :: [a]
log1Series = a
0 a -> [a] -> [a]
forall a. a -> [a] -> [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
et -> t -> t
forall a. Fractional a => a -> a -> a
/t
i) t -> [t] -> [t]
forall a. a -> [a] -> [a]
: t -> t -> [t]
go (t
it -> t -> t
forall 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 :: [a]
dyckSeries = [ Integer -> a
forall a. Num a => Integer -> a
fromInteger (Int -> Integer
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 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> [Integer]
forall a. a -> [a]
repeat Integer
0
coinSeries (Int
k:[Int]
ks) = [Integer]
xs where
  xs :: [Integer]
xs = (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) ([Int] -> [Integer]
coinSeries [Int]
ks) (Int -> Integer -> [Integer]
forall a. Int -> a -> [a]
replicate Int
k Integer
0 [Integer] -> [Integer] -> [Integer]
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' :: [(a, Int)] -> [a]
coinSeries' [] = a
1 a -> [a] -> [a]
forall a. a -> [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 = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) ([(a, Int)] -> [a]
forall a. Num a => [(a, Int)] -> [a]
coinSeries' [(a, Int)]
aks) (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
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 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ Integer -> [Integer]
forall a. a -> [a]
repeat Integer
0
  worker :: [Int] -> [Integer]
worker [] = [Integer]
series
  worker (Int
k:[Int]
ks) = [Integer]
xs where
    xs :: [Integer]
xs = (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) ([Int] -> [Integer]
worker [Int]
ks) (Int -> Integer -> [Integer]
forall a. Int -> a -> [a]
replicate Int
k Integer
0 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ [Integer]
xs)

convolveWithCoinSeries' :: Num a => [(a,Int)] -> [a] -> [a]
convolveWithCoinSeries' :: [(a, Int)] -> [a] -> [a]
convolveWithCoinSeries' [(a, Int)]
ks [a]
series1 = [(a, Int)] -> [a]
worker [(a, Int)]
ks where
  series :: [a]
series = [a]
series1 [a] -> [a] -> [a]
forall a. [a] -> [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 = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) ([(a, Int)] -> [a]
worker [(a, Int)]
aks) (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
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 = ([Integer] -> [Int] -> [Integer])
-> [Integer] -> [[Int]] -> [Integer]
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (([Int] -> [Integer] -> [Integer])
-> [Integer] -> [Int] -> [Integer]
forall a b c. (a -> b -> c) -> b -> a -> c
flip [Int] -> [Integer] -> [Integer]
convolveWithPSeries) [Integer]
forall a. Num a => [a]
unitSeries

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

convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]
convolveWithProductPSeries :: [[Int]] -> [Integer] -> [Integer]
convolveWithProductPSeries [[Int]]
kss [Integer]
ser = ([Integer] -> [Int] -> [Integer])
-> [Integer] -> [[Int]] -> [Integer]
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (([Int] -> [Integer] -> [Integer])
-> [Integer] -> [Int] -> [Integer]
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' :: [[(a, Int)]] -> [a] -> [a]
convolveWithProductPSeries' [[(a, Int)]]
akss [a]
ser = ([a] -> [(a, Int)] -> [a]) -> [a] -> [[(a, Int)]] -> [a]
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (([(a, Int)] -> [a] -> [a]) -> [a] -> [(a, Int)] -> [a]
forall a b c. (a -> b -> c) -> b -> a -> c
flip [(a, Int)] -> [a] -> [a]
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 [Integer]
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 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ Integer -> [Integer]
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 = (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) (Int -> Integer -> [Integer]
forall a. Int -> a -> [a]
replicate Int
k Integer
0 [Integer] -> [Integer] -> [Integer]
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' :: [(a, Int)] -> [a]
pseries' [(a, Int)]
aks = [(a, Int)] -> [a] -> [a]
forall a. Num a => [(a, Int)] -> [a] -> [a]
convolveWithPSeries' [(a, Int)]
aks [a]
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' :: [(a, Int)] -> [a] -> [a]
convolveWithPSeries' [(a, Int)]
aks [a]
series1 = [a]
ys where 
  series :: [a]
series = [a]
series1 [a] -> [a] -> [a]
forall a. [a] -> [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 = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
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 [Integer]
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 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ Integer -> [Integer]
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 -> (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(+) [Integer]
one [Integer]
two 
      Sign
Plus  -> (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
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 = Int -> Integer -> [Integer]
forall a. Int -> a -> [a]
replicate Int
k Integer
0 [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ [Integer]
ys
     
--------------------------------------------------------------------------------