NavList:

   

Reply

Earlier I explained the algorithm employed by P.J. Plaugher in his book
"The Standard C Library" (1992) to normalize the argument (radians) of
his sine / cosine function to the interval [-π/4, +π/4], i.e., -45° to
+45°. All that remains is to evaluate either of two polynomials, and
finally to reverse the sign if appropriate. The steps depend on the
quadrant of argument x, calculated in the previous message, and whether
the call of the function requested a sine or cosine.

That last is indicated by integer parameter which is 0 for sine, 1 for
cosine, and is added to the quadrant number. The result is an integer
whose last two bits control the switching logic. But that's not reliable
if quadrant is negative because a C language implementation is allowed
to represent negative integers in whatever format is natural for the
hardware. It may be ones complement or twos complement. That means the
bit pattern isn't predictable.

The solution is a variable of type "unsigned integer". The conversion
from a signed to unsigned integer is predictable regardless of machine
architecture. So you first compute quadrant in double precision floating
point, round that to the nearest signed integer, convert to an unsigned
integer (which I call "selector"), and add the cosine indicator.

Plaugher provides a polynomial for sine and another for cosine, valid in
the range [-π/4, π/4] radians, i.e., -45° to +45°. Each can be used for
the other function. For instance (in degrees), suppose x = 60 and sine
is requested. After subtracting whole multiples of 90 to put it in the
interval [-45, +45], x = -30. Then evaluate the cosine polynomial. Since
cos(-30) = sin(60), that gives the sine.

The polynomials are 14th degree, but only the even powers are used
(except one 1st degree term in the sine poly). So begin by initializing
double precision variable x2 to the square of x:

double x2 = x * x;

Then the last bit of the selector variable is tested to choose the
correct poly. (The C language & operator performs a bit by bit AND of
its operands.)

if ((selector & 1) == 0)  // use sine polynomial

x = (x2 * (x2 * (x2 * (x2 * (x2 * (x2 * (x2 *
-0.000000000000764723
+ 0.000000000160592578)
- 0.000000025052108383)
+ 0.000002755731921890)
- 0.000198412698412699)
+ 0.008333333333333372)
- 0.166666666666666667)
+ 1.0) * x;

else    // use cosine polynomial

x = (x2 * (x2 * (x2 * (x2 * (x2 * (x2 * (x2 *
-0.000000000011470879
+ 0.000000002087712071)
- 0.000000275573192202)
+ 0.000024801587292937)
- 0.001388888888888893)
+ 0.041666666666667325)
- 0.500000000000000000)
+ 1.0);

The polynomials may look unorthodox because they've been factored as
much as possible. For instance, consider the third degree expression

ax^3 + bx^2 + cx + d

Factor out x from the first three terms:

x(ax^2 + bx + c) + d

Factor out x again from the first two terms:

x(x(ax + b) + c) + d

Indicate multiplication explicitly with the * symbol and put the
coefficients on separate lines:

x * (x * (x *
a
+ b)
+ c)
+ d

One coefficient per line makes it easier to input and check the values.

After evaluation of the sin or cos polynomial it may be necessary to
reverse the sign of the result. That's determined by testing the 2nd bit
from last in the selector variable.

if ((selector & 2) == 0)
return x;
else
return -x;

I evaluated sine and cosine at 10 random angles. Root of the mean
squared error was 1.0 counts in the 16th decimal place. The (assumed)
correct values came from the Windows calculator accessory, which is
accurate to 32 decimal places.

What I presented was a simplified version of the Plaugher code. It omits
his special handling for extremely large or small values, as well as
"not a number" arguments.

I think the main thing we can take away from all this is gratitude that
we don't have to code our own transcendental functions! It requires
mathematical facility plus knowledge of obscure corners of your
programming language.

--
Paul Hirose
sofajpl.com

   

Reply