Numerically Efficient Approximation of cos(s)

Question

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.

Answer 1

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: https://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.

Answer 2

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.

Answer 3

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++];
  (* Padé approximant *)
  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.

Answer 4

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

Answer 5

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++];
  (* Padé approximant *)
  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]
               }}]