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

--

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²), and we reduce the range to positive values:

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

Substitute y for x²:

\[\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;

Relative error

Searching for relative error instead:

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

Using the same method as for absolute error, we get:

\[\max_{x \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert}{\bigg\lvert\dfrac{\sin(y)}{y\sqrt{y}}\bigg\rvert}} = E\] 
static real myfun(real const &y)
{
    real x = sqrt(y);
    return (sin(x) - x) / (x * y);
}

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

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

a0 = -1.666666666666666587374325845020415990185e-1;
a1 = +8.333333333333133768001243698120735518527e-3;
a2 = -1.984126984109960366729319073763957206143e-4;
a3 = +2.755731915499171528179303925040423384803e-6;
a4 = -2.505209340355388148617179634180834358690e-8;
a5 = +1.605725287696319345779134635418774782711e-10;
a6 = -7.535968124281960435283756562793611388136e-13;