{-# LANGUAGE RebindableSyntax #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Data.YAP.PowerSeries
-- Copyright   :  (c) Ross Paterson 2021
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  [email protected]
-- Stability   :  provisional
-- Portability :  portable
--
-- An example instance of the algebraic classes: formal power series.
--
-- * M. Douglas McIlroy, "Power series, power serious",
--   /Journal of Functional Programming/ 9:3 (1999) 325-337.
--   <https://www.cs.dartmouth.edu/~doug/powser.html>
--
-- * M. Douglas McIlroy, "The music of streams",
--   /Information Processing Letters/ 77 (2001) 189-195.
--
-----------------------------------------------------------------------------

module Data.YAP.PowerSeries (
    -- * Formal power series
    PowerSeries,
    constant,
    fromPolynomial,
    fromCoefficients,
    fromDerivatives,

    -- * Queries
    atZero,
    coefficients,
    derivatives,
    flatten,
    toPolynomial,
    modulus,
    order,
    distance,

    -- * Approximations
    approximations,
    approximationsWith,
    padeApproximants,

    -- * Composition
    identity,
    compose,
    inverse,
    inverseSimple,
    -- ** Special compositions
    (.*),
    (.^),

    -- * Other transforms
    divX,
    mulX,
    shift,
    recipSimple,
    divSimple,

    -- * Special series
    recipOneMinus,
    expS,
    logOnePlus,
    powerOnePlus,
    euler,
    -- ** Trigonometric functions
    sinS,
    cosS,
    secS,
    tanS,
    asinS,
    atanS,
    -- ** Hyperbolic functions
    sinhS,
    coshS,
    sechS,
    tanhS,
    asinhS,
    atanhS,

    -- * Continued fractions
    stieltjes,
    jacobi,

    -- * Sequence transforms
    -- ** Binomial transform
    binomialTransform,
    inverseBinomialTransform,

    -- ** Lambert transform
    lambertTransform,
    inverseLambertTransform,

    -- ** Euler transform
    eulerTransform,
    inverseEulerTransform,

    -- * Ordinary generating functions
    -- ** Number sequences
    -- $ogfs

    -- ** Linear recurrences
    linearRecurrence,

    -- ** Formal languages
    GradedLanguage,
    string,

    -- * Polynomial sequences
    -- $polynomial_sequences
    constants,
    monomials,
    riordan,
    delta,

    -- * Probability generating functions #PGF#
    -- $pgf
    bernoulliPGF,
    geometricPGF,
    poissonPGF,

  ) where

import Prelude.YAP
import Data.YAP.Algebra
import Data.YAP.FiniteSet
import qualified Data.YAP.Polynomial as Poly
import Data.YAP.Polynomial (Polynomial, RationalFunction)
import Data.YAP.Ratio
import Data.YAP.Utilities.List (longZipWith)

import Data.List (intercalate)

infixl 9 .^
infixl 9 .*

-- | Formal power series \(\sum_{n=0}^\infty a_n x^n\), represented by the
-- sequence of coefficients \(a_n\) of powers of the indeterminate \(x\).
newtype PowerSeries a = PS [a]
    -- ^ Coefficients, least significant first.

-- | Value of the function at zero
atZero :: (AdditiveMonoid a) => PowerSeries a -> a
atZero :: forall a. AdditiveMonoid a => PowerSeries a -> a
atZero (PS []) = a
forall a. AdditiveMonoid a => a
zero
atZero (PS (a
a:[a]
_)) = a
a

-- Horner view: coefficients of powers of x

-- | Discard the constant term and divide by the indeterminate
divX :: PowerSeries a -> PowerSeries a
divX :: forall a. PowerSeries a -> PowerSeries a
divX (PS []) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS []
divX (PS (a
_:[a]
as)) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS [a]
as

-- | Multiply by the indeterminate
mulX :: (AdditiveMonoid a) => PowerSeries a -> PowerSeries a
mulX :: forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
f = a -> PowerSeries a -> PowerSeries a
forall a. a -> PowerSeries a -> PowerSeries a
plusMulX a
forall a. AdditiveMonoid a => a
zero PowerSeries a
f

-- | Multiply by the indeterminate and add @a@
plusMulX :: a -> PowerSeries a -> PowerSeries a
plusMulX :: forall a. a -> PowerSeries a -> PowerSeries a
plusMulX a
a (PS [a]
as) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS (a
aa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
as)

-- | Shift the power series by \(k\) (which may be negative): multiply
-- by \(x^k\) and discard any terms with negative exponents.
shift :: (AdditiveMonoid a) => Int -> PowerSeries a -> PowerSeries a
shift :: forall a. AdditiveMonoid a => Int -> PowerSeries a -> PowerSeries a
shift Int
_ (PS []) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS []
shift Int
k (PS [a]
as)
  | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
k a
forall a. AdditiveMonoid a => a
zero [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
as)
  | Bool
otherwise = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (-Int
k) [a]
as)

-- | The coefficients of the power series, least significant first.
coefficients :: (AdditiveMonoid a) => PowerSeries a -> [a]
coefficients :: forall a. AdditiveMonoid a => PowerSeries a -> [a]
coefficients (PS [a]
as) = [a]
as [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
forall a. AdditiveMonoid a => a
zero

-- | Power series formed from a list of coefficients.
-- If the list is finite, the remaining coefficients are zero.
fromCoefficients :: (AdditiveMonoid a) => [a] -> PowerSeries a
fromCoefficients :: forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS

-- Maclaurin view: values of repeated derivatives at zero

instance (Semiring a) => Differentiable (PowerSeries a) where
    derivative :: PowerSeries a -> PowerSeries a
derivative (PS []) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS []
    derivative (PS (a
_:[a]
as)) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS ((Int -> a -> a) -> [Int] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> a -> a
forall b. ToInteger b => b -> a -> a
forall a b. (AdditiveMonoid a, ToInteger b) => b -> a -> a
atimes (Int -> [Int]
forall a. Semiring a => a -> [a]
from (Int
1::Int)) [a]
as)

instance (FromRational a) => Integrable (PowerSeries a) where
    integral :: PowerSeries a -> PowerSeries a
integral (PS [a]
as) = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS (a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:(a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Semiring a => a -> a -> a
(*) ((Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> a
forall {a}. FromRational a => Integer -> a
coeff (Integer -> [Integer]
forall a. Semiring a => a -> [a]
from Integer
1)) [a]
as)
      where
        coeff :: Integer -> a
coeff Integer
n = Rational -> a
forall a. FromRational a => Rational -> a
fromRational (Integer
1Integer -> Integer -> Rational
forall a.
(Eq a, StandardAssociate a, Euclidean a) =>
a -> a -> Ratio a
%Integer
n)

instance AdditiveFunctor PowerSeries where
    mapAdditive :: forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
mapAdditive a -> b
f (PS [a]
as) = [b] -> PowerSeries b
forall a. [a] -> PowerSeries a
PS ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
as)

-- | A list whose \(n\)th entry is the value of the \(n\)th derivative at zero.
derivatives :: (Semiring a) => PowerSeries a -> [a]
derivatives :: forall a. Semiring a => PowerSeries a -> [a]
derivatives PowerSeries a
p = (Integer -> a -> a) -> [Integer] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> a -> a
forall b. ToInteger b => b -> a -> a
forall a b. (AdditiveMonoid a, ToInteger b) => b -> a -> a
atimes [Integer]
factorials (PowerSeries a -> [a]
forall a. AdditiveMonoid a => PowerSeries a -> [a]
coefficients PowerSeries a
p)
  where
    factorials :: [Integer]
factorials = (Integer -> Integer -> Integer)
-> Integer -> [Integer] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Integer -> Integer -> Integer
forall a. Semiring a => a -> a -> a
(*) (Integer
1::Integer) [Integer
1..]

-- | Power series formed from a list of derivatives at zero.
-- If the list is finite, the remaining derivatives are zero.
fromDerivatives :: (FromRational a) => [a] -> PowerSeries a
fromDerivatives :: forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives [a]
as = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Semiring a => a -> a -> a
(*) ((Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> a
forall {a}. FromRational a => Integer -> a
coeff [Integer]
factorials) [a]
as)
  where
    coeff :: Integer -> a
coeff Integer
n = Rational -> a
forall a. FromRational a => Rational -> a
fromRational (Integer
1Integer -> Integer -> Rational
forall a.
(Eq a, StandardAssociate a, Euclidean a) =>
a -> a -> Ratio a
%Integer
n)
    factorials :: [Integer]
factorials = (Integer -> Integer -> Integer)
-> Integer -> [Integer] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Integer -> Integer -> Integer
forall a. Semiring a => a -> a -> a
(*) (Integer
1::Integer) [Integer
1..]

instance (AdditiveMonoid a) => AdditiveMonoid (PowerSeries a) where
    PS [a]
as + :: PowerSeries a -> PowerSeries a -> PowerSeries a
+ PS [a]
bs = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS ([a] -> [a] -> [a]
forall a. AdditiveMonoid a => [a] -> [a] -> [a]
add [a]
as [a]
bs)
    zero :: PowerSeries a
zero = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS []
    atimes :: forall b. ToInteger b => b -> PowerSeries a -> PowerSeries a
atimes b
n PowerSeries a
p
      | b
n b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== b
forall a. AdditiveMonoid a => a
zero = PowerSeries a
forall a. AdditiveMonoid a => a
zero
      | Bool
otherwise = (a -> a) -> PowerSeries a -> PowerSeries a
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive (b -> a -> a
forall b. ToInteger b => b -> a -> a
forall a b. (AdditiveMonoid a, ToInteger b) => b -> a -> a
atimes b
n) PowerSeries a
p

add :: (AdditiveMonoid a) => [a] -> [a] -> [a]
add :: forall a. AdditiveMonoid a => [a] -> [a] -> [a]
add = (a -> a -> a) -> [a] -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a] -> [a]
longZipWith a -> a -> a
forall a. AdditiveMonoid a => a -> a -> a
(+)

instance (AbelianGroup a) => AbelianGroup (PowerSeries a) where
    negate :: PowerSeries a -> PowerSeries a
negate = (a -> a) -> PowerSeries a -> PowerSeries a
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive a -> a
forall a. AbelianGroup a => a -> a
negate
    gtimes :: forall b.
(AbelianGroup b, ToInteger b) =>
b -> PowerSeries a -> PowerSeries a
gtimes b
n PowerSeries a
p
      | b
n b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== b
forall a. AdditiveMonoid a => a
zero = PowerSeries a
forall a. AdditiveMonoid a => a
zero
      | Bool
otherwise = (a -> a) -> PowerSeries a -> PowerSeries a
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive (b -> a -> a
forall b. (AbelianGroup b, ToInteger b) => b -> a -> a
forall a b.
(AbelianGroup a, AbelianGroup b, ToInteger b) =>
b -> a -> a
gtimes b
n) PowerSeries a
p

instance (Semiring a) => Semiring (PowerSeries a) where
    PS [a]
as * :: PowerSeries a -> PowerSeries a -> PowerSeries a
* PS [a]
bs = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS ([a] -> [a] -> [a]
forall a. Semiring a => [a] -> [a] -> [a]
mul [a]
as [a]
bs)
    one :: PowerSeries a
one = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant a
forall a. Semiring a => a
one
    fromNatural :: Natural -> PowerSeries a
fromNatural Natural
i = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant (Natural -> a
forall a. Semiring a => Natural -> a
fromNatural Natural
i)

-- convolution
mul :: (Semiring a) => [a] -> [a] -> [a]
mul :: forall a. Semiring a => [a] -> [a] -> [a]
mul [a]
_ [] = []
mul [a]
xs (a
y:[a]
ys) = (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> [a] -> [a]
mulOne [] [a]
xs
  where
    mulOne :: a -> [a] -> [a]
mulOne a
x [a]
zs = (a
xa -> a -> a
forall a. Semiring a => a -> a -> a
*a
y) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. AdditiveMonoid a => [a] -> [a] -> [a]
add ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
xa -> a -> a
forall a. Semiring a => a -> a -> a
*) [a]
ys) [a]
zs

