4 ms·

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

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

yes_or_gnome 9y ago
> 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. "This rule" seems a lot like the "skinny triangle" rule. I recently ran across it on Wikipedia after going through astrophysics articles. Specifically, to my understanding, the parsec could be calculated without using a trig function thanks to the skinny triangle rule. (To a certain level of error, of course. And, the parsec is now a defined value so the original trig-derived definition is invalid.) https://en.wikipedia.org/wiki/Skinny_triangle https://en.wikipedia.org/wiki/Skinny_triangle
- Sean1708 9y ago
  > "This rule" seems a lot like the "skinny triangle" rule This rule is exactly the Skinny Triangle Rule, or rather the part of the Skinny Triangle Rule which overlaps with the Small Angle Approximation[0]. [0]: https://en.wikipedia.org/wiki/Small-angle_approximation https://en.wikipedia.org/wiki/Small-angle_approximation
vanderZwan 9y ago
Do you have advice for which fast, good-enough approximations one can use instead in the situations where the precision is not needed? In many use-cases we would pre-multiply the value we throw into sin and cos with two times Pi anyway, so couldn't one skip the circle constant altogether and a write function like sin_one(t) and cos_one(t), where each call "implicitly multiplies" t by two times Pi? (perhaps, if constructed well, the function could even be more precise, since it avoids the rounding errors, however small, of multiplying by the pi constant)
- 0xfaded 9y ago
  In most cases a Chebyshev approximation will just work. There is advice out there which says to use Clenshaw's summation formula to sum the Chebyshev polynomials, however I've found that expanding out the polynomial coefficients to get something like a5 x^5 + a4 x^4 + a3 x^3 to be more accurate in practice. Scaling the polynomial, as you suggest, would simply require reevaluating the coefficients. For example, a5 (2pix)^5 would become 32 a5 pi^5 x^5. In the cos case where regions of interest approach zero, divide by another polynomial with the same zeros to create limits. cos(x)/(pi^2/4 - x^2) -> 1/pi as x -> pi/2. Then approximate that function, and evaluate approx(x) \* (pi^2/4 - x^2), so you are multiplying the approximated value by something close to zero. This way you get precision in a way similar to sin(x) = x when x small.
- stochastic_monk 9y ago
  I don't have advice on the specific `sin(x) / x`/`cos(x) / x^2`, but the Sleef vectorization library [0,1] supports operations with varying precision. (0.506, 1.0, and 3.5 ULP) If you're working in a case where you can afford slightly lower accuracy, the 3.5ULP method can be a very significant speedup, especially since you're vectorizing at the same time. I got dozens of times in speedup when testing just the 1.0 ULP bound when performing sin/cos on a large vector. (I don't remember the precise number and should just re-run the tests.) To generically dispatch this library, in case there's any interest, I wrapped it in some template metaprogramming in [2], so that the widest vectorization available is selected at compile-time when selecting a given operation. [0]: http://sleef.org/ http://sleef.org/ [1]: https://github.com/shibatch/sleef https://github.com/shibatch/sleef [2]: https://github.com/dnbaker/vec https://github.com/dnbaker/vec
jwilk 9y ago
> glibc is required to be accurate to machine precision Is this actually the case? This patch seems to suggest it's not: https://www.sourceware.org/ml/libc-alpha/2018-02/msg00336.html https://www.sourceware.org/ml/libc-alpha/2018-02/msg00336.ht...
- 0xfaded 9y ago
  It seems you are right. IEEE doesn't say anything, which I guess makes it implementation defined. However, in glibc, there are definitely code paths for looking up digits of pi.
deleted 9y ago
[deleted]

yes_or_gnome 9y ago

> 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. "This rule" seems a lot like the "skinny triangle" rule. I recently ran across it on Wikipedia after going through astrophysics articles. Specifically, to my understanding, the parsec could be calculated without using a trig function thanks to the skinny triangle rule. (To a certain level of error, of course. And, the parsec is now a defined value so the original trig-derived definition is invalid.) https://en.wikipedia.org/wiki/Skinny_triangle https://en.wikipedia.org/wiki/Skinny_triangle

Sean1708 9y ago
> "This rule" seems a lot like the "skinny triangle" rule This rule is exactly the Skinny Triangle Rule, or rather the part of the Skinny Triangle Rule which overlaps with the Small Angle Approximation[0]. [0]: https://en.wikipedia.org/wiki/Small-angle_approximation https://en.wikipedia.org/wiki/Small-angle_approximation

vanderZwan 9y ago

Do you have advice for which fast, good-enough approximations one can use instead in the situations where the precision is not needed? In many use-cases we would pre-multiply the value we throw into sin and cos with two times Pi anyway, so couldn't one skip the circle constant altogether and a write function like sin_one(t) and cos_one(t), where each call "implicitly multiplies" t by two times Pi? (perhaps, if constructed well, the function could even be more precise, since it avoids the rounding errors, however small, of multiplying by the pi constant)

0xfaded 9y ago
In most cases a Chebyshev approximation will just work. There is advice out there which says to use Clenshaw's summation formula to sum the Chebyshev polynomials, however I've found that expanding out the polynomial coefficients to get something like a5 x^5 + a4 x^4 + a3 x^3 to be more accurate in practice. Scaling the polynomial, as you suggest, would simply require reevaluating the coefficients. For example, a5 (2pix)^5 would become 32 a5 pi^5 x^5. In the cos case where regions of interest approach zero, divide by another polynomial with the same zeros to create limits. cos(x)/(pi^2/4 - x^2) -> 1/pi as x -> pi/2. Then approximate that function, and evaluate approx(x) \* (pi^2/4 - x^2), so you are multiplying the approximated value by something close to zero. This way you get precision in a way similar to sin(x) = x when x small.
stochastic_monk 9y ago
I don't have advice on the specific `sin(x) / x`/`cos(x) / x^2`, but the Sleef vectorization library [0,1] supports operations with varying precision. (0.506, 1.0, and 3.5 ULP) If you're working in a case where you can afford slightly lower accuracy, the 3.5ULP method can be a very significant speedup, especially since you're vectorizing at the same time. I got dozens of times in speedup when testing just the 1.0 ULP bound when performing sin/cos on a large vector. (I don't remember the precise number and should just re-run the tests.) To generically dispatch this library, in case there's any interest, I wrapped it in some template metaprogramming in [2], so that the widest vectorization available is selected at compile-time when selecting a given operation. [0]: http://sleef.org/ http://sleef.org/ [1]: https://github.com/shibatch/sleef https://github.com/shibatch/sleef [2]: https://github.com/dnbaker/vec https://github.com/dnbaker/vec

jwilk 9y ago

> glibc is required to be accurate to machine precision Is this actually the case? This patch seems to suggest it's not: https://www.sourceware.org/ml/libc-alpha/2018-02/msg00336.html https://www.sourceware.org/ml/libc-alpha/2018-02/msg00336.ht...

0xfaded 9y ago
It seems you are right. IEEE doesn't say anything, which I guess makes it implementation defined. However, in glibc, there are definitely code paths for looking up digits of pi.

deleted 9y ago

[deleted]