Implementing Advanced Math Functions • Hacker News Board

Re: Implementing Advanced Math Functions

#21

Post by stephencanon » Wed, Apr 25, 2012, 5:31 PM UTC

jvranish wrote: post #17

stephencanon wrote: post #10

Earlier quoted context omitted.

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.

Do you know of good resource for how to compute minimax polynomials that approximate transcendental functions? My Google foo appears to be failing me.

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.

Re: Implementing Advanced Math Functions

#22

Post by ArbitraryLimits » Wed, Apr 25, 2012, 6:38 PM UTC

stephencanon wrote: post #16

ArbitraryLimits wrote: post #15

Earlier quoted context omitted.

I thought the truncated Chebyshev series of an arbitrary function minimizes the L-inf error among all polynomials of the same degree, because Chebyshev polynomials have minimum L-inf norm among all monic polynomials of the same degree. Let's try this: Tell me how you would construct a minimax approximation according to your definition?

The truncated Chebyshev series does not generally minimize the L-inf error among all polynomials of the same degree (unless the function being approximated is a polynomial of degree n+1). Constructing a minimax approximation is typically done via the Remez exchange algorithm, which is iterative, fussy, and prone to convergence failures. However, none of those matter when designing an offline approximation as is the c…

Thank you, that was enlightening.

Re: Implementing Advanced Math Functions

#23

Post by ajross » Wed, Apr 25, 2012, 7:28 PM UTC

jensnockert wrote: post #5

adrianN wrote: post #2

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

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/

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.

Re: Implementing Advanced Math Functions

#24

Post by balsam » Wed, Apr 25, 2012, 9:10 PM UTC

stephencanon wrote: post #21

jvranish wrote: post #17

Earlier quoted context omitted.

Do you know of good resource for how to compute minimax polynomials that approximate transcendental functions? My Google foo appears to be failing me.

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.

Mathematica, if you have access to it, can give you the minimax polynomial in one line.