instance (Ring a) => Ring (PowerSeries a) where
    fromInteger :: Integer -> PowerSeries a
fromInteger Integer
i = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant (Integer -> a
forall a. Ring a => Integer -> a
fromInteger Integer
i)

instance (FromRational a) => FromRational (PowerSeries a) where
    fromRational :: Rational -> PowerSeries a
fromRational Rational
x = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant (Rational -> a
forall a. FromRational a => Rational -> a
fromRational Rational
x)

-- | Units have non-zero leading coefficients.
-- Standard associates are monomials of the form \(x^n\).
-- 'stdUnit' and 'stdRecip' are not defined on zero.
instance (Eq a, Field a) => StandardAssociate (PowerSeries a) where
    stdAssociate :: PowerSeries a -> PowerSeries a
stdAssociate = PowerSeries a -> PowerSeries a
forall a. (Eq a, Semiring a) => PowerSeries a -> PowerSeries a
stdAssociateStream
    stdUnit :: PowerSeries a -> PowerSeries a
stdUnit = (PowerSeries a -> Bool)
-> (PowerSeries a -> PowerSeries a)
-> PowerSeries a
-> PowerSeries a
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) (a -> Bool) -> (PowerSeries a -> a) -> PowerSeries a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero) PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX
    stdRecip :: PowerSeries a -> PowerSeries a
stdRecip PowerSeries a
f = PowerSeries a
1 PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Field a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
`exactDiv` PowerSeries a -> PowerSeries a
forall a. StandardAssociate a => a -> a
stdUnit PowerSeries a
f

stdAssociateStream :: (Eq a, Semiring a) => PowerSeries a -> PowerSeries a
stdAssociateStream :: forall a. (Eq a, Semiring a) => PowerSeries a -> PowerSeries a
stdAssociateStream PowerSeries a
f
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
f a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
forall a. AdditiveMonoid a => a
zero = PowerSeries a
forall a. Semiring a => a
one
  | Bool
otherwise = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX (PowerSeries a -> PowerSeries a
forall a. (Eq a, Semiring a) => PowerSeries a -> PowerSeries a
stdAssociateStream (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f))

-- | @'mod' x y@ is non-zero if @y@ has more leading zeros than @x@.
-- See also 'modulus', for the remainder as a polynomial.
instance (Eq a, Field a) => Euclidean (PowerSeries a) where
    divMod :: PowerSeries a -> PowerSeries a -> (PowerSeries a, PowerSeries a)
divMod PowerSeries a
f PowerSeries a
g = ([a] -> PowerSeries a)
-> (PowerSeries a, [a]) -> (PowerSeries a, PowerSeries a)
forall a b. (a -> b) -> (PowerSeries a, a) -> (PowerSeries a, b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
divModStream PowerSeries a
f PowerSeries a
g)
    div :: PowerSeries a -> PowerSeries a -> PowerSeries a
div = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divStream
    mod :: PowerSeries a -> PowerSeries a -> PowerSeries a
mod PowerSeries a
f PowerSeries a
g = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (PowerSeries a -> PowerSeries a -> [a]
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> [a]
modStream PowerSeries a
f PowerSeries a
g)

    euclideanNorm :: PowerSeries a -> Natural
euclideanNorm (PS [a]
as) = Int -> Natural
forall a. Euclidean a => a -> Natural
euclideanNorm ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0) [a]
as))

divModStream :: (Eq a, Field a) =>
    PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
divModStream :: forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
divModStream PowerSeries a
f PowerSeries a
g
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = (PowerSeries a
q, PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
fa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
r)
  | Bool
otherwise = (PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Field a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
exactDiv PowerSeries a
f PowerSeries a
g, [])
  where
    (PowerSeries a
q, [a]
r) = PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> (PowerSeries a, [a])
divModStream (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g)

divStream :: (Eq a, Field a) => PowerSeries a -> PowerSeries a -> PowerSeries a
divStream :: forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divStream PowerSeries a
f PowerSeries a
g
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
(Eq a, Field a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divStream (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g)
  | Bool
otherwise = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Field a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
exactDiv PowerSeries a
f PowerSeries a
g

modStream :: (Eq a, AdditiveMonoid a) => PowerSeries a -> PowerSeries a -> [a]
modStream :: forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> [a]
modStream PowerSeries a
f PowerSeries a
g
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall a. AdditiveMonoid a => a
zero = PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
fa -> [a] -> [a]
forall a. a -> [a] -> [a]
:PowerSeries a -> PowerSeries a -> [a]
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> [a]
modStream (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g)
  | Bool
otherwise = []

-- | Reduce a power series \( \sum_{n=0} p_n(x) x^n \) with power series
-- coefficients to a simple power series.
flatten :: (AdditiveMonoid a) => PowerSeries (PowerSeries a) -> PowerSeries a
flatten :: forall a.
AdditiveMonoid a =>
PowerSeries (PowerSeries a) -> PowerSeries a
flatten (PS [PowerSeries a]
ps) = (PowerSeries a -> PowerSeries a -> PowerSeries a)
-> PowerSeries a -> [PowerSeries a] -> PowerSeries a
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
AdditiveMonoid a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
add_shifted PowerSeries a
forall a. AdditiveMonoid a => a
zero [PowerSeries a]
ps
  where
    add_shifted :: PowerSeries a -> PowerSeries a -> PowerSeries a
add_shifted PowerSeries a
p PowerSeries a
q = a -> PowerSeries a -> PowerSeries a
forall a. a -> PowerSeries a -> PowerSeries a
plusMulX (PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
p) (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
p PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => a -> a -> a
+ PowerSeries a
q)

-- | Polynomial with the first \(n\) coefficients, equivalent to the
-- remainder of division by \(x^n\).
toPolynomial :: (AdditiveMonoid a) => Int -> PowerSeries a -> Polynomial a
toPolynomial :: forall a. AdditiveMonoid a => Int -> PowerSeries a -> Polynomial a
toPolynomial Int
n (PS [a]
as) = [a] -> Polynomial a
forall a. [a] -> Polynomial a
Poly.fromCoefficients (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n [a]
as)

-- | The remainder of division by a non-zero power series is a polynomial.
modulus :: (Eq a, AdditiveMonoid a) =>
    PowerSeries a -> PowerSeries a -> Polynomial a
modulus :: forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Polynomial a
modulus PowerSeries a
f PowerSeries a
g = [a] -> Polynomial a
forall a. [a] -> Polynomial a
Poly.fromCoefficients (PowerSeries a -> PowerSeries a -> [a]
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> [a]
modStream PowerSeries a
f PowerSeries a
g)

-- only valid if g0 is non-zero
exactDiv :: (Field a) => PowerSeries a -> PowerSeries a -> PowerSeries a
exactDiv :: forall a.
Field a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
exactDiv PowerSeries a
f PowerSeries a
g = PowerSeries a
q
  where
    q :: PowerSeries a
q = a -> PowerSeries a -> PowerSeries a
forall a. a -> PowerSeries a -> PowerSeries a
plusMulX (a
f0a -> a -> a
forall a. Semifield a => a -> a -> a
/a
g0) (a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant (a -> a
forall a. DivisionSemiring a => a -> a
recip a
g0) PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a -> a
- PowerSeries a
q PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g))
    f0 :: a
f0 = PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
f
    g0 :: a
g0 = PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g

-- | Power series representing a constant value \(c\)
constant :: (AdditiveMonoid a) => a -> PowerSeries a
constant :: forall a. AdditiveMonoid a => a -> PowerSeries a
constant a
a = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS [a
a]

-- | Polynomial considered as a power series
fromPolynomial :: (AdditiveMonoid a) => Polynomial a -> PowerSeries a
fromPolynomial :: forall a. AdditiveMonoid a => Polynomial a -> PowerSeries a
fromPolynomial Polynomial a
p = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (Polynomial a -> [a]
forall a. Polynomial a -> [a]
Poly.unsafeCoefficients Polynomial a
p)

-- | \({1 \over 1-x} = 1 + x + x^2 + \ldots \) (geometric series),
-- converges for \(|x| < 1\).
-- 'derivatives' yields the factorials.
recipOneMinus :: (Semiring a) => PowerSeries a
recipOneMinus :: forall a. Semiring a => PowerSeries a
recipOneMinus = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (a -> [a]
forall a. a -> [a]
repeat a
forall a. Semiring a => a
one)

-- | \(e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots\),
-- converges everywhere.  All derivatives are 1.
expS :: (FromRational a) => PowerSeries a
expS :: forall a. FromRational a => PowerSeries a
expS = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives (a -> [a]
forall a. a -> [a]
repeat a
forall a. Semiring a => a
one)

