I was reading about Babbage's Difference engine the other day, and stumbled upon a very interesting application of the forward differences method. It turns out that if we get a sequence generated by a polynomial, under certain conditions we can find the generating polynomial from just a few elements in the sequence.

For example, 0, 1, 5, 12, 22, 35, 51... is a sequence known as the pentagonal numbers, and we can use this technique to figure out that the polynomial generates it [1].

## Notation

Let's start with some mathematical notation. We'll call the underlying function generating the sequence . In our example, , , , and so on.

- The
*first difference*is the sequence , , etc. - The
*second difference*is the sequence , etc. - In general, the k-th difference is: .
- As a starting condition in the induction of differences, we can say that itself is the 0-th difference.

## Difference table

We can construct the *difference table* for our sequence and observe its
properties. In a difference table, the first column is *n*, which runs from 0
to whatever number of elements we have for the sequence. The second column is
the values of the sequence at these *n*. Then come the first difference, the
second difference and so on. For our sample sequence we get the table:

Notice how at some point the differences become all-zero! We'll soon see why.

Obviously, we can construct such a table from only the column of - that's what we just did! A more interesting observation is that if we accept that all differences (columns) are 0 from a certain point, we can also construct this table from just the first row! For example, with and in hand, we know ; with and , we know etc. Try it!

## Inferring the polynomial's degree from the table

**Claim:** if has degree *k*, then the *k*-th difference column
in the table is constant.

**Proof:** this is a general *k*-th degree polynomial with coefficients
:

By definition of the first difference, if we expand the polynomial form and
perform the subtraction per power of *n*:

Using the binomial theorem we know that:

Therefore:

Now if we look at the sum we got for again, we'll notice
that in each term, the *j*-th power of *n* gets canceled out. This means that
the highest power of *n* in is going to be *k-1*.

We can similarly show that in , the highest power
of *n* is going to be *k-2*. Therefore, the *k*-th difference
will be constant .

Two observations for extra credit:

- Note that the claim goes one way -
*if*is a*k*-th degree polynomial, the*k*-th difference is constant. What we observe in the table is the*k*-th difference is constant, so can we infer that the function is a*k*-th degree polynomial? Not in the general case! The sequence could be generated by some higher-degree polynomial, or by a completely different kind of function. That said, since we assume is polynomial and are seeking to find the simplest (lowest degree) one, this inference is valid. - Did you notice the equivalence to derivatives of polynomials? The
*k*-th difference of a polynomial of degree*k*is constant... but the same is true for the*k*-th derivative! This is not by chance - since we have a discrete domain, differences play largely the same role as derivatives for continuous functions. This isn't a rigorous proof - but think about the definition of derivatives (the one with the limit) - what do you get when you take (also sometimes called*h*) and set it to 1?

## Finding the coefficients with the Newton polynomial

Polynomial interpolation can fit any N points
(with distinct *x* values) with a N-1 degree polynomial. One way of finding such
a polynomial was discovered by Isaac Newton and is called the
Newton polynomial.

Our problem of finding a polynomial that generates a given set of points can
be reduced to this interpolation problem, since we've just figured out the
degree of the generating polynomial! Looking at the difference table, we've
found when the difference becomes constant, and that gives us the polynomial's
degree *k*. So all we need is the first *k+1* points.

Here's how to develop the Newton polynomial from scratch; we'll start with the
first few coefficients and will then generalize for any *k*.

The Newton polynomial for our set of forward differences can be expressed as follows:

This polynomial is constructed in a clever way; notice that for any *p*, when
we calculate all the elements starting with will be
multiplied by zero and vanish. This helps us determine this polynomial's
coefficients in a gradual manner.

We'll start with the point , substituting it into the Newton polynomial:

This gives us the first coefficient . Next, let's look at :

Since we know that , we can infer that . Another way to express that is .

Continuing to :

We've substituted the values of and that we've found earlier; let's solve for now:

The last line's numerator is - by definition - :

We can keep going with this (feel free to do as an exercise), but the emerging generalization is that:

And a concise way to write Newton's polynomial for forward differences is:

Where is:

Note that we only use the differences for , meaning that we need just the first row of the difference table! Let's try it for the pentagonal numbers example.

First, we've determined that is a constant, so the degree
of the polynomial is 2. We only need to calculate until *i=2* in the sum:

Substituting the values from the difference table we have for *n=0*, we get:

Which is exactly what we expected!

## Another example

Let's work through another example, taking the sequence -8, -12, -6, 16, 60... We'll start by constructing the difference table:

The difference appears to be constant, so we can generate this sequence with a degree 3 polynomial. Let's use the first line of the table to construct the Newton polynomial for it.

Verifying that this polynomial generates our sequence as its first 5 elements is an easy exercise.

Note that given 5 elements, we can always find a 4th-degree polynomial fitting it. Here we found a 3rd-degree one, though, leveraging the technique of differences. This becomes more acute if we have more elements in the sequence - using differences it's often possible to find significantly simpler polynomials.

## Recap

For an integer sequence, if this sequence is generated by a polynomial we can figure out which polynomial it is - given enough elements. We start by constructing a difference table and noticing if a column becomes constant from some point on. This tells us the degree of the generating polynomial. With that in hand, we can use the Newton polynomial to discover a polynomial that generates the sequence.

## Resources

- This blog post was the original inspiration. The post explains how Babbage's difference engine worked, and makes off-hand remarks like "for a 2nd degree polynomial, we only need up to the second difference to know everything". This follow-up by the same author mentioned going from sequences back to polynomials.
- The old Encyclopedia of integer sequences has an intriguing section 2.5, which unfortunately presents several lemmas with no proofs.
- This wiki page from brilliant.org is the best single resource, though its proof of constructing the Newton polynomial is lacking, IMHO.
- Wikipedia: Divided differences and Newton polynomial pages.
- Newton's forward differences YouTube lecture.
- Knuth covers this
*very briefly*in section 4.6.4 of Part 2 of TAOCP, relegating the derivation of the Newton polynomial to an exercise that has a terse solution.

## Appendix: Another way to find the polynomial's coefficients

Here's another way to discover the coefficients of the polynomial; the following discusses the coefficients of the highest power, but can be generalized to other powers as well. Using Newton's polynomial is simpler, though, so this is just an appendix for some extra practice in manipulating such polynomials.

**Claim:** if has degree *k*, its coefficient
is equal to . Since we've proven that
is a constant, we can find the precise
value of this way.

**Proof:** Let's go back to the sum formulation of from
the previous proof.

Expanding the binomial, we get:

The highest power of *n* here is *k-1*; its coefficient comes only from the
first term of the binomial expansion for *j=k*, and is equal to .
In order not to deal with long sums, let's just focus on the highest-degree
term in :

Where represents other elements with lower powers of *k*, so we
don't care about them for this discussion [2]. Let's move on to the next
difference:

We've just found the coefficient of the highest power of *n* in
.
It's clear that if we continue doing this, will be:

In other words, , meaning that we can know
the coefficient of the highest power of *n* in from the *k*-th
difference .

[1] | Naturally, there's an infinitude of potential functions that generate this sequence; it's more precise to say we're looking for the simplest (lowest degree) polynomial that would do this. |

[2] | In the following calculation, we're playing loose with the
s to represent "anything with lower powers of n
that we don't care about". |