(defn bulirsch-stoer
"Takes
- a (potentially lazy) sequence of `points` of the form `[x (f x)]` and
- a point `x` to interpolate
and generates a lazy sequence of approximations of `P(x)`. Each entry in the
return sequence incorporates one more point from `points` into the `P(x)`
estimate.
`P(x)` is rational function fit (using the Bulirsch-Stoer algorithm, of
similar style to Neville's algorithm described
in [[emmy.numerical.interpolate.polynomial]]) to every point in the
supplied sequence `points`.
\"The Bulirsch-Stoer algorithm produces the so-called diagonal rational
function, with the degrees of numerator and denominator equal (if m is even)
or with the degree of the denominator larger by one if m is odd.\" ~ Press,
Numerical Recipes, p105
The implementation follows [Equation 3.2.3 on on page 105 of
Press](http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf).
### Column
If you supply an integer for the third (optional) `column` argument,
`bulirsch-stoer` will return that /column/ offset the interpolation tableau
instead of the first row. This will give you a sequence of nth-order
polynomial approximations taken between point `i` and the next `n` points.
As a reminder, this is the shape of the tableau:
```
p0 p01 p012 p0123 p01234
p1 p12 p123 p1234 .
p2 p23 p234 . .
p3 p34 . . .
p4 . . . .
```
So supplying a `column` of `1` gives a sequence of 2-point approximations
between pairs of points; `2` gives 3-point approximations between successive
triplets, etc.
References:
- Stoer & Bulirsch, ['Introduction to Numerical Analysis'](https://www.amazon.com/Introduction-Numerical-Analysis-Applied-Mathematics/dp/144193006X)
- [PDF of the same reference](http://www.math.uni.wroc.pl/~olech/metnum2/Podreczniki/(eBook)%20Introduction%20to%20Numerical%20Analysis%20-%20J.Stoer,R.Bulirsch.pdf)
- Press's Numerical Recipes (p105), [Section 3.2](http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-2.pdf)"
([points x] (bulirsch-stoer points x nil))
([points x column]
(let [merge (bs-merge x)
tableau (pi/tableau-fn bs-prepare merge points)]
(bs-present
(if column
(nth tableau column)
(pi/first-terms tableau))))))