You can derive the formula for this using the Chinese Remainder Theorem.<p>Consider the problem of finding a polynomial p(x) such that p(x_i) = y_i for some finite set of (x_i) and (y_i). The classic approach to this is <i>Lagrange interpolation</i>, which can be seen as a special case of the Chinese Remainder Theorem. Essentially, you are trying to find a p which satisfies the following modular congruences:<p><pre><code> p(x) == y_i (mod (x - x_i))
</code></pre>
The constructive proof of CRT, specialised to this problem, yields an algorithm equivalent to Lagrange interpolation. To extend Lagrange interpolation so as to include derivative information y'_i at each x_i, change that congruence to:<p><pre><code> p(x) == y_i + y'_i (x - x_i) (mod (x - x_i)^2)
</code></pre>
And to give second derivative information, change the modular congruence to:<p><pre><code> p(x) == y_i + y'_i (x - x_i) + y''_i (x - x_i)^2 / 2! (mod (x - x_i)^3)
</code></pre>
And so on...