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