Bug in CubicSplineInterpolator integration function

[jamiescottie1]

The way the integral in the CubicSplineInterpolator is implemented is slightly off. We need to evaluate the integral $$\int_0^t (1-t')y_l + t'y_r + t'(1-t')((1-t')a + t'b)\mathrm{d}t'$$ which results in (see WolframAlpha where $l\to y_l$, $r\to y_r$, $k\to t$, $t\to t'$) $$\int_0^t\dots\mathrm{d}t' = \frac{t}{12}\left(at\left(3t^2-8t+6\right) + bt^2\left(4-3t\right) - 6ty_l + 6ty_r + 12y_l\right)$$

If we rewrite this as implemented in the current code, we get $$\int_0^t\dots\mathrm{d}t' = a\left(\frac{1}{4}t^4 - \frac{2}{3}t^3 + \frac{1}{2}t^2\right) + b\left(\frac{1}{3}t^3 - \frac{1}{4}t^4\right) + y_l\left(t - \frac{1}{2}t^2\right) + \frac{1}{2}y_rt$$

This is almost the same as what is implemented, with the only difference being the factor $\frac{1}{2}$ in front of the quadratic term $t^2$ that gets multiplied by $a$, which is missing in the current implementation.

---

Side note: For performance reasons the various integer powers should be replaced by explicit products, which are much more efficient than the general purpose pow.