Version 12 (modified by sam, 9 years ago) (diff)

--

# Trigonometric functions

Our research notes about implementation of fast trigonometric functions.

## Evaluating minimax polynomials

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.

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:

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:

Substitute y for x²:

Divide through by √y:

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

Divide through by y:

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:

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

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

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:

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

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;