Re: Minimax approximation for cosine

From: Glen Low (glenlow_at_pixelglow.com)
Date: 11/17/04 

Date: 17 Nov 2004 05:44:19 -0800

> For cos(x) on the interval [0, Pi/2], the following is a minimax
> approximation based on absolute error. The approximation is accurate
> to 24.35 bits, so the result can be off by one in the last place.
>
> P0 = 0.99999 99534 64
> P2 = -0.49999 90534 55
> P4 = 0.04166 35846 769
> P6 = -0.00138 53704 264
> P8 = 0.00002 31539 3167
>
> cos(x) = (((P8 * x^2) + P6) * x^2) + P4) * x^2) + P2) * x^2 + P0

Hmm... because it's based on absolute error it fares badly near Pi/2.
E.g. in Mathematica, using x=Pi/2-0.001, the relative error is
0.0000436172.

> Sometimes, it helps to break the P0 term into two halves and add them
> separately.
>
> Notice that cos(0) is not 1.

If it has to be, it has to be. The main consideration is even powers
only.

> I'll type in the coefficients for splitting the range [0, Pi/2] into
> two parts, but let me ask a couple of questions first:
>
> 1. Are you developing a "library quality" routine or a special purpose
> one?

A library quality implementation for SIMD hardware. It should not
differ by more than 5 ulps at most from a scalar library
implementation.

> 2. Is the range [0, Pi/2] the result of range reduction?

Yes it is. I've managed to do a good trigonometric range reduction for
sine where 1 ulp <= Pi, basically using 72 bits of Pi for IEEE 754
float (23 bits). In fact my first attempt at cos was just to throw in
a bias of 0.5 so that the range reduction would do the right thing but
alas the errors were pretty awful around Pi/2.

Cheers,
Glen Low, Pixelglow Software
www.pixelglow.com