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<p><pre><code>  sin(x) / x
  (1 - cos(x)) / x^2
</code></pre>
directly without incurring the numerical problems arising from the division.<p>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.<p>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.<p>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<p><pre><code>  sin(0.000001) ~= 0.1
</code></pre>
Then the absolute error is small (~0.1), but the result is still off by four orders of magnitude.<p>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 < 1e-16.<p>But now consider:<p><pre><code>  x = pi/2; cos(x) = 0
</code></pre>
The logical way to compute this is as<p><pre><code>  -sin(x - pi/2)
</code></pre>
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.<p>This makes computing cos(x > pi/4) about twice as expensive than cos(x < pi/4), but there is a certain irony since it is very unlikely x is known down to 7+ digits, so those CPU cycles are wasted. How many times have you seen x = sqrt(a<i>a + b</i>b + c*c); ?<p>On my 4 year old laptop:<p><pre><code>  time ./a.out 0.5   4.12s
  time ./a.out 1.5   7.31s
</code></pre>
On my raspberry pi:<p><pre><code>  time ./a.out 0.5   1m15.829s
  time ./a.out 1.5   1m47.522s</code></pre>