-- | \(\log (1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \ldots \),
-- converges for \(-1 < x \leq 1\)
logOnePlus :: (FromRational a) => PowerSeries a
logOnePlus :: forall a. FromRational a => PowerSeries a
logOnePlus = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral ([a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
1, -a
1]))

-- | \( (1 + x)^a = \sum_{n=0}^\infty \binom{a}{n} x^n \) (binomial series),
-- converges for \(|x| < 1\)
powerOnePlus :: (FromRational a) => a -> PowerSeries a
powerOnePlus :: forall a. FromRational a => a -> PowerSeries a
powerOnePlus a
a = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives ((a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Semiring a => a -> a -> a
(*) a
1 ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. AbelianGroup a => a -> a -> a
subtract a
1) a
a))

-- | \(\sin x\), converges everywhere (but faster for smaller \(x\))
sinS :: (FromRational a) => PowerSeries a
sinS :: forall a. FromRational a => PowerSeries a
sinS = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
0, a
1, a
0, -a
1])

-- | \(\cos x\), converges everywhere (but faster for smaller \(x\))
cosS :: (FromRational a) => PowerSeries a
cosS :: forall a. FromRational a => PowerSeries a
cosS = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
1, a
0, -a
1, a
0])

-- | \(\tan x\), converges for \(|x| < {\pi \over 2}\)
tanS :: (FromRational a) => PowerSeries a
tanS :: forall a. FromRational a => PowerSeries a
tanS = PowerSeries a
forall a. FromRational a => PowerSeries a
sinS PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a
forall a. FromRational a => PowerSeries a
secS

-- | \(1 \over \cos x\), converges for \(|x| < {\pi \over 2}\)
secS :: (FromRational a) => PowerSeries a
secS :: forall a. FromRational a => PowerSeries a
secS = PowerSeries a -> PowerSeries a
forall a. Ring a => PowerSeries a -> PowerSeries a
recipSimple PowerSeries a
forall a. FromRational a => PowerSeries a
cosS

-- | \(\sin^{-1} x\), converges for \(|x| \leq 1\)
asinS :: (FromRational a) => PowerSeries a
asinS :: forall a. FromRational a => PowerSeries a
asinS = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral (PowerSeries a -> PowerSeries a) -> PowerSeries a -> PowerSeries a
forall a b. (a -> b) -> a -> b
$ a -> PowerSeries a
forall a. FromRational a => a -> PowerSeries a
powerOnePlus (-a
0.5) PowerSeries a -> a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> a -> PowerSeries a
.* (-a
1) PowerSeries a -> Int -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a
.^ Int
2

-- | \(\tan^{-1} x\), converges for \(|x| \leq 1\)
atanS :: (FromRational a) => PowerSeries a
atanS :: forall a. FromRational a => PowerSeries a
atanS = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral ([a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
1, a
0, -a
1, a
0]))

-- | \(\sinh x\), converges everywhere (but faster for smaller \(x\))
sinhS :: (FromRational a) => PowerSeries a
sinhS :: forall a. FromRational a => PowerSeries a
sinhS = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
forall a. AdditiveMonoid a => a
zero, a
forall a. Semiring a => a
one])

-- | \(\cos x\), converges everywhere (but faster for smaller \(x\))
coshS :: (FromRational a) => PowerSeries a
coshS :: forall a. FromRational a => PowerSeries a
coshS = [a] -> PowerSeries a
forall a. FromRational a => [a] -> PowerSeries a
fromDerivatives ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
forall a. Semiring a => a
one, a
forall a. AdditiveMonoid a => a
zero])

-- | \(\tanh x\), converges for \(|x| < {\pi \over 2}\)
tanhS :: (FromRational a) => PowerSeries a
tanhS :: forall a. FromRational a => PowerSeries a
tanhS = PowerSeries a
forall a. FromRational a => PowerSeries a
sinhS PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a
forall a. FromRational a => PowerSeries a
sechS

-- | \(1 \over \cosh x\), converges for \(|x| < {\pi \over 2}\)
sechS :: (FromRational a) => PowerSeries a
sechS :: forall a. FromRational a => PowerSeries a
sechS = PowerSeries a -> PowerSeries a
forall a. Ring a => PowerSeries a -> PowerSeries a
recipSimple PowerSeries a
forall a. FromRational a => PowerSeries a
coshS

-- | \(\sinh^{-1} x\), converges for \(|x| \leq 1\)
asinhS :: (FromRational a) => PowerSeries a
asinhS :: forall a. FromRational a => PowerSeries a
asinhS = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral (PowerSeries a -> PowerSeries a) -> PowerSeries a -> PowerSeries a
forall a b. (a -> b) -> a -> b
$ a -> PowerSeries a
forall a. FromRational a => a -> PowerSeries a
powerOnePlus (-a
0.5) PowerSeries a -> Int -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a
.^ Int
2

-- | \(\tanh^{-1} x\), converges for \(|x| < 1\)
atanhS :: (FromRational a) => PowerSeries a
atanhS :: forall a. FromRational a => PowerSeries a
atanhS = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral ([a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
cycle [a
forall a. Semiring a => a
one, a
forall a. AdditiveMonoid a => a
zero]))

