Mail Archives: djgpp-workers/2000/05/12/17:13:31
Dieter Buerssner wrote:
> I see, but I am not sure, whether I agree. Unfortunately, this makes the task
> much more tedious.
Yes, I know. I spent perhaps 1/3 of the time actually developing algorithms
that would calculate the functions accurately, and 2/3 of the time in making
the error handling bullet-proof and fast.
> I think the range reduction for expl() is not that difficult, because of the
> limited range of the argument. When expl() and expm1l() are implemented,
> sinhl(), coshl() and tanhl() are not that difficult either, when an maximum
> errror of say 2ULPs = 2^-62 is acceptable.
I haven't looked at Moshier's code lately, but the inherent problem with range
reduction for functions like expl() is that the range reduction needs to be
done to higher precision than that of the final result, because of the
magnfication of the relative error of the argument. It turns out that for
expl(), with an output significand of 64 bits, and an exponent range of +/-
16384 (or 15 bits), you need about 15 bits of extra precision in the range
reduction, or about 64+15=79. The 64-bit precision of the coprocessor was
adequate for range reduction in a double routine, but not long double.
> Sin(), cos() and tan() are more difficult. My adapted range reduction code
> (from Stephen Moshier) seems to produce correct results up to 2^64. (I tested
> it against 100 decimal digit code for linear and exponentially scaled random
> arguments.) When really a larger range needs to be supported, it gets much
> more difficult.
It strikes me that it's difficult for any argument where fabs(x) >> 1. K. B.
Williams pointed me to some technical articles showing how to compute a
correctly-rounded double result for any double argument:
http://www.validgh.com/arg.ps
but it seemed like an unreasonable amount of extra execution time for
"workhorse" code that might well be used where speed is critical. For long
double routines, you also have to hard-code 1/pi to approximately 16000 bits,
so even the code size becomes an issue. You might find the article
interesting, though; they had some clever tricks to reduce the order of
complexity of the multiple-precision computation.
> The most difficult is IHMO powl(), and I don't have a good solution
> for that.
powl() is definitely one of the most difficult functions to compute. If one
carries internal results to 80 bits, normal arguments can be handled without
undue difficulty, but there are about 50 exception conditions that must be
handled -- powl(inf, 0.), powl(0., NaN), powl(-1., inf), and so forth. (Let's
see... is infinity an integer? ;-) My own code doesn't even get them
completely correct, in that the signs of infinities are not correctly handled,
but the problem exhausted me.
> From my testing, the floating point functions built into the FPU are normally
> better than 2^-63 (when called with an argument correctly reduced to the
> supported range). When this is an acceptable error (which I think it is), the
> functions won't be slow. When this is not an acceptable error, I can't do it
> :-(
Assuming that you find some good way to do range reduction, the coprocessor's
internal accuracy seems adequate.
> I do have gamma[l], lgamma[l], erf[l], erfc[l] and the Bessel functions
> (again adapted from Steven Moshier's code).
Those should help. I suspect that Stephen Moshier is actually half a dozen
people -- I don't see how a single person could have written all the code he
has.
-Eric Rudd
rudd AT cyberoptics DOT com
- Raw text -