Implementing cosine in C from scratch
web.eecs.utk.edu
web.eecs.utk.edu
Actually minimizing error in the floating point domain is a really cool subject. Intuitively, you want the most accuracy when you're close to zero and the mantissa is 10e-37. When you're in the region with mantissa 1, +- 10e-10 means a lot less.
One "trick" is to approximate a different polynomial, and then multiply by x, since x itself is naturally distributed across the floating point domain.
E.g, instead of approximating sin(x), approximate sin(x)/x distributing the error evenly between +-pi/4. Then multiply by x to get sin(x).
For cos, cos(x) gets small and difficult at +-pi/2. If you look at libc, you'll find lots of if statements that apply different hacks to preserve accuracy depending on how close to zero the result is.
E. g., you can get your approximation both smooth and continuous by adding the target derivatives in +-PI/2 into the Vandermonde matrix. The overall error may become higher with the same polynomial order but the approximation will be smooth.
I measured this for a contemporary Intel processor some maybe ten years ago and found out that 6th degree (actually, its 3 for x^2 since you can throw away all the even coefficients) makes sense performance-wise, yes the error is greater but that's the point of the trade-off. But even 8 or 10 is already doubtful.
Of course, there are architectures other than Intel and it might make sense to use higher order polynomials for better accuracy there.
If I remember correctly, Pade approximants would be more accurate with less work. Instead of using a polynomial, you use the ratio of two polynomials. They generally work better than Chebyshev since they're a larger class of approximants.
So, instead of dividing two polynomials of degree n, you could compute one of degree 2n.
I wrote with "Pade approximants would be more accurate with less work." for a reason.
If they have identical complexity they cannot have twice the degrees of freedom. You're confusing ideas here.
First, there are a lot of ways to evaluate Pade than a naive compute two polynomials and then divide them. There's quite a large literature on evaluation methods, just as there are for polynomails (and the naive polynomial evaluation method is generally terrible - at least consider Horner's method).
Second, the Chebyshev poly of deg d can simply be the Pade poly P(n) / 1 with P(n) of degree d. Thus Checbyshev is a special case of Pade at the same cost. Thus it's a smaller class, in the pure mathematical sense that one space of functions is a strict subset of the other.
So if you want to check all approximations of fixed computational complexity, vary the top and bottom degrees appropriately to check a strictly larger space of approximants for fixed computational cost. You'll usually find that Pade, for the same computational cost, makes a better approximant. It's simply richer.
What's more, once you have computed the various powers of x, you can reuse them (unless you're using Horner's method or similar, which is generally better).
If you're really looking for minmal cost computation, there is another entire literature on decomposition of functions to minimal operations, and finding the optimal one is NP-complete (but completely doable for these low order approximations). Knuth TAOCP has an intro to the lit under (if I recall) "addition chains" for evaluating polynomials. There is parallel theory for Pade.
Thus Pade is a strictly better method. It completely encompasses Chebyshev, and since I've used it a lot, I have found it always better, for the same computational work.
It simply let's you optimize over a larger class of approximants.
This works fast, doesn't use the cache, and performs reasonably well precision-wise. It also has very small error in +/-PI/2 so, unlike with the Taylor's series, the continuation is virtually seamless.
Moreover, since cosine is symmetrical, you can throw away the even coefficients completely and it will become even faster.
You may want to optimize for something else than minimizing the maximum error, though, as it may have glaring errors. For example, for sin(x), it feels worse to have sin(0) not be exactly equal to zero than to have sin(π/4) not be exactly equal to ½√2.
It may be hard to come up with a good better optimization criterion, though, and it may be easier to avoid ‘howler’ errors by making your approximation good enough.
[1] https://news.ycombinator.com/item?id=15266331
Similar code with optional double precision support and more options: https://www.geometrictools.com/GTE/Mathematics/CosEstimate.h
If you don’t like C++, both of them are fairly simple to port to C.
I didn’t benchmark it though. I’m curious how it would do.