Power Series

This namespace contains an implementation of two data types:

• Series, which represents a generic infinite series of arbitrary values, and
• PowerSeries, a series that represents a power series in a single variable; in other words, a series where the nth entry is interpreted as the coefficient of Loading...:

[a b c d ...] Loading...

We'll proceed by building up implementations of the arithmetic operations +, -, *, / and a few others on bare Clojure lazy sequences, and then introduce deftype wrappers so that we can install these types into the generic dispatch system.

The implementation follows Doug McIlroy's beautiful paper, "Power Series, Power Serious". Doug also has a 10-line version in Haskell on his website.

Okay, let's proceed, in roughly the same order as the paper.

Sequence Operations

A 'series' is an infinite sequence of numbers, represented by Clojure's lazy sequence. First, a function ->series that takes some existing sequence, finite or infinite, and coerces it to an infinite seq by concatenating it with an infinite sequence of zeros. (We use g/zero-like so that everything plays nicely with generic arithmetic.)

This works as expected:

The core observation we'll use in the following definitions (courtesy of McIlroy) is that a power series Loading... in a variable Loading...:

Decomposes into a head element Loading... plus a tail series, multiplied by Loading...:

We'll use this observation to derive the more complicated sequence operations below.

Negation

To negate a series, negate each element:

Addition

We can derive series addition by expanding the series Loading... and Loading... into head and tail and rearranging terms:

This is particularly straightforward in Clojure, where map already merges sequences elementwise:

A constant is a series with its first element populated, all zeros otherwise. To add a constant to another series we only need add it to the first element. Here are two versions, constant-on-left vs constant-on-right:

Subtraction

Subtraction comes for free from the two definitions above:

Subtract a constant from a sequence by subtracting it from the first element:

To subtract a sequence from a constant, subtract the first element as before, but negate the tail of the sequence:

Multiplication

What does it mean to multiply two infinite sequences? As McIlroy notes, multiplication is where the lazy-sequence-based approach really comes into its own.

First, the simple cases of multiplication by a scalar on either side of a sequence:

To multiply sequences, first recall from above that we can decompose each sequence Loading... and Loading... into a head and tail.

Mutiply the expanded representations out and rearrange terms:

Loading... appears on the left and the right, so use an inner function that closes over Loading... to simplify matters, and rewrite the above definition in Clojure:

This works just fine on two infinite sequences:

NOTE This is also called the "Cauchy Product" of the two sequences: https://en.wikipedia.org/wiki/Cauchy_product The description on the Wikipedia page has complicated index tracking that simply doesn't come in to play with the stream-based approach. Amazing!

Division

The quotient Loading... of Loading... and Loading... should satisfy:

From McIlroy, first expand out Loading..., Loading... and one instance of Loading...:

Look at just the constant terms and note that Loading....

Consider the terms multiplied by Loading... and solve for Loading...:

There are two special cases to consider:

• If Loading..., Loading... can only succeed if Loading...; in this case, $Q = {F_1 \over G1}$, from the larger formula above.
• If Loading..., Loading...

Encoded in Clojure:

A simple example shows success:

Reciprocal

We could generate the reciprocal of Loading... by dividing Loading... by Loading.... Page 21 of an earlier paper by McIlroy gives us a more direct formula.

We want Loading... such that Loading.... Expand Loading...:

Solve for R:

A recursive definition is no problem in the stream abstraction:

This definition of invert matches the more straightforward division implementation:

An example:

Division of a constant by a series comes easily from our previous multiplication definitions and invert:

It's not obvious that this works:

But we can recover the initial series:

To divide a series by a constant, divide each element of the series:

Division by a constant undoes multiplication by a constant:

Functional Composition

To compose two series Loading... and Loading... means to create a new series Loading.... Derive this by substituting Loading... for Loading... in the expansion of Loading...:

For the stream-based calculation to work, we need to be able to calculate the head element and attach it to an infinite tail; unless Loading... above the head element depends on Loading..., an infinite sequence.

If Loading... the calculation simplifies:

In Clojure, using an inner function that captures Loading...:

Composing Loading... should square all Loading...s, and give us a sequence of only even powers:

Reversion

The functional inverse of a power series Loading... is a series Loading... that satisfies Loading....

Following McIlroy, we expand Loading... (substituting Loading... for Loading...) and one occurrence of Loading...:

Just like in the composition derivation, in the general case the head term depends on an infinite sequence. Set Loading... to address this:

For this to work, the constant Loading... must be 0 as well, hence

Loading...

This works as an implementation because Loading.... Loading... is allowed to reference Loading... thanks to the stream-based approach:

An example, inverting a series starting with 0:

Series Calculus

Derivatives of power series are simple and mechanical:

Implies that all entries shift left by 1, and each new entry gets multiplied by its former index (i.e., its new index plus 1).

The definite integral Loading... is similar. To take the anti-derivative of each term, move it to the right by appending a constant term onto the sequence and divide each element by its new position:

With a custom constant term:

By default, the constant term is 0:

Exponentiation

Exponentiation of a power series by some integer is simply repeated multiplication. The implementation here is more efficient the iterating seq:*, and handles negative exponent terms by inverting the original sequence.

We can use expt to verify that Loading... expands to Loading...:

Square Root of a Series

The square root of a series Loading... is a series Loading... such that Loading.... We can find this using our calculus methods from above:

or

When the head term of Loading... is nonzero, i.e., Loading..., the head of Loading... must be Loading... for the multiplication to work out.

Integrate both sides:

One optimization appears if the first two terms of Loading... vanish, i.e., Loading.... In this case Loading....

Here it is in Clojure:

And a test that we can recover the naturals:

We can maintain precision of the first element is the square of a rational number:

We get a correct result if the sequence starts with 0, 0:

Examples

Power series computations can encode polynomial computations. Encoding Loading... as a power series returns the correct result:

Encoding Loading... returns the power series Loading... which sums to that value in its region of convergence:

Loading... is the derivative of the above series:

Various Power Series

With the above primitives we can define a number of series with somewhat astonishing brevity.

Loading... is its own derivative, so Loading...:

This bare definition is enough to generate the power series for Loading...:

Loading... and Loading... afford recursive definitions. Loading... and $D(cos) = -sin$, so (with appropriate constant terms added) on:

tangent and secant come easily from these:

Reversion lets us define arcsine from sine:

These two are less elegant, perhaps:

The hyperbolic trig functions are defined in a similar way:

https://en.wikipedia.org/wiki/Mercator_series

Generating Functions

Catalan numbers

These are a few more examples from McIlroy's "Power Serious" paper, presented here without context. (If you have the energy to write about these, please feel free and send us a PR!)

ordered trees...

The catalan numbers again!

See here for the recurrence relation: https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula