Version 6 (modified by sam, 11 years ago) (diff)

use \max instead of \forall

Trigonometric functions

Our research notes about implementation of fast trigonometric functions.

Minimax polynomial for sin(x)

Absolute error

Suppose we want to approximate sin(x) on [-π/2; π/2] with a polynomial P(x) such that the absolute error is never more than E:

\[\max_{x \in [-\pi/2, \pi/2]}{\big\vert\sin(x) - P(x)\big\vert} = E\]

We know sin(x) is an odd function, so instead we look for a polynomial Q(x) such that P(x) = xQ(x²):

\[\max_{x \in [-\pi/2, \pi/2]}{\big\vert\sin(x) - xQ(x^2)\big\vert} = E\]

Substitute y for x² and reduce the range to positive values:

\[\max_{x \in [0, \pi^2/4]}{\big\lvert\sin(\sqrt{y}) - \sqrt{y}Q(y)\big\rvert} = E\]

Divide through by √y:

\[\max_{x \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})}{\sqrt{y}} - Q(y)\bigg\rvert}{\dfrac{1}{|\sqrt{y}|}}} = E\]

If we want to force the asymptotic behaviour at x=0, we substitute Q(y) with 1+yR(y):

\[\max_{x \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})}{\sqrt{y}} - 1 - yR(y)\bigg\rvert}{\dfrac{1}{|\sqrt{y}|}}} = E\]

Divide through by y:

\[\max_{x \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert}{\dfrac{1}{|y\sqrt{y}|}}} = E\]

We then use the following code: 

static real myfun(real const &y)
{
    real x = sqrt(y);
    return (sin(x) - x) / (x * y);
}

static real myerr(real const &y)
{
    return re(y * sqrt(y));
}

RemezSolver<6> solver;
solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 15);

These are the resulting R, Q and P:

\begin{eqnarray*} R(y) & = & a_0 + a_1 y + a_2 y^2 + a_3 y^3 + a_4 y^4 + a_5 y^5 + a_6 y^6 \\ Q(y) & = & 1 + a_0 y + a_1 y^2 + a_2 y^3 + a_3 y^4 + a_4 y^5 + a_5 y^6 + a_6 y^7 \\ P(x) & = & x + a_0 x^3 + a_1 x^5 + a_2 x^7 + a_3 x^9 + a_4 x^{11} + a_5 x^{13} + a_6 x^{15} \\ \end{eqnarray*}

With the following coefficients: 

a0 = -1.66666666666658080941942898789420724e-1;
a1 = +8.33333333326271609442503773834687308e-3;
a2 = -1.98412698200591143928364634696492885e-4;
a3 = +2.75573160733868922065738227278330896e-6;
a4 = -2.50518513021429359590028300127165228e-8;
a5 = +1.60472959182597740337401201006549498e-10;
a6 = -7.36458957326227991327065122848667046e-13;