-- | Euler function (<https://oeis.org/A010815 OEIS A010815>), defined as
-- \[
--	\phi(x) = \prod_{k=1}^\infty (1-x^k)
--	= \sum_{k = -\infty}^\infty (-1)^k x^{k(3k+1) \over 2}
-- \]
euler :: (Ring a) => PowerSeries a
euler :: forall a. Ring a => PowerSeries a
euler =
    [(a, Int)] -> PowerSeries a
forall a. AdditiveMonoid a => [(a, Int)] -> PowerSeries a
expand [(Int -> a
forall a. Ring a => Int -> a
sign Int
k, Int
kInt -> Int -> Int
forall a. Semiring a => a -> a -> a
*(Int
3Int -> Int -> Int
forall a. Semiring a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. AdditiveMonoid a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Euclidean a => a -> a -> a
`div` Int
2) | Int
k <- [Int
0..]] PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => a -> a -> a
+
    [(a, Int)] -> PowerSeries a
forall a. AdditiveMonoid a => [(a, Int)] -> PowerSeries a
expand [(Int -> a
forall a. Ring a => Int -> a
sign Int
k, Int
kInt -> Int -> Int
forall a. Semiring a => a -> a -> a
*(Int
3Int -> Int -> Int
forall a. Semiring a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. AbelianGroup a => a -> a -> a
-Int
1) Int -> Int -> Int
forall a. Euclidean a => a -> a -> a
`div` Int
2) | Int
k <- [Int
1..]]

-- expand a list of coefficients and indices (in ascending order of index)
expand :: (AdditiveMonoid a) => [(a, Int)] -> PowerSeries a
expand :: forall a. AdditiveMonoid a => [(a, Int)] -> PowerSeries a
expand = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients ([a] -> PowerSeries a)
-> ([(a, Int)] -> [a]) -> [(a, Int)] -> PowerSeries a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [(a, Int)] -> [a]
forall {a}. AdditiveMonoid a => Int -> [(a, Int)] -> [a]
aux Int
0
  where
    aux :: Int -> [(a, Int)] -> [a]
aux Int
_ [] = []
    aux Int
n ((a
a, Int
i):[(a, Int)]
ais) = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
iInt -> Int -> Int
forall a. AbelianGroup a => a -> a -> a
-Int
n) a
forall a. AdditiveMonoid a => a
zero [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [(a, Int)] -> [a]
aux (Int
iInt -> Int -> Int
forall a. AdditiveMonoid a => a -> a -> a
+Int
1) [(a, Int)]
ais

-- (-1)^k
sign :: (Ring a) => Int -> a
sign :: forall a. Ring a => Int -> a
sign Int
k
  | Int -> Bool
forall a. ToInteger a => a -> Bool
odd Int
k = -a
1
  | Bool
otherwise = a
1

-- | The infinite list of evaluations of truncations of the power series.
approximations :: (Semiring a) => PowerSeries a -> a -> [a]
approximations :: forall a. Semiring a => PowerSeries a -> a -> [a]
approximations = (a -> a -> a) -> PowerSeries a -> a -> [a]
forall a b c.
(AdditiveMonoid a, Semiring b, AdditiveMonoid c) =>
(a -> b -> c) -> PowerSeries a -> b -> [c]
approximationsWith a -> a -> a
forall a. Semiring a => a -> a -> a
(*)

-- | The infinite list of evaluations of truncations of the power series,
-- given a function to combine coefficients and powers.
approximationsWith :: (AdditiveMonoid a, Semiring b, AdditiveMonoid c) =>
    (a -> b -> c) -> PowerSeries a -> b -> [c]
approximationsWith :: forall a b c.
(AdditiveMonoid a, Semiring b, AdditiveMonoid c) =>
(a -> b -> c) -> PowerSeries a -> b -> [c]
approximationsWith a -> b -> c
f PowerSeries a
as b
x =
    (c -> c -> c) -> c -> [c] -> [c]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl c -> c -> c
forall a. AdditiveMonoid a => a -> a -> a
(+) c
forall a. AdditiveMonoid a => a
zero ((a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> b -> c
f (PowerSeries a -> [a]
forall a. AdditiveMonoid a => PowerSeries a -> [a]
coefficients PowerSeries a
as) ((b -> b) -> b -> [b]
forall a. (a -> a) -> a -> [a]
iterate (b -> b -> b
forall a. Semiring a => a -> a -> a
*b
x) b
forall a. Semiring a => a
one))

-- | The list of Padé approximants of \(A(x)\) of order \([m/n]\) such
-- that \(k = m+n\) in order of increasing \(n\).
-- Each approximant is a ratio of polynomials \(P(x)\) and \(Q(x)\)
-- of degrees \(m\) and \(n\) respectively, such that
--
-- \[
-- A(x) \equiv {P(x) \over Q(x)} \pmod{x^{m+n+1}}
-- \]
--
-- The list is obtained by applying the extended Euclidean algorithm
-- to the polynomials \(x^{k+1}\) and \(A_k(x) = A(x) \mod x^{k+1}\)
-- (i.e. selecting the terms up to \(x^k\)), which yields a sequence
-- of equations
--
-- \[
-- P_i(x) = S_i(x) x^{k+1} + Q_i(x) A_k(x)
-- \]
--
padeApproximants :: (Eq a, Field a) =>
    Int -> PowerSeries a -> [RationalFunction a]
padeApproximants :: forall a.
(Eq a, Field a) =>
Int -> PowerSeries a -> [RationalFunction a]
padeApproximants Int
k PowerSeries a
f = [Polynomial a
pPolynomial a -> Polynomial a -> RationalFunction a
forall a.
(Eq a, StandardAssociate a, Euclidean a) =>
a -> a -> Ratio a
%Polynomial a
q | (Polynomial a
_, Polynomial a
p, Polynomial a
_, Polynomial a
q) <- Polynomial a
-> Polynomial a
-> [(Polynomial a, Polynomial a, Polynomial a, Polynomial a)]
forall a.
(Eq a, AbelianGroup a, StandardAssociate a, Euclidean a) =>
a -> a -> [(a, a, a, a)]
extendedEuclid Polynomial a
xk1 Polynomial a
truncation]
  where
    truncation :: Polynomial a
truncation = PowerSeries a -> PowerSeries a -> Polynomial a
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Polynomial a
modulus PowerSeries a
f (Polynomial a -> PowerSeries a
forall a. AdditiveMonoid a => Polynomial a -> PowerSeries a
fromPolynomial Polynomial a
xk1)
    xk1 :: Polynomial a
xk1 = Polynomial a
forall a. Semiring a => Polynomial a
Poly.identityPolynomial a -> Int -> Polynomial a
forall a b. (Semiring a, ToInteger b) => a -> b -> a
^(Int
kInt -> Int -> Int
forall a. AdditiveMonoid a => a -> a -> a
+Int
1)

-- | Degree of the first non-zero coefficient.
-- Undefined if the argument is zero.
order :: (Eq a, AdditiveMonoid a) => PowerSeries a -> Int
order :: forall a. (Eq a, AdditiveMonoid a) => PowerSeries a -> Int
order PowerSeries a
f = PowerSeries a -> PowerSeries a -> Int
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Int
countSame PowerSeries a
f PowerSeries a
forall a. AdditiveMonoid a => a
zero

-- | Metric on power series (undefined if the arguments are equal)
distance :: (Eq a, AdditiveMonoid a) =>
    PowerSeries a -> PowerSeries a -> Rational
distance :: forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Rational
distance PowerSeries a
f PowerSeries a
g = Rational
2Rational -> Int -> Rational
forall a b.
(DivisionSemiring a, AbelianGroup b, ToInteger b) =>
a -> b -> a
^^(Int -> Int
forall a. AbelianGroup a => a -> a
negate (PowerSeries a -> PowerSeries a -> Int
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Int
countSame PowerSeries a
f PowerSeries a
g))

countSame :: (Eq a, AdditiveMonoid a) => PowerSeries a -> PowerSeries a -> Int
countSame :: forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Int
countSame PowerSeries a
f PowerSeries a
g
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
f a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g = PowerSeries a -> PowerSeries a -> Int
forall a.
(Eq a, AdditiveMonoid a) =>
PowerSeries a -> PowerSeries a -> Int
countSame (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g) Int -> Int -> Int
forall a. AdditiveMonoid a => a -> a -> a
+ Int
1
  | Bool
otherwise = Int
0

from :: (Semiring a) => a -> [a]
from :: forall a. Semiring a => a -> [a]
from = (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. AdditiveMonoid a => a -> a -> a
+a
forall a. Semiring a => a
one)

-- | Identity function, i.e. the indeterminate
identity :: (Semiring a) => PowerSeries a
identity :: forall a. Semiring a => PowerSeries a
identity = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
forall a. Semiring a => a
one

-- | composition \(f(g(x))\). This is only defined if \(g_0\) is 0.
compose :: (Eq a, Semiring a) =>
    PowerSeries a -> PowerSeries a -> PowerSeries a
compose :: forall a.
(Eq a, Semiring a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
compose PowerSeries a
f PowerSeries a
g
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
g a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
forall a. AdditiveMonoid a => a
zero =
    [Char] -> PowerSeries a
forall a. HasCallStack => [Char] -> a
error [Char]
"compose: second series has non-zero constant term"
  | Bool
otherwise = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX PowerSeries a
f (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g)

-- | Special composition \(f(x g(x))\).
composeMulX :: (Semiring a) =>
    PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX :: forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX PowerSeries a
f PowerSeries a
g = a -> PowerSeries a -> PowerSeries a
forall a. a -> PowerSeries a -> PowerSeries a
plusMulX (PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
f) (PowerSeries a
g PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) PowerSeries a
g)

-- | pre-inverse under 'compose':
--
-- * @'compose' f ('inverse' f) = 'identity'@
--
-- This is only defined if \(f_0\) is zero and \(f_1\) is non-zero.
inverse :: (Eq a, Field a) => PowerSeries a -> PowerSeries a
inverse :: forall a. (Eq a, Field a) => PowerSeries a -> PowerSeries a
inverse PowerSeries a
f
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero PowerSeries a
f a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0 = [Char] -> PowerSeries a
forall a. HasCallStack => [Char] -> a
error [Char]
"inverse: zero constant term"
  | PowerSeries a -> a
forall a. AdditiveMonoid a => PowerSeries a -> a
atZero (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = [Char] -> PowerSeries a
forall a. HasCallStack => [Char] -> a
error [Char]
"inverse: non-zero linear term"
  | Bool
otherwise = PowerSeries a -> PowerSeries a
forall a. Field a => PowerSeries a -> PowerSeries a
unsafeInverse PowerSeries a
f

-- assumes \(f_0\) is zero and \(f_1\) is non-zero.
unsafeInverse :: (Field a) => PowerSeries a -> PowerSeries a
unsafeInverse :: forall a. Field a => PowerSeries a -> PowerSeries a
unsafeInverse PowerSeries a
f = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
r
  where
    r :: PowerSeries a
r = PowerSeries a
1 PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Field a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
`exactDiv` PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
f) PowerSeries a
r

-- | Special case of 'inverse' that is only valid if the constant term
-- is 0 and the coefficient of the linear term is 1.
inverseSimple :: (Eq a, Ring a) => PowerSeries a -> PowerSeries a
inverseSimple :: forall a. (Eq a, Ring a) => PowerSeries a -> PowerSeries a
inverseSimple PowerSeries a
p = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
r
  where
    r :: PowerSeries a
r = PowerSeries a -> PowerSeries a
forall a. Ring a => PowerSeries a -> PowerSeries a
recipSimple (PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
p) PowerSeries a
r)

-- | Maps a function \(f(x)\) to \(f(a x)\), equivalent to composition
-- with @a*'identity'@.
(.*) :: (Semiring a) => PowerSeries a -> a -> PowerSeries a
PS [a]
as .* :: forall a. Semiring a => PowerSeries a -> a -> PowerSeries a
.* a
a = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Semiring a => a -> a -> a
(*) [a]
as ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
*a
a) a
forall a. Semiring a => a
one))

-- | Maps a function \(f(x)\) to \(f(x^k)\) for positive \(k\), equivalent
-- to composition with @'identity'^k@.
(.^) :: (AdditiveMonoid a) => PowerSeries a -> Int -> PowerSeries a
PS [a]
as .^ :: forall a. AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a
.^ Int
k
  | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = [Char] -> PowerSeries a
forall a. HasCallStack => [Char] -> a
error [Char]
"non-positive power"
  | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS [a]
as
  | Bool
otherwise = [a] -> PowerSeries a
forall a. [a] -> PowerSeries a
PS ([a] -> [[a]] -> [a]
forall a. [a] -> [[a]] -> [a]
intercalate (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
kInt -> Int -> Int
forall a. AbelianGroup a => a -> a -> a
-Int
1) a
forall a. AdditiveMonoid a => a
zero) ((a -> [a]) -> [a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map a -> [a]
forall a. a -> [a]
forall (f :: * -> *) a. Applicative f => a -> f a
pure [a]
as))

-- | Special case of 'recip' that is only valid if the constant term is 1.
recipSimple :: (Ring a) => PowerSeries a -> PowerSeries a
-- recipSimple p = compose recipOneMinus (1 - p)
recipSimple :: forall a. Ring a => PowerSeries a -> PowerSeries a
recipSimple PowerSeries a
p = PowerSeries a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> PowerSeries a
recipOneMinusOf (PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a
negate (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
p))

-- | @recipOneMinusOf p = 1/(1 - 'identity'*p)) = 1 + p + p^2 + ...@,
-- but with a less constrained type.
recipOneMinusOf :: (Semiring a) => PowerSeries a -> PowerSeries a
-- recipOneMinusOf p = compose recipOneMinus (mulX p)
recipOneMinusOf :: forall a. Semiring a => PowerSeries a -> PowerSeries a
recipOneMinusOf PowerSeries a
p = PowerSeries a
q
  where
    q :: PowerSeries a
q = a -> PowerSeries a -> PowerSeries a
forall a. a -> PowerSeries a -> PowerSeries a
plusMulX a
forall a. Semiring a => a
one (PowerSeries a
pPowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
*PowerSeries a
q)

