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 < 1e-16.
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 < 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(aa + 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
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
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 < 1e-16.
But now consider:
The logical way to compute this is as 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 < 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(aa + bb + c*c); ?
On my 4 year old laptop:
On my raspberry pi: