
Implementing Advanced Math Functions - silentbicycle
http://spin.atomicobject.com/2012/04/24/implementing-advanced-math-functions/
======
stephencanon
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
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
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
Do you know of good resource for how to compute minimax polynomials that
approximate transcendental functions? My Google foo appears to be failing me.

~~~
stephencanon
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.

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

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

~~~
jensnockert
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/>

~~~
ajross
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
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
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/> [1][http://www.axiom-developer.org/axiom-
website/documentation.h...](http://www.axiom-developer.org/axiom-
website/documentation.html)

------
jeremysalwen
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>

------
evincarofautumn
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);
        }

