yap-examples-0.1: examples of the algebraic classes in the yap package
Copyright(c) Ross Paterson 2021
LicenseBSD-style (see the file LICENSE)
Maintainer[email protected]
Stabilityprovisional
Portabilityportable
Safe HaskellNone
LanguageHaskell2010

Data.YAP.PowerSeries

Description

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

Formal power series

data PowerSeries a Source #

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\).

Instances

Instances details
AdditiveFunctor PowerSeries Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

mapAdditive :: (AdditiveMonoid a, AdditiveMonoid b) => (a -> b) -> PowerSeries a -> PowerSeries b #

Semiring a => Differentiable (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

derivative :: PowerSeries a -> PowerSeries a #

FromRational a => Integrable (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

integral :: PowerSeries a -> PowerSeries a #

AbelianGroup a => AbelianGroup (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

(-) :: PowerSeries a -> PowerSeries a -> PowerSeries a #

negate :: PowerSeries a -> PowerSeries a #

gtimes :: (AbelianGroup b, ToInteger b) => b -> PowerSeries a -> PowerSeries a #

AdditiveMonoid a => AdditiveMonoid (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

(+) :: PowerSeries a -> PowerSeries a -> PowerSeries a #

zero :: PowerSeries a #

atimes :: ToInteger b => b -> PowerSeries a -> PowerSeries a #

(Eq a, Field a) => Euclidean (PowerSeries a) Source #

mod x y is non-zero if y has more leading zeros than x. See also modulus, for the remainder as a polynomial.

Instance details

Defined in Data.YAP.PowerSeries

Methods

div :: PowerSeries a -> PowerSeries a -> PowerSeries a #

mod :: PowerSeries a -> PowerSeries a -> PowerSeries a #

divMod :: PowerSeries a -> PowerSeries a -> (PowerSeries a, PowerSeries a) #

euclideanNorm :: PowerSeries a -> Natural #

FromRational a => FromRational (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

fromRational :: Rational -> PowerSeries a #

Ring a => Ring (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

fromInteger :: Integer -> PowerSeries a #

Semiring a => Semiring (PowerSeries a) Source # 
Instance details

Defined in Data.YAP.PowerSeries

Methods

(*) :: PowerSeries a -> PowerSeries a -> PowerSeries a #

one :: PowerSeries a #

fromNatural :: Natural -> PowerSeries a #

rescale :: PowerSeries a -> PowerSeries a -> (PowerSeries a, PowerSeries a, PowerSeries a -> PowerSeries a) #

(Eq a, Field a) => StandardAssociate (PowerSeries a) Source #

Units have non-zero leading coefficients. Standard associates are monomials of the form \(x^n\). stdUnit and stdRecip are not defined on zero.

Instance details

Defined in Data.YAP.PowerSeries

Methods

stdAssociate :: PowerSeries a -> PowerSeries a #

stdUnit :: PowerSeries a -> PowerSeries a #

stdRecip :: PowerSeries a -> PowerSeries a #

constant :: AdditiveMonoid a => a -> PowerSeries a Source #

Power series representing a constant value \(c\)

fromPolynomial :: AdditiveMonoid a => Polynomial a -> PowerSeries a Source #

Polynomial considered as a power series

fromCoefficients :: AdditiveMonoid a => [a] -> PowerSeries a Source #

Power series formed from a list of coefficients. If the list is finite, the remaining coefficients are zero.

fromDerivatives :: FromRational a => [a] -> PowerSeries a Source #

Power series formed from a list of derivatives at zero. If the list is finite, the remaining derivatives are zero.

Queries

atZero :: AdditiveMonoid a => PowerSeries a -> a Source #

Value of the function at zero

coefficients :: AdditiveMonoid a => PowerSeries a -> [a] Source #

The coefficients of the power series, least significant first.

derivatives :: Semiring a => PowerSeries a -> [a] Source #

A list whose \(n\)th entry is the value of the \(n\)th derivative at zero.

flatten :: AdditiveMonoid a => PowerSeries (PowerSeries a) -> PowerSeries a Source #

Reduce a power series \( \sum_{n=0} p_n(x) x^n \) with power series coefficients to a simple power series.

toPolynomial :: AdditiveMonoid a => Int -> PowerSeries a -> Polynomial a Source #

Polynomial with the first \(n\) coefficients, equivalent to the remainder of division by \(x^n\).

modulus :: (Eq a, AdditiveMonoid a) => PowerSeries a -> PowerSeries a -> Polynomial a Source #

The remainder of division by a non-zero power series is a polynomial.

order :: (Eq a, AdditiveMonoid a) => PowerSeries a -> Int Source #

Degree of the first non-zero coefficient. Undefined if the argument is zero.

distance :: (Eq a, AdditiveMonoid a) => PowerSeries a -> PowerSeries a -> Rational Source #

Metric on power series (undefined if the arguments are equal)

Approximations

approximations :: Semiring a => PowerSeries a -> a -> [a] Source #

The infinite list of evaluations of truncations of the power series.

approximationsWith :: (AdditiveMonoid a, Semiring b, AdditiveMonoid c) => (a -> b -> c) -> PowerSeries a -> b -> [c] Source #

The infinite list of evaluations of truncations of the power series, given a function to combine coefficients and powers.

padeApproximants :: (Eq a, Field a) => Int -> PowerSeries a -> [RationalFunction a] Source #

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) \]

Composition

identity :: Semiring a => PowerSeries a Source #

Identity function, i.e. the indeterminate

compose :: (Eq a, Semiring a) => PowerSeries a -> PowerSeries a -> PowerSeries a Source #

composition \(f(g(x))\). This is only defined if \(g_0\) is 0.

inverse :: (Eq a, Field a) => PowerSeries a -> PowerSeries a Source #

pre-inverse under compose:

This is only defined if \(f_0\) is zero and \(f_1\) is non-zero.

inverseSimple :: (Eq a, Ring a) => PowerSeries a -> PowerSeries a Source #

Special case of inverse that is only valid if the constant term is 0 and the coefficient of the linear term is 1.

Special compositions

(.*) :: Semiring a => PowerSeries a -> a -> PowerSeries a infixl 9 Source #

Maps a function \(f(x)\) to \(f(a x)\), equivalent to composition with a*identity.

(.^) :: AdditiveMonoid a => PowerSeries a -> Int -> PowerSeries a infixl 9 Source #

Maps a function \(f(x)\) to \(f(x^k)\) for positive \(k\), equivalent to composition with identity^k.

Other transforms

divX :: PowerSeries a -> PowerSeries a Source #

Discard the constant term and divide by the indeterminate

mulX :: AdditiveMonoid a => PowerSeries a -> PowerSeries a Source #

Multiply by the indeterminate

shift :: AdditiveMonoid a => Int -> PowerSeries a -> PowerSeries a Source #

Shift the power series by \(k\) (which may be negative): multiply by \(x^k\) and discard any terms with negative exponents.

recipSimple :: Ring a => PowerSeries a -> PowerSeries a Source #

Special case of recip that is only valid if the constant term is 1.

divSimple :: (Eq a, Ring a) => PowerSeries a -> PowerSeries a -> PowerSeries a Source #

Special case of division that is only valid if the constant term of the divisor is 1.

Special series

recipOneMinus :: Semiring a => PowerSeries a Source #

\({1 \over 1-x} = 1 + x + x^2 + \ldots \) (geometric series), converges for \(|x| < 1\). derivatives yields the factorials.

expS :: FromRational a => PowerSeries a Source #

\(e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots\), converges everywhere. All derivatives are 1.

logOnePlus :: FromRational a => PowerSeries a Source #

\(\log (1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \ldots \), converges for \(-1 < x \leq 1\)

powerOnePlus :: FromRational a => a -> PowerSeries a Source #

\( (1 + x)^a = \sum_{n=0}^\infty \binom{a}{n} x^n \) (binomial series), converges for \(|x| < 1\)

euler :: Ring a => PowerSeries a Source #

Euler function (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} \]

Trigonometric functions

sinS :: FromRational a => PowerSeries a Source #

\(\sin x\), converges everywhere (but faster for smaller \(x\))

cosS :: FromRational a => PowerSeries a Source #

\(\cos x\), converges everywhere (but faster for smaller \(x\))

secS :: FromRational a => PowerSeries a Source #

\(1 \over \cos x\), converges for \(|x| < {\pi \over 2}\)

tanS :: FromRational a => PowerSeries a Source #

\(\tan x\), converges for \(|x| < {\pi \over 2}\)

asinS :: FromRational a => PowerSeries a Source #

\(\sin^{-1} x\), converges for \(|x| \leq 1\)

atanS :: FromRational a => PowerSeries a Source #

\(\tan^{-1} x\), converges for \(|x| \leq 1\)

Hyperbolic functions

sinhS :: FromRational a => PowerSeries a Source #

\(\sinh x\), converges everywhere (but faster for smaller \(x\))

coshS :: FromRational a => PowerSeries a Source #

\(\cos x\), converges everywhere (but faster for smaller \(x\))

sechS :: FromRational a => PowerSeries a Source #

\(1 \over \cosh x\), converges for \(|x| < {\pi \over 2}\)

tanhS :: FromRational a => PowerSeries a Source #

\(\tanh x\), converges for \(|x| < {\pi \over 2}\)

asinhS :: FromRational a => PowerSeries a Source #

\(\sinh^{-1} x\), converges for \(|x| \leq 1\)

atanhS :: FromRational a => PowerSeries a Source #

\(\tanh^{-1} x\), converges for \(|x| < 1\)

Continued fractions

stieltjes :: Semiring a => [a] -> PowerSeries a Source #

A Stieltjes fraction (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

Expand
  • The Catalan numbers \( C(n) = {(2n)! \over n! (n+1)!} \) (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,...

\[ (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 (OEIS A000110):
>>> coefficients $ stieltjes $ concat [[1, n] | n <- [1..]]
[1,1,2,5,15,52,203,877,4140,21147,115975,678570,4213597,27644437,...

jacobi :: Semiring a => [a] -> [a] -> PowerSeries a Source #

A Jacobi fraction (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

Expand
  • Using the above identity, yet another expression for the Catalan numbers (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 (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 (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 (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 (OEIS A006789):
>>> coefficients $ jacobi [1..] (repeat 1)
[1,1,2,5,14,43,143,509,1922,7651,31965,139685,636712,3020203,...

Sequence transforms

Binomial transform

binomialTransform :: Semiring a => PowerSeries a -> PowerSeries a Source #

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.

Examples

Expand
  • Maximal number of pieces obtained by cutting a circle with \(n\) lines (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 (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 (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,...

inverseBinomialTransform :: Ring a => PowerSeries a -> PowerSeries a Source #

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.

Lambert transform

lambertTransform :: AdditiveMonoid a => PowerSeries a -> PowerSeries a Source #

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.

Examples

Expand
  • 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,...

inverseLambertTransform :: Ring a => PowerSeries a -> PowerSeries a Source #

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.

Examples

Expand
>>> 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\) (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,...

Euler transform

eulerTransform :: (Eq a, FromRational a) => PowerSeries a -> PowerSeries a Source #

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

Expand

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 (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 (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 (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,...

inverseEulerTransform :: (Eq a, FromRational a) => PowerSeries a -> PowerSeries a Source #

The inverse of eulerTransform.

Ordinary generating functions

Number sequences

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}\) (OEIS A001008 over 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 (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 (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} \) (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 (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.)
>>> 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 (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,...

>>> 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 (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,...

Linear recurrences

linearRecurrence :: Ring a => [a] -> [a] -> PowerSeries a Source #

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

Expand

\[ \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,...

\[ \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,...

\[ \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,...

\[ \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,...

\[ \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,...

\[ \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,...

Formal languages

type GradedLanguage a = PowerSeries (FiniteLanguage a) Source #

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

Examples

Expand
  • 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"],...

string :: Ord a => [a] -> GradedLanguage a Source #

A language consisting of a single string.

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, 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 (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 (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 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 (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],...

constants :: AdditiveMonoid a => PowerSeries a -> PowerSeries (Polynomial a) Source #

Promotion of a power series to a power series of constant polynomials

monomials :: (Eq a, Semiring a) => PowerSeries (Polynomial a) Source #

\(1 \over 1 - xt\), the generating function for the sequence of monomials \(\{x^n\}\)

riordan :: Semiring a => PowerSeries a -> PowerSeries a -> PowerSeries (Polynomial a) Source #

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

Expand
  • The sequence of polynomials with binomial coefficients (Pascal's triangle, 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],...

>>> 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],...

delta :: Semiring a => [a] -> [a] -> PowerSeries (Polynomial a) Source #

Deléham's delta operator (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

Expand
>>> 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 (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 (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)\) (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],...

Probability generating functions

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.

bernoulliPGF :: Ring a => a -> PowerSeries a Source #

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.

geometricPGF :: Ring a => a -> PowerSeries a Source #

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.

poissonPGF :: Floating a => a -> PowerSeries a Source #

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)} \).