NavList:

   

Reply

In the old days applications didn't always run on a computer able to
calculate sines etc. in hardware. Remember when a deluxe machine had a
387 on the motherboard? So, if you called sin() the compiler would
include a routine which used the primitive operators of the language (+,
-, etc.) to calculate the function. A book by P.J. Plaugher ("The
Standard C Library," 1992) provides C language source code for a full
set of math functions. I thought it would be interesting to look at his
sine and cosine function.

The first operation is subtract a whole multiple of π/2 (90°) to reduce
argument x (radians) to the interval [-π/4, +π/4], i.e., -45° to +45°.
The whole multiple is the 90° quadrant in which x lies. Quadrant 0 is
-45 to +45, quadrant 1 is 45 to 135, etc. Quadrant number is computed by
multiplying x by 2/π (in degrees, equivalent to division by 90). Round
the result to the nearest whole number.

However, the obvious method to remove whole multiples of π/2 loses
accuracy, especially if x is large. It's easier to see this if we
imagine a large x and a calculator with only 6 digits. Suppose x =
899.188 radians. That angle is in quadrant 572, and (π/2 * 572) =
898.498. When you subtract that from x, the result is 0.690. The last
digit is wrong by 2.5 units.

A more accurate result is possible with our 6-digit calculator if we
split π/2 into components c1 = 1.57 and c2 = 7.96327e-4. Their sum is
π/2. To remove 572 quadrants from x:

x = 899.188 - 572 * c1 - 572 * c2
x = 899.188 - 898.040 - 0.455499
x = 0.692501

That's more accurate. It could be argued that 6 decimal places is fake
precision since x is known to only .001. Well, we can't improve the
precision of the argument, but we can try to extract the maximum from
what has been supplied.

The computer implementation by Plaugher uses the same principle, but
with the 64-bit double precision floating point format. Constant c1 is
defined in a non-obvious way:

c1 = 3294198.0 / 2097152.0

The numerator is an integer of 22 bits. The denominator is 2 to the 21st
power. Thus, to form the constant c1, simply shift the binary point of
the numerator 21 places to the left. That gives the first 22 bits of
π/2. The other bits are 0.

If you try to initialize c1 with a numeric literal (1.570...) that exact
bit pattern is not guaranteed since a decimal fraction doesn't have an
exact binary representation (except special cases such as 0.25).

Constant c2 = 3.139164786504813217e-7. It supplies the remaining digits
of π/2. I.e., c1 + c2 = π/2 to high precision. Then Plaugher does this
to reduce angle x to the interval [-π/4, +π/4]:

x = (x - quadrant * c1) - quadrant * c2

He doesn't explain the parentheses, but I think their purpose is to
force the order of evaluation. If you change the order of the
subtractions in my 6-digit example there's a loss of accuracy.

But parentheses do not guarantee order of evaluation in the C language.
As explained by Harbison and Steele ("C, A Reference Manual," 1991), a
compiler is allowed to disregard parentheses and re-arrange expressions
to a mathematically equivalent sequence. For example, if your code says
(a + b) + (c + d) it might execute, right to left, as a + d + b + c.
They say the only way to be sure of evaluation order is use temporary
variables. Thus, I added a double precision temporary variable to the
Plaugher code:

double temp = x - quadrant * c1;
x = temp - quadrant * c2;

I have explained the adjustment of the argument at length to make clear
why it's a bad idea to remove whole circles to "help" a trig function.
It knows how to do that, and probably does a better job than you.

In the second part I'll show Plaugher's sine and cosine expressions.

--
Paul Hirose
sofajpl.com

   

Reply