From: K. F. <kfr...@gm...> - 2011-05-03 13:32:30
|
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 |