-- | Special case of division that is only valid if the constant term
-- of the divisor is 1.
divSimple :: (Eq a, Ring a) => PowerSeries a -> PowerSeries a -> PowerSeries a
divSimple :: forall a.
(Eq a, Ring a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divSimple PowerSeries a
p PowerSeries a
q = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus PowerSeries a
p (PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a
negate (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
q))

-- | @divOneMinus p q = p/(1 - 'identity'*q) = p*(1 + q + q^2 + ...)@,
-- but with a less constrained type.
divOneMinus :: (Semiring a) => PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus :: forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus PowerSeries a
as PowerSeries a
bs = PowerSeries a
q
  where
    q :: PowerSeries a
q = PowerSeries a
as PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => a -> a -> a
+ PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX (PowerSeries a
bsPowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
*PowerSeries a
q)

{- | A Stieltjes fraction
(<https://dlmf.nist.gov/3.10#Px1 Digital Library of Mathematical Functions 3.10>)
generated by a sequence \(a_n\) has the form

\[
    \cfrac{1}{1 - \cfrac{a_0 x}{1 - \cfrac{a_1 x}{1 - \cdots}}}
\]

=== __Examples__

* The Catalan numbers \( C(n) = {(2n)! \over n! (n+1)!} \)
(<https://oeis.org/A000108 OEIS A000108>):

>>> coefficients $ stieltjes $ repeat 1
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440,...

* As the coefficients are repeated, the above is equivalent to a recursive
definition \( C = {1 \over 1-xC} \):

>>> let catalan = recipOneMinus `compose` mulX catalan
>>> coefficients catalan
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440,...

* Double factorials of odd numbers
(<https://oeis.org/A001147 OEIS A001147>):

\[ (2n-1)!! = { (2n)! \over 2^n n! } = 1.3.5. \ldots .(2n-1) \]

>>> coefficients $ stieltjes [1..]
[1,1,3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575,...

* /Bell numbers/, the number of ways to partition \(n\) elements
(<https://oeis.org/A000110 OEIS A000110>):

>>> coefficients $ stieltjes $ concat [[1, n] | n <- [1..]]
[1,1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,...

-}
stieltjes :: (Semiring a) => [a] -> PowerSeries a
stieltjes :: forall a. Semiring a => [a] -> PowerSeries a
stieltjes [a]
as = PowerSeries a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> PowerSeries a
recipOneMinusOf ((PowerSeries a -> PowerSeries a -> PowerSeries a)
-> PowerSeries a -> [PowerSeries a] -> PowerSeries a
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus PowerSeries a
forall a. AdditiveMonoid a => a
zero ((a -> PowerSeries a) -> [a] -> [PowerSeries a]
forall a b. (a -> b) -> [a] -> [b]
map a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant [a]
as))

{- | Deléham's delta operator (<https://oeis.org/A084938 OEIS A084938>)
takes two sequences \(\{r_n\}\) and \(\{s_n\}\) and forms the Stieltjes
fraction defined by the sequence of polynomials \(\{r_n + s_n x\}\),
yielding the ordinary generating function of a polynomial sequence.

=== __Examples__

* The Catalan or Narayana triangle
(<https://oeis.org/A090181 OEIS A090181>):

>>> coefficients $ delta (cycle [0, 1]) (cycle [1, 0])
[fromCoefficients [1],
 fromCoefficients [0,1],
 fromCoefficients [0,1,1],
 fromCoefficients [0,1,3,1],
 fromCoefficients [0,1,6,6,1],
 fromCoefficients [0,1,10,20,10,1],...

* The sequence of rising factorial polynomials with coefficients the
unsigned Stirling numbers of the first kind
(<https://oeis.org/A132393 OEIS A132393>), the number of permutations
of \(n\) elements that have exactly \(k\) cycles.

>>> coefficients $ delta (0:concat [[n, n] | n <- [1..]]) (cycle [1, 0])
[fromCoefficients [1],
 fromCoefficients [0,1],
 fromCoefficients [0,1,1],
 fromCoefficients [0,2,3,1],
 fromCoefficients [0,6,11,6,1],
 fromCoefficients [0,24,50,35,10,1],...

* The sequence of Touchard polynomials with coefficients the Stirling
numbers of the second kind (<https://oeis.org/A048993 OEIS A048993>),
the number of ways to partition a set of \(n\) elements into \(k\)
non-empty subsets:

>>> coefficients $ delta (concat [[0, n] | n <- [1..]]) (cycle [1, 0])
[fromCoefficients [1],
 fromCoefficients [0,1],
 fromCoefficients [0,1,1],
 fromCoefficients [0,1,3,1],
 fromCoefficients [0,1,7,6,1],
 fromCoefficients [0,1,15,25,10,1],...

* The sequence of Eulerian polynomials with coefficients the Eulerian
numbers \(A(n,k)\) (<https://oeis.org/A123125 OEIS A123125>), the number
of permutations of the numbers 1 to \(n\) in which exactly \(k\) elements
are greater than the previous element.

>>> coefficients $ delta (concat [[n, 0] | n <- [1..]]) (concat [[0, n] | n <- [1..]])
[fromCoefficients [1],
 fromCoefficients [1],
 fromCoefficients [1,1],
 fromCoefficients [1,4,1],
 fromCoefficients [1,11,11,1],
 fromCoefficients [1,26,66,26,1],...

-}
delta :: (Semiring a) => [a] -> [a] -> PowerSeries (Polynomial a)
delta :: forall a. Semiring a => [a] -> [a] -> PowerSeries (Polynomial a)
delta [a]
r [a]
s =
    [Polynomial a] -> PowerSeries (Polynomial a)
forall a. Semiring a => [a] -> PowerSeries a
stieltjes ((Polynomial a -> Polynomial a -> Polynomial a)
-> [Polynomial a] -> [Polynomial a] -> [Polynomial a]
forall a. (a -> a -> a) -> [a] -> [a] -> [a]
longZipWith Polynomial a -> Polynomial a -> Polynomial a
forall a. AdditiveMonoid a => a -> a -> a
(+) ((a -> Polynomial a) -> [a] -> [Polynomial a]
forall a b. (a -> b) -> [a] -> [b]
map a -> Polynomial a
forall a. a -> Polynomial a
Poly.constant [a]
r) ((a -> Polynomial a) -> [a] -> [Polynomial a]
forall a b. (a -> b) -> [a] -> [b]
map a -> Polynomial a
forall a. AdditiveMonoid a => a -> Polynomial a
Poly.linear [a]
s))

-- A085880 = delta (repeat one) (repeat one)

{- | A Jacobi fraction
(<https://dlmf.nist.gov/3.10#Px3 Digital Library of Mathematical Functions 3.10>)
generated by sequences \(a_n\) and \(b_n\) has the form

\[
    \cfrac{1}{1 - a_0 x - \cfrac{b_0 x^2}{1 - a_1 x - \cfrac{b_1 x^2}{1 - a_2 x - \cdots}}}
\]

The even part of a Stieltjes fraction is an equivalent Jacobi fraction:

\[
    \cfrac{1}{1 - \cfrac{a_0 x}{1 - \cfrac{a_1 x}{1 - \cfrac{a_2 x}{1 - \cfrac{a_3 x}{1 - \cfrac{a_4 x}{1 - \cdots}}}}}}
    =
    \cfrac{1}{1 - a_0 x - \cfrac{a_0 a_1 x^2}{1 - (a_1+a_2) x - \cfrac{a_2 a_3 x^2}{1 - (a_3+a_4) x - \cdots}}}
\]

=== __Examples__

* Using the above identity, yet another expression for the
/Catalan numbers/ (<https://oeis.org/A000108 OEIS A000108>) is:

>>> coefficients $ jacobi (1:repeat 2) (repeat 1)
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440,...

* Another instance of this identity yields a Jacobi fraction for the
/Bell numbers/ (<https://oeis.org/A000110 OEIS A000110>):

>>> coefficients $ jacobi [1..] [1..]
[1,1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,...

* /Motzkin numbers/, the number of ways to draw non-intersecting chords
joining \(n\) points on a circle
(<https://oeis.org/A001006 OEIS A001006>):

>>> coefficients $ jacobi (repeat 1) (repeat 1)
[1,1,2,4,9,21,51,127,323,835,2188,5798,15511,41835,113634,310572,...

* As the coefficients are repeated, the above is equivalent to a recursive
definition \( M = {1 \over 1-x(1+xM)} \):

>>> let motzkin = recipOneMinus `compose` mulX (1 + mulX motzkin)
>>> coefficients motzkin
[1,1,2,4,9,21,51,127,323,835,2188,5798,15511,41835,113634,310572,...

* The number of involutions on a set of \(n\) elements
(<https://oeis.org/A000085 OEIS A000085>):

>>> coefficients $ jacobi (repeat 1) [1..]
[1,1,2,4,10,26,76,232,764,2620,9496,35696,140152,568504,...

* /Bessel numbers/, the number of non-overlapping partitions of \(n\)
elements into equivalence classes (<https://oeis.org/A006789 OEIS A006789>):

>>> coefficients $ jacobi [1..] (repeat 1)
[1,1,2,5,14,43,143,509,1922,7651,31965,139685,636712,3020203,...

-}
jacobi :: (Semiring a) => [a] -> [a] -> PowerSeries a
jacobi :: forall a. Semiring a => [a] -> [a] -> PowerSeries a
jacobi [a]
as [a]
bs = ((a, a) -> PowerSeries a -> PowerSeries a)
-> PowerSeries a -> [(a, a)] -> PowerSeries a
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (a, a) -> PowerSeries a -> PowerSeries a
forall {a}. Semiring a => (a, a) -> PowerSeries a -> PowerSeries a
fraction PowerSeries a
forall a. AdditiveMonoid a => a
zero ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([a]
as [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
forall a. AdditiveMonoid a => a
zero) (a
forall a. Semiring a => a
onea -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
bs))
  where
    fraction :: (a, a) -> PowerSeries a -> PowerSeries a
fraction (a
a, a
b) PowerSeries a
p = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus (a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant a
b) (a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant a
a PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => a -> a -> a
+ PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
p)

{- |

If \(A(x)\) is the generating function of the
sequence \( \{a_n\} \), then \( \frac{1}{1-x} A\left( \frac{x}{1-x} \right) \)
is the generating function of the sequence \( \{b_n\} \) where

\[
    b_n = \sum_{k=0}^{n} \binom{n}{k} a_k
\]

For the same transform on sequences from exponential generating functions,
see "Data.YAP.PowerSeries.Maclaurin#g:binomial_transform".

=== __Examples__

* Maximal number of pieces obtained by cutting a circle with \(n\) lines
(<https://oeis.org/A000124 OEIS A000124>).

>>> coefficients $ binomialTransform $ fromCoefficients [1,1,1]
[1,2,4,7,11,16,22,29,37,46,56,67,79,92,106,121,137,154,172,191,211,232,...

* Maximal number of pieces obtained by cutting a sphere with \(n\) planes
(<https://oeis.org/A000125 OEIS A000125>).

>>> coefficients $ binomialTransform $ fromCoefficients [1,1,1,1]
[1,2,4,8,15,26,42,64,93,130,176,232,299,378,470,576,697,834,988,1160,...

* Maximal number of regions obtained by the straight lines joining \(n\)
points around a circle (<https://oeis.org/A000127 OEIS A000127>).

>>> coefficients $ binomialTransform $ fromCoefficients [1,1,1,1,1]
[1,2,4,8,16,31,57,99,163,256,386,562,794,1093,1471,1941,2517,3214,...

-}
binomialTransform :: (Semiring a) => PowerSeries a -> PowerSeries a
binomialTransform :: forall a. Semiring a => PowerSeries a -> PowerSeries a
binomialTransform PowerSeries a
a = PowerSeries a
forall a. Semiring a => PowerSeries a
recipOneMinus PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX PowerSeries a
a PowerSeries a
forall a. Semiring a => PowerSeries a
recipOneMinus

{- |
The inverse of `binomialTransform`, namely

\[
    a_n = \sum_{k=0}^{n} (-1)^{n-k} \binom{n}{k} b_k
\]

is defined by

\[
	A(x) = \frac{1}{1+x} B\left( {x \over 1+x} \right)
\]

For the same transform on sequences from exponential generating functions,
see "Data.YAP.PowerSeries.Maclaurin#g:binomial_transform".

-}
inverseBinomialTransform :: (Ring a) => PowerSeries a -> PowerSeries a
inverseBinomialTransform :: forall a. Ring a => PowerSeries a -> PowerSeries a
inverseBinomialTransform PowerSeries a
a = PowerSeries a
recipOnePlus PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
composeMulX PowerSeries a
a PowerSeries a
recipOnePlus
  where
    recipOnePlus :: PowerSeries a
recipOnePlus = PowerSeries a
forall a. Semiring a => PowerSeries a
recipOneMinus PowerSeries a -> a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> a -> PowerSeries a
.* (-a
1)

{- |
This transform maps a power series \( A(x) = \sum_{n=1}^\infty a_n x^n\)
(with a zero constant term)
to the Lambert series with the same coefficients:

\[
	B(x)
	= \sum_{n=1}^\infty A(x^n)
	= \sum_{n=1}^\infty a_n \frac{x^n}{1-x^n}
	= \sum_{n=1}^\infty b_n x^n
\]

where \(b_n = \sum_{k \mid n} a_k\).
For the same sum-of-divisors transform on
sequences from Dirichlet generating functions, see
"Data.YAP.DirichletSeries#g:sum_of_divisors_transform".

=== __Examples__

* Number of divisors of \(n\) is generated by the Lambert transform
of \(x \over 1-x\):

>>> coefficients $ lambertTransform $ mulX recipOneMinus
[0,1,2,2,3,2,4,2,4,3,4,2,6,2,4,4,5,2,6,2,...

* Sum of divisors of \(n\) is generated by the Lambert transform
of \(x \over (1-x)^2\):

>>> coefficients $ lambertTransform $ mulX $ derivative recipOneMinus
[0,1,3,4,7,6,12,8,15,13,18,12,28,14,24,24,31,18,39,20,...

-}
lambertTransform :: (AdditiveMonoid a) => PowerSeries a -> PowerSeries a
lambertTransform :: forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
lambertTransform PowerSeries a
p = PowerSeries (PowerSeries a) -> PowerSeries a
forall a.
AdditiveMonoid a =>
PowerSeries (PowerSeries a) -> PowerSeries a
flatten (PowerSeries (PowerSeries a) -> PowerSeries a)
-> PowerSeries (PowerSeries a) -> PowerSeries a
forall a b. (a -> b) -> a -> b
$ PowerSeries (PowerSeries a) -> PowerSeries (PowerSeries a)
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX (PowerSeries (PowerSeries a) -> PowerSeries (PowerSeries a))
-> PowerSeries (PowerSeries a) -> PowerSeries (PowerSeries a)
forall a b. (a -> b) -> a -> b
$ [PowerSeries a] -> PowerSeries (PowerSeries a)
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients [PowerSeries a
p' PowerSeries a -> Int -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a
.^ Int
k | Int
k <- [Int
1..]]
  where
    p' :: PowerSeries a
p' = PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
p

{- |
The inverse of 'lambertTransform' computes the Möbius transform of the
generated sequence, and can be computed from

\[
	A(x) = B(x) - \sum_{n=2}^\infty A(x^n)
\]

provided \(B(x)\) has a zero constant term.
For the same transform on sequences from Dirichlet generating functions,
see "Data.YAP.DirichletSeries#g:moebius_transform".

=== __Examples__

* The Möbius function \( \mu(n) \) (<https://oeis.org/A008683 OEIS A008683>,
see "Data.YAP.DirichletSeries") has the ordinary generating function whose
Lambert transform is \(x\):

>>> coefficients $ inverseLambertTransform identity
[0,1,-1,-1,0,-1,1,-1,0,0,1,-1,0,-1,1,1,0,-1,0,-1,0,1,1,-1,0,...

* The number of numbers less than \(n\) and coprime with \(n\)
(<https://oeis.org/A000010 OEIS A000010>) is the Euler totient
function \( \varphi(n) \), which has the ordinary generating function
whose Lambert transform is \(x \over 1-x\):

>>> coefficients $ inverseLambertTransform $ mulX $ derivative recipOneMinus
[0,1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,...

-}
inverseLambertTransform :: (Ring a) => PowerSeries a -> PowerSeries a
inverseLambertTransform :: forall a. Ring a => PowerSeries a -> PowerSeries a
inverseLambertTransform PowerSeries a
q = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX PowerSeries a
p'
  where
    p' :: PowerSeries a
p' = PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
q PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a -> a
- PowerSeries (PowerSeries a) -> PowerSeries a
forall a.
AdditiveMonoid a =>
PowerSeries (PowerSeries a) -> PowerSeries a
flatten ([PowerSeries a] -> PowerSeries (PowerSeries a)
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (PowerSeries a
0 PowerSeries a -> [PowerSeries a] -> [PowerSeries a]
forall a. a -> [a] -> [a]
: [PowerSeries a
p' PowerSeries a -> Int -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a
.^ Int
k | Int
k <- [Int
2..]]))

{- |
A power series \(B(x)\) with \(b_0 = 1\) is the Euler transform of
a series \(A(x)\) with \(a_0 = 0\) if

\[
	B(x)
	= \prod_{n=1}^\infty {1 \over (1 - x^i)^{a_i}}
	= \exp \left( \sum_{n=1}^\infty {A(x^n) \over n} \right)
\]

This can be expressed using a Lambert transform to an intermediate sequence:

\[
	C(x)
	= x {d \over dx} \log(B(x))
	= \sum_{n=1}^\infty x^n (A^\prime(x^n))
\]

The product form implies that if the coefficients of \(A(x)\) are
non-negative integers, so are the coefficients of \(B(x)\), but it doesn't
seem possible to prove this to the type system.

Let \(P\) be a property of unlabelled graphs that holds of a graph if and
only if it holds for each of the connected components of that graph.
Then if \(a_n\) is the number of connected graphs on \(n\) nodes with
property \(P\), then \(b_n\) is the number of all graphs on \(n\) nodes
with property \(P\).

=== __Examples__

Since \(n\)-node rooted unlabelled trees are in one-to-one correspondence
with \(n-1\)-node forests of rooted unlabelled trees, we can define
the ordinary generating function for the number of \(n\)-node rooted
unlabelled trees (<https://oeis.org/A000081 OEIS A000081>):

>>> let rootedTrees = mulX (eulerTransform rootedTrees)
>>> map Data.Ratio.numerator $ coefficients rootedTrees
[0,1,1,2,4,9,20,48,115,286,719,1842,4766,12486,32973,87811,235381,...

Counting unrooted trees (<https://oeis.org/A000055 OEIS A000055>) is
more obscure:

>>> let trees = 1 + rootedTrees + 0.5*(rootedTrees .^ 2 - rootedTrees^2)
>>> map Data.Ratio.numerator $ coefficients trees
[1,1,1,1,2,3,6,11,23,47,106,235,551,1301,3159,7741,19320,48629,123867,...

Then the ordinary generating function for the number of acyclic undirected
graphs on \(n\) nodes (<https://oeis.org/A005195 OEIS A005195>) is the
Euler transform of @trees@:

>>> map Data.Ratio.numerator $ coefficients $ eulerTransform $ trees-1
[1,1,2,3,6,10,20,37,76,153,329,710,1601,3658,8599,20514,49905,122963,...

-}
eulerTransform :: (Eq a, FromRational a) => PowerSeries a -> PowerSeries a
eulerTransform :: forall a. (Eq a, FromRational a) => PowerSeries a -> PowerSeries a
eulerTransform PowerSeries a
a =
    PowerSeries a
forall a. FromRational a => PowerSeries a
expS PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
(Eq a, Semiring a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
`compose` PowerSeries a -> PowerSeries a
forall a. FromRational a => PowerSeries a -> PowerSeries a
divN (PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
lambertTransform (PowerSeries a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> PowerSeries a
mulN PowerSeries a
a))

{- |
The inverse of 'eulerTransform'.
-}
inverseEulerTransform ::
    (Eq a, FromRational a) => PowerSeries a -> PowerSeries a
inverseEulerTransform :: forall a. (Eq a, FromRational a) => PowerSeries a -> PowerSeries a
inverseEulerTransform PowerSeries a
b =
    PowerSeries a -> PowerSeries a
forall a. FromRational a => PowerSeries a -> PowerSeries a
divN (PowerSeries a -> PowerSeries a
forall a. Ring a => PowerSeries a -> PowerSeries a
inverseLambertTransform (PowerSeries a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> PowerSeries a
mulN (PowerSeries a
forall a. FromRational a => PowerSeries a
logOnePlus PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
(Eq a, Semiring a) =>
PowerSeries a -> PowerSeries a -> PowerSeries a
`compose` (PowerSeries a
bPowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a -> a
-PowerSeries a
1))))

-- | Multiply the \(n\)th term by \(n\).
mulN :: (Semiring a) => PowerSeries a -> PowerSeries a
mulN :: forall a. Semiring a => PowerSeries a -> PowerSeries a
mulN PowerSeries a
p = PowerSeries a -> PowerSeries a
forall a. AdditiveMonoid a => PowerSeries a -> PowerSeries a
mulX (PowerSeries a -> PowerSeries a
forall a. Differentiable a => a -> a
derivative PowerSeries a
p)

-- | Divide the \(n\)th term by \(n\).
divN :: (FromRational a) => PowerSeries a -> PowerSeries a
divN :: forall a. FromRational a => PowerSeries a -> PowerSeries a
divN PowerSeries a
p = PowerSeries a -> PowerSeries a
forall a. Integrable a => a -> a
integral (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
p)

{- $ogfs

A power series \( \sum_{k=0} a_k x^k \) may be viewed as an ordinary
generating function for the sequence \( \{ a_k \} \).  This sequence
may be extracted using 'coefficients'.

=== Examples

* The generating function of \( \{ {1 \over k} \} \) is

>>> coefficients $ integral recipOneMinus
[0 % 1,1 % 1,1 % 2,1 % 3,1 % 4,1 % 5,1 % 6,1 % 7,1 % 8,...

* Multiplying a generating function by \( {1 \over 1-x} \) generates a
sequence of partial sums.
An example is the Harmonic numbers \(H_n = \sum_{k=1}^n {1 \over k}\)
(<https://oeis.org/A001008 OEIS A001008> over
<https://oeis.org/A002805 OEIS A002805>):

>>> coefficients $ recipOneMinus * integral recipOneMinus
[0 % 1,1 % 1,3 % 2,11 % 6,25 % 12,137 % 60,49 % 20,363 % 140,...

* The Mertens function (<https://oeis.org/A002321 OEIS A002321>) consists of
partial sums of the Möbius function:

>>> coefficients $ recipOneMinus * inverseLambertTransform identity
[0,1,0,-1,-1,-2,-1,-2,-2,-2,-1,-2,-2,-3,-2,-1,-1,-2,-2,-3,...

* The Catalan numbers (<https://oeis.org/A000108 OEIS A000108>)
are generated by the solution of the equation \( C = 1 + x C^2 \),
namely the function \( 1 - \sqrt{1-4x} \over 2x \):

>>> let catalan = one + mulX (catalan*catalan)
>>> coefficients catalan
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440,...

* Yet another way to construct the generating function of the Catalan
numbers is to observe that the above equation is equivalent to requiring
that \(C-1\) is a solution to \( x = { y \over {(1+y)}^2 } \), which
can be obtained by inversion:

>>> coefficients $ 1 + inverseSimple (mulX (derivative recipOneMinus .* (-1)))
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440,...

* The derivative of \( 1 - \sqrt{1-4x} \over 2 \)
is \( 1 \over \sqrt{1 - 4x} \), whose coefficients are the
central binomial coefficients \( \binom{2n}{n} = {(2n!) \over {n!}^2} \)
(<https://oeis.org/A000984 OEIS A000984>):

>>> coefficients $ derivative $ mulX catalan
[1,2,6,20,70,252,924,3432,12870,48620,184756,705432,2704156,...

* The number of ways to partition \(n\) unlabelled elements
(<https://oeis.org/A000041 OEIS A000041>) is generated by the reciprocal
of 'euler'.
(This is the number of terms in the \(n\)th Bell polynomial:
see "Data.YAP.FiniteMap#g:Bell_polynomials".)

>>> coefficients $ recipSimple euler
[1,1,2,3,5,7,11,15,22,30,42,56,77,101,135,176,231,297,385,490,...

* The number of ways to partition \(n\) unlabelled elements into parts
of distinct sizes or equivalently the number of partitions into odd parts
(<https://oeis.org/A000009 OEIS A000009>) is generated by

\[
    \prod_{k=1}^{\infty} \frac{1}{1 - x^{2k-1}}
    =
    {
    \prod_{k=1}^{\infty} 1 - x^{2k}
    \over
    \prod_{k=1}^{\infty} 1 - x^k
    }
\]

>>> coefficients $ euler .^ 2 * recipSimple euler
[1,1,1,2,2,3,4,5,6,8,10,12,15,18,22,27,32,38,46,54,64,76,89,104,...

* The Ramanujan tau function (<https://oeis.org/A000594 OEIS A000594>):

>>> coefficients $ mulX (euler^24)
[0,1,-24,252,-1472,4830,-6048,-16744,84480,-113643,-115920,534612,...

* The ratios of successive terms of Stern's diatomic series
  (<https://oeis.org/A002487 OEIS A002487>) define an enumeration of
  the non-negative rationals.  The generating function of the sequence
  satisfies \(x A(x) = (1 + x + x^2)A(x^2)\), which, with a little
  manipulation, yields

>>> let sternAux = 1 + mulX (1 + fromCoefficients [1,1,1]*sternAux .^ 2)
>>> coefficients $ mulX (1 + mulX sternAux)
[0,1,1,2,1,3,2,3,1,4,3,5,2,5,3,4,1,5,4,7,3,8,5,7,2,7,5,8,3,7,4,5,...

-}

{- |
@'linearRecurrence' a b@, where @a@ and @b@ are lists of length \(k\),
yields an ordinary generating function for the sequence \( \{c_n\} \)
defined by the recurrence

\[
    \begin{align}
        c_0 & = b_0 \\
        & \vdots \\
        c_{k-1} & = b_{k-1} \\
        c_{n+k} & = a_0 c_{n+k-1} + \cdots + a_{k-1} c_n
    \end{align}
\]

The generating function for this sequence is

\[ c(z) = {b(z) (1 - z a(z)) \mod z^k \over 1 - z a(z)} \]

=== __Examples__

* Fibonacci numbers (<https://oeis.org/A000045 OEIS A000045>):

\[
    \begin{align}
        F_0 & = 0 \\
        F_1 & = 1 \\
        F_{n+2} & = F_{n+1} + F_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [1,1] [0,1]
[0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,...

* Lucas numbers (<https://oeis.org/A000032 OEIS A000032>):

\[
    \begin{align}
        L_0 & = 2 \\
        L_1 & = 1 \\
        L_{n+2} & = L_{n+1} + L_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [1,1] [2,1]
[2,1,3,4,7,11,18,29,47,76,123,199,322,521,843,1364,2207,3571,5778,...

* Pell numbers (<https://oeis.org/A000129 OEIS A000129>):

\[
    \begin{align}
        P_0 & = 0 \\
        P_1 & = 1 \\
        P_{n+2} & = 2 P_{n+1} + P_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [2,1] [0,1]
[0,1,2,5,12,29,70,169,408,985,2378,5741,13860,33461,80782,195025,...

* Pell-Lucas numbers (<https://oeis.org/A002203 OEIS A002203>):

\[
    \begin{align}
        Q_0 & = Q_1 = 2 \\
        Q_{n+2} & = 2 Q_{n+1} + Q_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [2,1] [2,2]
[2,2,6,14,34,82,198,478,1154,2786,6726,16238,39202,94642,228486,...

* Padovan sequence (<https://oeis.org/A000931 OEIS A000931>):

\[
    \begin{align}
        P_0 & = P_1 = P_2 = 1 \\
        P_{n+3} & = P_{n+1} + P_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [0,1,1] [1,1,1]
[1,1,1,2,2,3,4,5,7,9,12,16,21,28,37,49,65,86,114,151,...

* Perrin sequence (<https://oeis.org/A001608 OEIS A001608>):

\[
    \begin{align}
        P_0 & = 3 \\
        P_1 & = 0 \\
        P_2 & = 2 \\
        P_{n+3} & = P_{n+1} + P_n
    \end{align}
\]

>>> coefficients $ linearRecurrence [0,1,1] [3,0,2]
[3,0,2,3,2,5,5,7,10,12,17,22,29,39,51,68,90,119,158,209,...

-}
linearRecurrence :: (Ring a) => [a] -> [a] -> PowerSeries a
linearRecurrence :: forall a. Ring a => [a] -> [a] -> PowerSeries a
linearRecurrence [a]
recurrence [a]
initial = PowerSeries a -> PowerSeries a -> PowerSeries a
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
divOneMinus PowerSeries a
r PowerSeries a
q
  where
    k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
recurrence) ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial)
    p :: PowerSeries a
p = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients [a]
initial
    q :: PowerSeries a
q = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients [a]
recurrence
    r :: PowerSeries a
r = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
k (PowerSeries a -> [a]
forall a. AdditiveMonoid a => PowerSeries a -> [a]
coefficients (PowerSeries a
pPowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
*(PowerSeries a
1 PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. AbelianGroup a => a -> a -> a
- PowerSeries a
forall a. Semiring a => PowerSeries a
identityPowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
*PowerSeries a
q))))

{- |
Infinite formal languages: the \(n\)th element of the sequence is
a finite set of strings of length \(n\).

=== __Examples__

* Addition is union of languages:

>>> coefficients $ string "ab" + string "c"
[fromList [],
 fromList ["c"],
 fromList ["ab"],
 fromList [],
 fromList [],...

* Multiplication is concatenation:

>>> coefficients $ (string "ab" + string "c") * (string "d" + string "ef")
[fromList [],
 fromList [],
 fromList ["cd"],
 fromList ["abd","cef"],
 fromList ["abef"],
 fromList [],...

* Kleene star is a sum of powers:

>>> let star l = recipOneMinus `compose` mulX (divX l)
>>> coefficients $ star $ string "ab" + string "c"
[fromList [""],
 fromList ["c"],
 fromList ["ab","cc"],
 fromList ["abc","cab","ccc"],
 fromList ["abab","abcc","cabc","ccab","cccc"],...

-}
type GradedLanguage a = PowerSeries (FiniteLanguage a)

-- | A language consisting of a single string.
string :: (Ord a) => [a] -> GradedLanguage a
string :: forall a. Ord a => [a] -> GradedLanguage a
string [a]
s = [FiniteLanguage a] -> PowerSeries (FiniteLanguage a)
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients (Int -> FiniteLanguage a -> [FiniteLanguage a]
forall a. Int -> a -> [a]
replicate ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
s) FiniteLanguage a
forall a. AdditiveMonoid a => a
zero [FiniteLanguage a] -> [FiniteLanguage a] -> [FiniteLanguage a]
forall a. [a] -> [a] -> [a]
++ [[a] -> FiniteLanguage a
forall a. a -> FiniteSet a
singleton [a]
s])

{- $polynomial_sequences
A polynomial sequence is a sequence of polynomials \(\{p_n(x)\}\)
where each polynomial \(p_n\) has degree at most \(n\), so that
the coefficients form a triangle.  Such sequences can be defined
by ordinary generating functions of the form

\[
\sum_{n=0}^\infty p_n(x) t^n
\]

=== Examples

For convenience, define

>>> let x = Data.YAP.Polynomial.identity

The sequence of polynomials with binomial coefficients
(Pascal's triangle, <https://oeis.org/A007318 OEIS A007318>)
has the generating
function \( {1 \over 1 - (1+x)t} = \sum_{n=0}^\infty {(1+x)}^n t^n \):

>>> coefficients $ compose recipOneMinus $ mulX $ constant $ 1 + x
[fromCoefficients [1],
 fromCoefficients [1,1],
 fromCoefficients [1,2,1],
 fromCoefficients [1,3,3,1],
 fromCoefficients [1,4,6,4,1],
 fromCoefficients [1,5,10,10,5,1],...

Chebyshev polynomials of the first kind
(<https://oeis.org/A053120 OEIS A053120>)
satisfy \( \cos (n \theta) = T_n(\cos \theta) \).
They are defined by the linear recurrence

\[
    \begin{align}
        T_0(x) & = 1 \\
        T_1(x) & = x \\
        T_{n+2}(x) & = 2 x\,T_{n+1}(x) - T_n(x)
    \end{align}
\]

or equivalently by the generating function \(1 - xt \over 1 - 2xt + t^2\).

>>> coefficients $ linearRecurrence [2*x, -1] [1, x]
[fromCoefficients [1],
 fromCoefficients [0,1],
 fromCoefficients [-1,0,2],
 fromCoefficients [0,-3,0,4],
 fromCoefficients [1,0,-8,0,8],
 fromCoefficients [0,5,0,-20,0,16],...

Chebyshev polynomials of the second kind
(<https://oeis.org/A053117 OEIS A053117>)
satisfy \( \sin\big((n + 1)\theta\big) = U_n(\cos \theta) \sin \theta \).
They are defined by the linear recurrence

\[
    \begin{align}
        U_0(x) & = 1 \\
        U_1(x) & = 2 x \\
        U_{n+2}(x) & = 2 x\,U_{n+1}(x) - U_n(x)
    \end{align}
\]

or equivalently by the generating function \(1 \over 1 - 2xt + t^2\).

>>> coefficients $ linearRecurrence [2*x, -1] [1, 2*x]
[fromCoefficients [1],
 fromCoefficients [0,2],
 fromCoefficients [-1,0,4],
 fromCoefficients [0,-4,0,8],
 fromCoefficients [1,0,-12,0,16],
 fromCoefficients [0,6,0,-32,0,32],...

The sequence of Legendre polynomials
(numerators <https://oeis.org/A100258 OEIS A100258>)
has the generating function \( 1 \over \sqrt{1 - 2xt + t^2} \):

>>> coefficients $ powerOnePlus (-0.5) `compose` fromCoefficients [0, -2*x, 1]
[fromCoefficients [1 % 1],
 fromCoefficients [0 % 1,1 % 1],
 fromCoefficients [(-1) % 2,0 % 1,3 % 2],
 fromCoefficients [0 % 1,(-3) % 2,0 % 1,5 % 2],
 fromCoefficients [3 % 8,0 % 1,(-15) % 4,0 % 1,35 % 8],
 fromCoefficients [0 % 1,15 % 8,0 % 1,(-35) % 4,0 % 1,63 % 8],...

The triangle of partition numbers (<https://oeis.org/A008284 OEIS A008284>)
has generating function

\[
	\exp\left( \prod_{n=1}^\infty \frac{x^n}{n} \frac{t^n}{1-t^n} \right) - 1
	= \sum_{n=0}^\infty \left( \sum_{k=0}^\infty T(n, k) x^k \right) t^n
\]

where \( T(n, k) \) is the number of partitions of \(n\) in which the
greatest is \(k\):

>>> map (mapAdditive Data.Ratio.numerator) $ coefficients $ compose expS (lambertTransform (integral (divX monomials))) - 1
[fromCoefficients [],
 fromCoefficients [0,1],
 fromCoefficients [0,1,1],
 fromCoefficients [0,1,1,1],
 fromCoefficients [0,1,2,1,1],
 fromCoefficients [0,1,2,2,1,1],
 fromCoefficients [0,1,3,3,2,1,1],...

-}

-- | Promotion of a power series to a power series of constant polynomials
constants :: (AdditiveMonoid a) => PowerSeries a -> PowerSeries (Polynomial a)
constants :: forall a.
AdditiveMonoid a =>
PowerSeries a -> PowerSeries (Polynomial a)
constants PowerSeries a
f = (a -> Polynomial a) -> PowerSeries a -> PowerSeries (Polynomial a)
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive a -> Polynomial a
forall a. a -> Polynomial a
Poly.constant PowerSeries a
f

-- | \(1 \over 1 - xt\), the generating function for the sequence of
-- monomials \(\{x^n\}\)
monomials :: (Eq a, Semiring a) => PowerSeries (Polynomial a)
monomials :: forall a. (Eq a, Semiring a) => PowerSeries (Polynomial a)
monomials = PowerSeries (Polynomial a) -> PowerSeries (Polynomial a)
forall a. Semiring a => PowerSeries a -> PowerSeries a
recipOneMinusOf (Polynomial a -> PowerSeries (Polynomial a)
forall a. AdditiveMonoid a => a -> PowerSeries a
constant Polynomial a
forall a. Semiring a => Polynomial a
Poly.identity)

{- | The Riordan array generated by series \(f\) and \(g\), with \(g(0) = 0\),
is a polynomial sequence with ordinary generating function

\[
    {f(t) \over 1 - x g(t)} =
        f(t) \left( 1 + x g(t) + (x g(t))^2 + \cdots \right)
\]

which can be constructed as
@`constants` f * `compose` `monomials` (`constants` g)@.

Setting \(x = 1\) yields \( {f(t) \over 1 - g(t)} \), the ordinary
generating function of the sequence of sums of the coefficients of
the rows.

Since \(g(t)\) is required have a constant term of zero, and thus
has the form \(g(t) = x h(t)\), some authors use the equivalent definition

\[
    {f(t) \over 1 - x t h(t)} =
        f(t) \left( 1 + x t h(t) + (x t h(t))^2 + \cdots \right)
\]

=== __Examples__

* The sequence of polynomials with binomial coefficients
(Pascal's triangle, <https://oeis.org/A007318 OEIS A007318>):

>>> coefficients $ riordan recipOneMinus (mulX recipOneMinus)
[fromCoefficients [1],
 fromCoefficients [1,1],
 fromCoefficients [1,2,1],
 fromCoefficients [1,3,3,1],
 fromCoefficients [1,4,6,4,1],
 fromCoefficients [1,5,10,10,5,1],...

* The Catalan triangle (<https://oeis.org/A033184 OEIS A033184>):

>>> coefficients $ riordan catalan (mulX catalan)
[fromCoefficients [1],
 fromCoefficients [1,1],
 fromCoefficients [2,2,1],
 fromCoefficients [5,5,3,1],
 fromCoefficients [14,14,9,4,1],
 fromCoefficients [42,42,28,14,5,1],...

-}
riordan :: (Semiring a) =>
    PowerSeries a -> PowerSeries a -> PowerSeries (Polynomial a)
riordan :: forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries (Polynomial a)
riordan PowerSeries a
f PowerSeries a
g =
    (a -> Polynomial a) -> PowerSeries a -> PowerSeries (Polynomial a)
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive a -> Polynomial a
forall a. a -> Polynomial a
Poly.constant PowerSeries a
f PowerSeries (Polynomial a)
-> PowerSeries (Polynomial a) -> PowerSeries (Polynomial a)
forall a.
Semiring a =>
PowerSeries a -> PowerSeries a -> PowerSeries a
`divOneMinus` (a -> Polynomial a) -> PowerSeries a -> PowerSeries (Polynomial a)
forall a b.
(AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> PowerSeries a -> PowerSeries b
forall (f :: * -> *) a b.
(AdditiveFunctor f, AdditiveMonoid a, AdditiveMonoid b) =>
(a -> b) -> f a -> f b
mapAdditive a -> Polynomial a
forall a. AdditiveMonoid a => a -> Polynomial a
Poly.linear (PowerSeries a -> PowerSeries a
forall a. PowerSeries a -> PowerSeries a
divX PowerSeries a
g)

-- This is the form used by
-- * Luschny, discarding constant term of g(t)
-- * Riordan Group book, requiring constant term of g(t) is zero
--
-- Sprugnoli uses f(t) \over 1 - t x g(t), so defined for any series,
-- and corresponding to the Bell transform rather than a Sheffer sequence.

{- $pgf
A power series whose coefficients sum to 1 may encode the probabilites of
each possible value of a natural number-valued random variable \(X\),
i.e. \( p_k = \Pr(X = k) \).
Such a series is then the probability generating function \(E[z^X]\).

* `one` is the PGF for the constant random variable with value 0.

* `identity` is the PGF for the constant random variable with value 1.

* If @pX@ is the PGF for the random variable \(X\), then the PGF for
  \(c X\) for natural number \(c\) is @pX .^ c@.

* The product of PGFs for independent random variables \(X\) and \(Y\)
  is the PGF for \(X+Y\).

* The sum of the coefficients of the derivative of a PGF is the expected
  value of the random variable.
-}

-- | The distribution of a single Bernoulli trial with probability \(p\)
-- of success has generating function \((1-p) + p z\).
--
-- Powers of this distribution yield the binomial distribution.
bernoulliPGF :: (Ring a) => a -> PowerSeries a
bernoulliPGF :: forall a. Ring a => a -> PowerSeries a
bernoulliPGF a
p = [a] -> PowerSeries a
forall a. AdditiveMonoid a => [a] -> PowerSeries a
fromCoefficients [a
1a -> a -> a
forall a. AbelianGroup a => a -> a -> a
-a
p, a
p]

-- | The geometric distribution, of the number of Bernoulli trials with
-- probability \(p\) before the first success,
-- has generating function \( {p \over 1 - (1-p)z} \).
--
-- Powers of this distribution yield the negative binomial or Pascal
-- distribution.
geometricPGF :: (Ring a) => a -> PowerSeries a
geometricPGF :: forall a. Ring a => a -> PowerSeries a
geometricPGF a
p = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant a
p PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* (PowerSeries a
forall a. Semiring a => PowerSeries a
recipOneMinus PowerSeries a -> a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> a -> PowerSeries a
.* (a
1a -> a -> a
forall a. AbelianGroup a => a -> a -> a
-a
p))

-- | The Poisson distribution with parameter \(\lambda\), of the number
-- of events with mean rate \(\lambda\) occurring in an interval,
-- has generating function \( e^{\lambda(z-1)} \).
poissonPGF :: (Floating a) => a -> PowerSeries a
poissonPGF :: forall a. Floating a => a -> PowerSeries a
poissonPGF a
l = a -> PowerSeries a
forall a. AdditiveMonoid a => a -> PowerSeries a
constant (a -> a
forall a. Floating a => a -> a
exp a
l) PowerSeries a -> PowerSeries a -> PowerSeries a
forall a. Semiring a => a -> a -> a
* PowerSeries a
forall a. FromRational a => PowerSeries a
expS PowerSeries a -> a -> PowerSeries a
forall a. Semiring a => PowerSeries a -> a -> PowerSeries a
.* a
l