# Code dissection: polynomial extrapolation

| 1 Comment | No TrackBacks

Today's dissection is a bit mathematical but rather less so than the previous one. Again, I managed to save a bit (only 5% this time) when explaining it caused me to spot things I'd previously missed.

The challenge was to take as input `k a1 a2 … ak m`, fit the unique polynomial of degree `k-1` through the implicit points `(i, ai)`, and then extrapolate to the next `m` integers.

The key mathematical technique we'll be using is the finite difference operator. If we have a polynomial `p(x)` of degree `n` (and assuming that it's not `p(x) = 0`) then `Δp(x) = p(x+1) - p(x)` is a polynomial of degree `n-1`. (This can easily be demonstrated by the binomial theorem). So by repeated application we find that `Δnp(x)` is a polynomial of degree 0, i.e. a constant. Now, extrapolating constant functions is really easy, so we can use that to extrapolate the linear function `Δn-1p(x)`, and by extrapolating that we can extrapolate the quadratic, etc. By way of example, here's the difference table for `1 2 4 7`.

```1   2   4   7
1   2   3
1   1```

We can then extend the trailing edge from the bottom up to get the next value, 11:

```1   2   4   7   11
1   2   3   4
1   1   1```

So here's the full, undissected, code:

`~])\({[{.@-\}*])\}*;]-1%){0\{+.}/p]}*;`

Let's take this step by step.

Firstly the `~` evaluates the string we got from stdin, which happens to be in a very nice format for us. That leaves the integers `k a1 a2 … ak m` on the stack, and `]` gathers them into an array. `)\` pulls `m` from the end of the array and moves it to the bottom of the stack. We're now ready for the first of two main loops.

`({[{.@-\}*])\}*` pulls `k` from the start of the array (if it weren't there we'd have to do `.,`) and then executes `[{.@-\}*])\` that many times. The eventual goal of this first stage is to get the trailing edge of the difference table. This code looks quite strange: we have an array on the stack, so that `*` is a fold of the block over the array, but why wrap a fold in `[ ]`? The answer is twofold: unlike map (`%`), the values we produce won't be pulled off the stack after each execution of the block; and secondly, fold is designed for binary operations, so we save initialising the difference with `0`. The block itself is very simple: `.@-\` turns `x y` into `(y-x) y`. So we get the array of differences, i.e. the next row of the table; with a final element which is the last element of the row we were working from - exactly what we wanted for the trailing edge. `)\` extracts it and brings the array of differences back to the top for the next iteration of the outer loop.

Then there's a little bit of glue: `;]-1%)` discards the (now empty) array of differences; gathers the entire stack into an array; reverses it; and splits out the last element, which was previously the bottom of the stack, the value `m` which we dumped there so long ago. The effect of the reversal is that the array contains the trailing edge from the bottom up, which is the order we want for iteration. And that leads us into our second main loop. `{0\{+.}/p]}*` executes `0\{+.}/p]` `m` times, starting out each time with a single array, the trailing edge of the difference table, on the stack. This time we do need to supply an initial value and use a foreach rather than a fold because otherwise our "trailing edge" will get one element shorter each time round the loop. So we push a `0`, flip it behind the array, and for each element in the array we add and then duplicate, effectively mapping the array to a cumulative total of the array which is the next trailing edge. The duplicate of the final (i.e. top) element is precisely the value we need to print, and then we gather the entire stack back up into an array.

There's one final thing left to do: `;` to pop that last trailing edge and avoid it being printed. And there we have it.

## 1 Comment

Then there's a little bit of glue: ;]-1%) discards the (now empty) array of differences; gathers the entire stack into an array; reverses it; and splits out the last element, which was previously the bottom of the stack, the value m which we dumped there so long ago. The effect of the reversal is that the array contains the trailing edge from the bottom up, which is the order we want for iteration. And that leads us into our second main loop. {0\{+.}/p]}* executes 0\{+.}/p] m times, starting out each time with a single array, the trailing edge of the difference table, on the stack. This time we do need to supply an initial value and use a foreach rather than a fold because otherwise our "trailing edge" will get one element shorter each time round the loop. So we push a 0, flip it behind the array, and for each element in the array we add and then duplicate, effectively mapping the array to a cumulative total of the array which is the next trailing edge. The duplicate of the final (i.e. top) element is precisely the value we need to print, and then we gather the entire stack back up into an array.

### Pages

• Subscribe to this blog's feed