> It is not generally known that to provide seven- decimal accuracy of the sine function, allowing third-order interpolation, only 15 entries are required [*] ...

That got me thinking: is that so? how do you do that?

Luckily, scipy has all sorts of interpolation tools, so one can just experiment a little.

Third-order interpolation -> that means cubic splines. How many points do you need for a cubic spline approximation of the sine function to have 7 decimal place accuracy (in other words the error should be about 0.5e-8? Well, the sine function being anti-symmetric and periodic, you need to only approximate it on the interval [0, pi/2]. It turns out you need 42 points, not including the ends of the interval (where we don't need to tabulate the values of sine; it is simply 0 at 0 and 1 at pi/2).

import numpy as np
from scipy.interpolate import CubicSpline
xs = np.linspace(0,np.pi/2,44,endpoint=True)
sines = np.sin(xs)
spline=CubicSpline(xs, sines)
ts = np.linspace(0,1,10000,endpoint=True)\*np.pi/2
print("Max Error is ", np.max(np.abs(spline(ts)-np.sin(ts))))
>> Max Error is 5.02852541828247e-08

But, we can do better. If we tabulate values of sin, then we have for free the values of the cosine at the complementary angles ( cos(x) = sin(pi/2-x). But the derivative of sine is the cosine, so we can use a cubic spline that matches both the values and the derivatives at a number of points. This is called a Cubic Hermite Spline. Turns out we need 23 points.

from scipy.interpolate import CubicHermiteSpline
xs = np.linspace(0,np.pi/2,25,endpoint=True)
sines = np.sin(xs)
cosines = sines[::-1] #note that we reuse the sines,
#we don't invoke np.cos
spline=CubicHermiteSpline(xs, sines, dydx=cosines)
ts = np.linspace(0,1,10000,endpoint=True)*np.pi/2
print("Max Error is ", np.max(np.abs(spline(ts)-np.sin(ts))))

>> Max Error is 4.775709550042251e-08

Now that we got to do Hermite interpolation, why stop here? The second derivative of the sine function is the sine with the sign flipped (oh, man, how do you avoid this pun?). So, you don't need to tabulate any new values. In scipy you can use BPoly to get this higher order spline. Now, you only need to tabulate 4 values of the sine function.

xs = np.linspace(0,np.pi/2,6,endpoint=True)
sines = np.sin(xs)
cosines = sines[::-1]
spline = BPoly.from_derivatives(xs, np.vstack((sines,cosines,-sines)).T)
ts =np.linspace(0,np.pi/2,10000)
print("Max Error is ", np.max(np.abs(spline(ts)-
np.sin(ts))))

>> Max Error is 2.057940906574629e-08

Well, we can go all the way now. If we only use the values of the sine function and its first 4 derivatives only at the end of the [0,pi/2] interval, we can find a 9th degree polynomial that matches them.

spline = BPoly.from_derivatives([0,np.pi/2],
[[0,1,0,-1,0], [1,0,-1,0,1]])
PSpline = PPoly.from_bernstein_basis(spline)
P = Polynomial(PSpline.c.T[0][::-1])
print("Max Error is ", np.max(np.abs(P(ts)-
np.sin(ts))))

>>Max Error is 1.700470619869776e-08

And here's that 9th degree polynomial, if you are curious:

x −0.16666666666666785 x^3+8.881784197001252e-15 x^4+0.008331651428065356 x^5+5.167900831715144e-06 x^6−0.000204626778274708 x^7+3.5525642158307225e-06 x^8+1.8946044087406189e-06 x^9