From: James K B. <bea...@ve...> - 2011-05-03 16:37:55
|
There have been special cases in scientific libraries for quite a long time. The special case for atan(1.0) is used by many programmers to extract the system stored value for pi/4, knowing that that the math library (or co-processor) wouldnt compute it instead of simply retrieving the known result. The simple fact that a function like pow(*,*) will have an accuracy requirement that allows a bit or two of numerical error while the nature of the special case sqrt(*) allows full mantissa accuracy to within a possible toggle of the last bit, is not the same thing. The arc tangent is computed by libraries by expansion about the argument points 0 and 1, and identities like arctan(|x|) = pi/2 - arctan(1/|x|) when |x|>1 and such. But this doesn't apply to pow(*,*) because, other than extraction of the exponents of the two arguments to limit the range of the logarithm and exponential in the identity a^b = exp(b*ln(a)) there is no breaking up of things with a natural boundary at b=0.5, although special cases for b=0 and 1 are often taken, and sometimes checks are made to determine whether or not b is an integer. The usual process is breaking up a into exponent and mantissa, a = (2^ea)*ma then computation of the logarithm of the result, ln(a^b) = lr = b*ln(a) = b*(eb*ln(2)+ln(ma)) = ilr + mlr, |mlr|<1 and finally exponentiation of the log result, a^b = (e^ilr)*(e^mlr). There is no obvious advantage for checking for any particular value of b other than 0 and 1 and possibly integers. Thus it seems to me, the obvious thing to do is to accept the specified accuracy of pow(*,*) and not do special cases other than 0, 1, and possibly integers. One might stretch things a bit by evaluating what we would call pow(a^2, b/2) but that would not likely survive an analysis of speed and accuracy as compared to pow(a,b). James K Beard -----Original Message----- From: K. Frank [mailto:kfr...@gm...] Sent: Tuesday, May 03, 2011 9:32 AM To: mingw64; MinGW Users List Subject: Re: [Mingw-users] [Mingw-w64-public] Math library discrepancies that surprised me. Hello Kai and Everybody! A couple of comments about whether to use special cases for pow (x, y) for things like y = 0.5 or y = integer... On Tue, May 3, 2011 at 6:58 AM, Kai Tietz wrote: > 2011/5/3 Peter Rockett <p.r...@sh...>: >> On 03/05/2011 10:11, Kai Tietz wrote: >> ... >> The reason is not introducing a discontinuity into pow() around 0.5. >> Forcing the two functions to agree may seem like a good idea but in >> reality, you are creating another much more insidious problem. You are >> embedding a potentially disastrous property into pow(). >> >>> It might be that I haven't seen the forbidden >>> thing on the spec. The sqrt special-case has at least one advantage >>> over the log2/exp2 variant. First it is faster, as just one FPU >>> command is used to calculate sqrt. Secondly, is provides better result >>> and compatible one to gcc's internal used soft-floating-point library. >>> As I don't see in spec that this special casing shouldn't be done, I >>> think we will keep it on mingw-w64's side. >> ... > I am aware of this what you are telling me, but there is in pow > function already special-casing for x^y for y with integer values in > range for y >= INT_MIN and y <= INT_MAX. As here more precise > calculation via powi is used. This makes that results of pow having > already discontinuity. > So pow functions doesn't have the requirement of having continuity. > ... I'm not aware that either the C standard or the IEEE-754 standard specifies one behavior or the other. (I think that they don't, but I could be wrong.) I believe that is has been commonplace over the years to use these special cases for pow. (This is just the sense I have from personal recollection and oral tradition.) I think speed optimization has been the reason, rather than improved accuracy or agreement with things like sqrt, but perhaps it's been a mixture of those. Of course, just because there is precedent, doesn't mean it's the right way to go. If I had to choose between pow (x, 0.5) agreeing with sqrt (x) and pow (x, y) being monotone in y, I think I would choose the latter (i.e., choose not to use the special cases). But there are good arguments on both sides. One solution that would be attractive to me would be to implement pow "exactly" (i.e., bankers' rounding to the nearest representable value), and then use all the special cases you want as a speed optimization. Now the special cases won't change the value of pow (assuming sqrt is properly implemented), so the special cases won't introduce any non-monotonicity or "discontinuities" into pow. What I don't know is how hard or expensive this is. I'm pretty sure it's not too hard to calculate exp "exactly." But I think that calculating log is harder. Perhaps James Beard knows whether it's practical to calculate pow "exactly," and maybe even how to do it. (Another approach, which I very much don't recommend, is to use pow to calculate sqrt. You get reduced speed and reduced accuracy for sqrt, but at least you do get sqrt and pow to agree without introducing non-monotonicity into pow.) Lastly, I have to believe that this is a well-understood and "solved" problem in the numerical-analysis community (of which I am not a part). Perhaps it would make sense to find out how the experts who have doing this for fifty-plus years deal with this issue. > Regards, > Kai And thanks to the mingw and mingw64 teams who have invested the time and effort to actually write the code we're so breezily discussing. My comments are in no way meant to be a criticism of those efforts. I'm just offering up one man's opinion from the user community in the hope that it might be helpful. Bets regards. K. Frank ---------------------------------------------------------------------------- -- WhatsUp Gold - Download Free Network Management Software The most intuitive, comprehensive, and cost-effective network management toolset available today. Delivers lowest initial acquisition cost and overall TCO of any competing solution. http://p.sf.net/sfu/whatsupgold-sd _______________________________________________ MinGW-users mailing list Min...@li... This list observes the Etiquette found at http://www.mingw.org/Mailing_Lists. We ask that you be polite and do the same. Disregard for the list etiquette may cause your account to be moderated. _______________________________________________ You may change your MinGW Account Options or unsubscribe at: https://lists.sourceforge.net/lists/listinfo/mingw-users Also: mailto:min...@li...?subject=unsubscribe |