Date: Fri, 12 May 2000 16:52:49 -0500 From: Eric Rudd Subject: Re: Math functions To: djgpp-workers AT delorie DOT com Message-id: <391C7D31.838DE634@cyberoptics.com> Organization: CyberOptics MIME-version: 1.0 X-Mailer: Mozilla 4.72 [en] (Win95; U) Content-type: text/plain; charset=us-ascii Content-transfer-encoding: 7bit X-Accept-Language: en References: <200005121903 DOT PAA29374 AT delorie DOT com> Reply-To: djgpp-workers AT delorie DOT com 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