Numerically Efficient Approximation of cos(s)

I have an application where I need to run $\cos(s)$ (and $\operatorname{sinc}(s) = \sin(s)/s$) a large number of times and is measured to be a bottleneck in my application.

I don't need every last digit of accuracy: $10^{-10}$ over an input range of $[-15^\circ < s < 15^{\circ}]$ should be sufficient (This is the limit of the input sensor data).

I have implemented a simple Taylor approximation, but would like to ask:

• (a) Is there a more efficient approximation?
• (b) If Taylor is the most efficient, is there a better way to implement it?

// Equation Coefficients
double s_2  = 0.25 * omega_ebe.squaredNorm();
double s_4  = s_2 * s_2;
double s_6  = s_4 * s_2;
double s_8  = s_6 * s_2;
double s_10 = s_8 * s_2;

double cos_coef  = 1.0 - s_2 /  2.0 + s_4 /  24.0 - s_6 /   720.0 + s_8 /  40320.0 - s_10 /  3628800.0;
double sinc_coef = 0.5 - s_2 / 12.0 + s_4 / 240.0 - s_6 / 10080.0 + s_8 / 725760.0 - s_10 / 79833600.0;


EDIT: I haven't forgotten to select an answer! I'm going to code a few of the up and run them on target (an embedded PowerPC and an embedded ARM) to see how they perform.

-
It's been a while for me, but perhaps you might want to look at the CORDIC algorithm? – user4143 Jan 9 '12 at 7:52
What about pre-calculating a (sufficiently large) table and using linear interpolation? – Potato Jan 9 '12 at 7:53
I suggest using the small angle approximations (see here: en.wikipedia.org/wiki/Small-angle_approximation) for $|s| < 10$ and using the Taylor approximation for $|s| \geq 10$. Potato's idea is also a very good one, which you may want to implement for the case $|s| \geq 10$. – user18063 Jan 9 '12 at 7:59
Of course, the small angle approximation is Taylor's approximation with a small number of terms. You would essentially increase the number of terms once you get passed the threshold $|s| = 10$. – user18063 Jan 9 '12 at 8:07
You can compute $a+bx+cx^2+dx^3$ as $((dx+c)x+b)x+a$, which avoids precomputing powers at all giving you fewer multiplications/divisions. Other than that, Taylor seems OK in the range requested. – fedja Jan 9 '12 at 12:31

This is probably more a question for StackOverflow because how efficient algorithms are depends on the processor architecture, you have to value tables against direct calculations and programmers are really good at that.

It is not surprising, then, that there is already a question that handles your problem quite well: http://stackoverflow.com/questions/2088194/fast-sin-cos-using-a-pre-computed-translation-array. Most solutions don't yield their accuracy but you can check it afterwards and raise the accuracy if it is not high enough.

-
Even better, I'd now ask the question on the new Computational Science page :) – Damien Jan 18 '12 at 12:30
The referenced solution has this absolute gem of a paper. – Damien Jan 19 '12 at 23:18
@Damien this is a real gem, thank you! – faleichik Feb 5 '12 at 9:18
I have accepted this answer, largely because of the excellent solution referenced. Thank you to the other contributors (especially @Peter Taylor) who have also given excellent answers. – Damien Feb 5 '12 at 12:49

You can improve on the Taylor/Padé idea. What you can do is scale your argument by repeated halving; say, $z^\prime=\dfrac{z}{2^n}$, where $n$ is chosen such that the error term in either of the Taylor or Padé approximants is tiny enough. Having done so, you can iterate the double-angle identity (that is, $\cos\,2z=2\cos^2 z-1$) $n$ times on the result of the Taylor or Padé approximant evaluated at $z^\prime$. The result is the cosine value you need.

If you need a sine value as well, you will have to couple the double angle identity for the cosine with the double angle identity for the sine, $\sin\,2z=2\sin\,z\cos\,z$.

A Mathematica demonstration:

testCos[z_?InexactNumberQ] := Module[{co, za = z, k = 0},
(* repeated halving *)
While[Abs[za] > 1/4, za /= 2; k++];
co = (1 + za^2*(-115/252 + (313*za^2)/15120))/(
1 + za^2*(11/252 + (13*za^2)/15120));
(* double-angle identity *)
Do[co = 2 co^2 - 1, {k}];
co]

