Sin and Cos: The Programmer's Pals
I'm half way through writing a blog post on how transcendental functions are computed in glibc, and more importantly how to extend the methods to compute things like  sin(x) / x
  (1 - cos(x)) / x^2

directly without incurring the numerical problems arising from the division.One thing to understand is that the stdlib math functions are accurate to machine precision, i.e. the closest representable value to the actual answer.For 32-bit floats, this is achievable using a 5 or 6 order polynomial approximation. However, there was one surprise I was previously unaware of.Obtaining the correct precision becomes difficult when the result is close to zero, since any absolute error is divided by the magnitude of the result. Silly example, if  sin(0.000001) ~= 0.1

Then the absolute error is small (~0.1), but the result is still off by four orders of magnitude.Sometimes this is not a problem. When x is small, sin(x) can be computed by simply returning x. This rule is valid until x^3/6 > precision ~= 1e-38, or x = 4e-16. So exact precision is obtainable by simply returning x if x But now consider:  x = pi/2; cos(x) = 0

The logical way to compute this is as  -sin(x - pi/2)

But this creates another problem, since x - pi/2 ~= 0, a catastrophic cancellation occurs. Now recall that glibc is required to be accurate to machine precision, which goes down to 1e-38. But x is order 1 and only represents 7 digits, down to 1e-6. So glibc is forced to scour lookup tables to recover the missing digits of pi.This makes computing cos(x > pi/4) about twice as expensive than cos(x a + bb + c*c); ?On my 4 year old laptop:  time ./a.out 0.5   4.12s
  time ./a.out 1.5   7.31s

On my raspberry pi:  time ./a.out 0.5   1m15.829s
  time ./a.out 1.5   1m47.522s