Table of Contents

  1. Other implementations
  2. Evaluating minimax polynomials
    1. Relative error on float ULPs
    2. Relative error on double precision ULPs
  3. Minimax polynomial for sin(x)
    1. Absolute error
    2. Relative error
    3. Another example

Trigonometric functions

Our research notes about implementation of fast trigonometric functions.

Other implementations

Evaluating minimax polynomials

Relative error on float ULPs

Let’s check how much an error on the last bit affects float variables, using relative error: 

for (float f = 1e-25f; f < 50.0f; f *= 1.001f)
{
    union { float f; uint32_t x; } u = { f }, v = u, w = u;
    ++v.x;
    --w.x;
    printf("%e: %e %e\n", u.f, (v.f - u.f) / u.f, (u.f - w.f) / u.f);
}

Plotting the results gives the answer: between 5.96e-8 and 1.19e-7.

Relative error on double precision ULPs

How about double? 

for (double d = 1e-25; d < 50.0; d *= 1.001)
{   
    union { double d; uint64_t x; } u = { d }, v = u, w = u;
    ++v.x;
    --w.x;
    printf("%e: %e %e\n", u.d, (v.d - u.d) / u.d, (u.d - w.d) / u.d);
}

Result: between 1.11e-16 and 2.22e-16.

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_{y \in [0, \pi^2/4]}{\big\lvert\sin(\sqrt{y}) - \sqrt{y}Q(y)\big\rvert} = E\]

Divide through by √y:

\[\max_{y \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_{y \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_{y \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, 40);

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.666666666666580938362041558393413542600e-1;
a1 = +8.333333333262715528278347301093116699226e-3;
a2 = -1.984126982005911547055378498482154233331e-4;
a3 = +2.755731607338689059680115311170244593349e-6;
a4 = -2.505185130214293461336864464029272975945e-8;
a5 = +1.604729591825977276746033955401272849354e-10;
a6 = -7.364589573262279656101943192883347174447e-13;

E = 1.098969630370672683831702893969063712485e-16;

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_{y \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, 40);

a0 = -1.666666666666618655330686951220600109093e-1;
a1 = +8.333333333285509269011390197274270193963e-3;
a2 = -1.984126982504390943378441670654599796999e-4;
a3 = +2.755731659890484005079622148869365590255e-6;
a4 = -2.505188017067512158000464673183918592305e-8;
a5 = +1.604809231079007402834515514014751626424e-10;
a6 = -7.373308642081174610234470417310065337523e-13;

E = 1.536616934979294924070319645582179175464e-16;

Another example

Trying to reproduce ​Mike Acton's experimental result:

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

Using P(x) = xQ(x²):

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

RemezSolver<4> solver;
solver.Run(real::R_1 >> 400, real::R_1, myfun, myfun, 40);

a0 = +1.570796318452974170937444182514099057668;
a1 = -6.459637106518004316895285560913086847441e-1;
a2 = +7.968967893119537554369879233657335437244e-2;
a3 = -4.673766402124326801818077428151491285886e-3;
a4 = +1.514849803879821510688050931907848092814e-4;

E = 5.310632770140865101487747414605249755889e-9;
Last modified 11 years ago Last modified on Dec 31, 2011, 12:11:00 AM