Plot[testCos[z] - Cos[z], {z, -Pi, Pi}, PlotRange -> All]


The Abs[za] > 1/4 criterion was determined by plotting the Padé approximant and determining the range where the difference between it and the cosine was "small enough". You will have to experiment with your particular application.

fedja's advice in the comments should be emphasized quite a bit: whether you take the Taylor or Padé route, it pays to put your polynomial or rational function in Horner form. Less effort is needed in the evaluation, and it can be more accurate than the expanded polynomial or rational function.

-
I like this method, since it generalizes nicely to evaluating the matrix-valued cosine. I have also used a similar route for numerically evaluating Jacobi and Weierstrass elliptic functions. – J. M. Jan 9 '12 at 14:06
(+1). Also the error plot looks like it is easy to approximate, that could be used to refine it even more. – Listing Jan 16 '12 at 14:57

Padé approximations are better than Taylor approximations for the same number of terms, but because you probably care more about average or maximum error over a range of interest rather than around a specific point, in practice different polynomials or rational functions are used. They're often derived using the Remez exchange algorithm.

For examples of functions used in practice see fdlibm or Boost Math.

I've also managed to dig out a paper I remembered reading on using genetic algorithms to evolve improvements to Padé approximations: Towards a Better Sine Wave, Matthew Streeter and Lee A. Becker, 2001 Genetic and Evolutionary Computation Conference Late Breaking Papers, pp398-404.

-

Here's a very important question: what is the range of your inputs? If they are of considerable size (say, ~2^100) then accurately taking modulo 2*pi (to reduce the range to where your Taylor series is accurate) will dominate the computational effort.

As for the kernel function, I would second the Remez algorithm for generating polynomials.

An easier method would be to generate the least squares polynomial for the interval you're interested in. It's basic linear algebra: you're projecting your vector (cosine) onto a set of basis vectors (the polynomials x^0...x^n), with the inner product(A,B) defined as the definite integral of A*B over the interval you're interested in. This will get you the polynomial with the least total square error from the cosine function of a given degree.

If you're interested in looking at some code, I have a few examples of this sort of thing for the sine function: https://github.com/jeremysalwen/vectrig

-

Since it seems that you need both the cosine and sine (cardinal) in your application, another method you can use is based on the identities

\begin{align*}\cos\,x&=\frac{1-\tan^2\frac{x}{2}}{1+\tan^2\frac{x}{2}}\\\sin\,x&=\frac{2\tan\frac{x}{2}}{1+\tan^2\frac{x}{2}}\end{align*}

and the double angle identity

$$\tan\,2x=\frac{2\tan\,x}{1-\tan^2 x}$$

The method is similar to the cosine method I gave earlier; halve the argument an appropriate number of times such that the error of the Taylor or Padé approximant evaluated at that reduced argument is small, evaluate the approximant, repeatedly apply the double angle formula, and transform into the sine and cosine using the appropriate identities.

Here is a Mathematica demonstration:

testCS[z_?InexactNumberQ] := Module[{ta, za = z/2, k = 0},
(* repeated halving *)
While[Abs[za] > 1/4, za /= 2; k++];
ta = za (1 + (za^2/105 - 1) (za/3)^2)/(1 + (za^2/7 - 4) (za/3)^2);
(* double-angle identity *)
Do[ta = 2 ta/(1 - ta^2), {k}];
(* cosine and sine *)
{(1 - ta^2)/(1 + ta^2), 2 ta/(1 + ta^2)}]

GraphicsGrid[{{
Plot[First[testCS[z]] - Cos[z], {z, -Pi/6, Pi/6},
PlotRange -> All],
Plot[Last[testCS[z]] - Sin[z], {z, -Pi/6, Pi/6},
PlotRange -> All]
}}]


-
While I like the theory and will have a play with it, it really need $sinc(s) = sin(s)/s$ (and trying to avoid divisions near zero). – Damien Jan 18 '12 at 12:33
@Damien: some modification would be needed for sine cardinal. Let me think about it... – J. M. Jan 18 '12 at 12:35