7 ms·

Implementing Advanced Math Functions

41 points by silentbicycle 15y ago

adrianN 15y ago
This article lacks a discussion about floating-point vs real arithmetic. Naively implementing these functions is bound to produce nasty rounding problems.
- silentbicycle 15y ago
  The previous post in this series is on fixed-point math: http://spin.atomicobject.com/2012/03/15/simple-fixed-point-math/ http://spin.atomicobject.com/2012/03/15/simple-fixed-point-m... .
- jensnockert 15y ago
  The problem isn't in rounding or floating-point, these implementations won't be good enough to get even close to have to consider the issue of having to round correctly. You really need to reduce the argument to a sane range, otherwise the `n' required even in infinite precision would be huge. If someone is interested in how to do it, check out http://www.netlib.org/fdlibm/ http://www.netlib.org/fdlibm/
  - ajross 15y ago
    It's worth pointing out that Netlib, along with glibc, uclibc, bionic and probably the BSDs too, is just using the SunPro code from the early 90's for its trancendental functions. I'm honestly not aware of anyone who isn't. And most of these are doing more than simply reducing N in range. They typically have a transform (or three) which maps N onto a function which is very close to linear, so they can use a 5-6 term polynomial approximation to get full precision across the whole range. It's all very clever, really. But it's terribly documented (and, at this point, ruthlessly forked). I'm honestly shocked there hasn't been an effort made to clean this code up, put it into a single project and get all the divergent libm implementations using it.
davidkellis 15y ago
This is the type of stuff I come to HN for. This is so refreshing to see after days worth of press releases about Apple, new gadgets, and other mainstream news regurgitation.
- zackzackzack 15y ago
  If you want to see some crazy math stuff, look into axiom[0]. It's an open source computer algebra system programmed in Common Lisp. The manuals for the project add up to be ~10,000 pages[1]. It's all done in the literate programming style and it's being actively developed. I am probably going to write an idea up in a month or so, but if a computer programmer wants to learn some advanced Math beyond first year calculus quickly, downloading axiom and working through the first tutorial is the way to go. [0]http://www.axiom-developer.org/ http://www.axiom-developer.org/ [1]http://www.axiom-developer.org/axiom-website/documentation.html http://www.axiom-developer.org/axiom-website/documentation.h...
stephencanon 15y ago
Unfortunately, Taylor series are just about the worst choice possible for implementing transcendental math functions, due to the highly non-uniform error distribution. You end up doing far to much work to get good accuracy near the edges of the domain of interest. Much better choices of polynomial approximations for almost any such situation are (in increasing order of goodness and difficulty to work with): Chebyshev series, Cathedory-Fejer, and Minimax.
- jvranish 15y ago
  If your goal is accuracy, you are absolutely correct. But in the instances where I've wanted to implement my own versions of transcendental math functions, accuracy was actually pretty low on my list of priorities. Often just a couple terms of a taylor series was accurate enough, fast, and easy to implement.
  - stephencanon 15y ago
    It's true regardless of what your goal is (unless your goal is to play with Taylor series approximations). Even if you're going purely for speed, you can often replace a third-order Taylor series with a second-order minimax polynomial, and deliver the same or better accuracy with significantly fewer operations.
    - jvranish 15y ago
      Do you know of good resource for how to compute minimax polynomials that approximate transcendental functions? My Google foo appears to be failing me.
      - jeremysalwen 15y ago
        Look up the remez algorithm
      - stephencanon 15y ago
        The standard approach to minimax approximations is called the "Remez Exchange Algorithm" (or Remez algorithm, or Remes algorithm). It can be fussy, however, and often requires a knowledgeable user to product good results, which is why Chebyshev and Cathedory-Fejer approximants are often used instead -- they are very nearly as good, and much easier to compute.
evincarofautumn 15y ago
I have the feeling that a naive polynomial approximation would give lower relative error than a naive Taylor series approximation—at least for the same number of clock cycles. Or am I just crazy? float sin(float const x) { float const y = (4 / pi) * x + (-4 / (pi * pi)) * x * abs(x); return y + 0.218 * (y * abs(y) - y); }
jeremysalwen 15y ago
I tried to write a library to calculate vectorized trig functions. Minimax or least-squares polynomials are much much better than Taylor series. But here's the thing: The majority of the computation required is not in the calculation of the kernel (i.e. f(x) for -pi/2<x<pi/2) but getting accurate reduction to this range for all possible floating point values. If you're using single precision numbers, you really only need the first 9 terms of the minimax polynomial (=x+x^3+x^5+x^7+x^9) to get within two ulps of accuracy. But doing extremely precise division of large numbers is not so quick. If you're interested, I have the test code here: https://github.com/jeremysalwen/vectrig https://github.com/jeremysalwen/vectrig

41 points by silentbicycle 15y ago

This article lacks a discussion about floating-point vs real arithmetic. Naively implementing these functions is bound to produce nasty rounding problems.

silentbicycle 15y ago
The previous post in this series is on fixed-point math: http://spin.atomicobject.com/2012/03/15/simple-fixed-point-math/ http://spin.atomicobject.com/2012/03/15/simple-fixed-point-m... .
jensnockert 15y ago
The problem isn't in rounding or floating-point, these implementations won't be good enough to get even close to have to consider the issue of having to round correctly. You really need to reduce the argument to a sane range, otherwise the `n' required even in infinite precision would be huge. If someone is interested in how to do it, check out http://www.netlib.org/fdlibm/ http://www.netlib.org/fdlibm/
- ajross 15y ago
  It's worth pointing out that Netlib, along with glibc, uclibc, bionic and probably the BSDs too, is just using the SunPro code from the early 90's for its trancendental functions. I'm honestly not aware of anyone who isn't. And most of these are doing more than simply reducing N in range. They typically have a transform (or three) which maps N onto a function which is very close to linear, so they can use a 5-6 term polynomial approximation to get full precision across the whole range. It's all very clever, really. But it's terribly documented (and, at this point, ruthlessly forked). I'm honestly shocked there hasn't been an effort made to clean this code up, put it into a single project and get all the divergent libm implementations using it.

davidkellis 15y ago

This is the type of stuff I come to HN for. This is so refreshing to see after days worth of press releases about Apple, new gadgets, and other mainstream news regurgitation.

zackzackzack 15y ago
If you want to see some crazy math stuff, look into axiom[0]. It's an open source computer algebra system programmed in Common Lisp. The manuals for the project add up to be ~10,000 pages[1]. It's all done in the literate programming style and it's being actively developed. I am probably going to write an idea up in a month or so, but if a computer programmer wants to learn some advanced Math beyond first year calculus quickly, downloading axiom and working through the first tutorial is the way to go. [0]http://www.axiom-developer.org/ http://www.axiom-developer.org/ [1]http://www.axiom-developer.org/axiom-website/documentation.html http://www.axiom-developer.org/axiom-website/documentation.h...

stephencanon 15y ago

Unfortunately, Taylor series are just about the worst choice possible for implementing transcendental math functions, due to the highly non-uniform error distribution. You end up doing far to much work to get good accuracy near the edges of the domain of interest. Much better choices of polynomial approximations for almost any such situation are (in increasing order of goodness and difficulty to work with): Chebyshev series, Cathedory-Fejer, and Minimax.

jvranish 15y ago
If your goal is accuracy, you are absolutely correct. But in the instances where I've wanted to implement my own versions of transcendental math functions, accuracy was actually pretty low on my list of priorities. Often just a couple terms of a taylor series was accurate enough, fast, and easy to implement.
- stephencanon 15y ago
  It's true regardless of what your goal is (unless your goal is to play with Taylor series approximations). Even if you're going purely for speed, you can often replace a third-order Taylor series with a second-order minimax polynomial, and deliver the same or better accuracy with significantly fewer operations.
  - jvranish 15y ago
    Do you know of good resource for how to compute minimax polynomials that approximate transcendental functions? My Google foo appears to be failing me.
    - jeremysalwen 15y ago
      Look up the remez algorithm
    - stephencanon 15y ago
      The standard approach to minimax approximations is called the "Remez Exchange Algorithm" (or Remez algorithm, or Remes algorithm). It can be fussy, however, and often requires a knowledgeable user to product good results, which is why Chebyshev and Cathedory-Fejer approximants are often used instead -- they are very nearly as good, and much easier to compute.

evincarofautumn 15y ago

I have the feeling that a naive polynomial approximation would give lower relative error than a naive Taylor series approximation—at least for the same number of clock cycles. Or am I just crazy? float sin(float const x) { float const y = (4 / pi) * x + (-4 / (pi * pi)) * x * abs(x); return y + 0.218 * (y * abs(y) - y); }

jeremysalwen 15y ago

I tried to write a library to calculate vectorized trig functions. Minimax or least-squares polynomials are much much better than Taylor series. But here's the thing: The majority of the computation required is not in the calculation of the kernel (i.e. f(x) for -pi/2<x<pi/2) but getting accurate reduction to this range for all possible floating point values. If you're using single precision numbers, you really only need the first 9 terms of the minimax polynomial (=x+x^3+x^5+x^7+x^9) to get within two ulps of accuracy. But doing extremely precise division of large numbers is not so quick. If you're interested, I have the test code here: https://github.com/jeremysalwen/vectrig https://github.com/jeremysalwen/vectrig