From: K. Frank <kfrank29.c@gm...>  20110430 00:51:07

Hi Greg! On Fri, Apr 29, 2011 at 7:38 PM, Greg Chicares <gchicares@...> wrote: > On 20110429 19:55Z, K. Frank wrote: >> >> By the way, could someone confirm that mingw does use msvcrt for >> sqrt and pow? > > http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src > > Looks to me like sqrt(double) calls into msvcrt. But pow() is a > different story. Danny implemented it in libmingwex because of > some qualityofimplementation problem with msvcrt; then people > complained that the libmingwex version was slower, and I can't > remember where it wound up, but it's all in the archives. Thanks for that  good to know. From the source link you gave: http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/pow.c?annotate=1.1.10.1&cvsroot=src the code for pow (the regular double pow (double, double)) has these comments about accuracy. 28: * ACCURACY: 29: * 30: * Relative error: 31: * arithmetic domain # trials peak rms 32: * IEEE 26,26 30000 4.2e16 7.7e17 33: * DEC 26,26 60000 4.8e17 9.1e18 34: * 1/26 < x < 26, with log(x) uniformly distributed. 35: * 26 < y < 26, y uniformly distributed. 36: * IEEE 0,8700 30000 1.5e14 2.1e15 37: * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. This commentary suggest that pow is not "exact." This would reasonably explain Arthur seeing pow not match sqrt. But it still leaves the puzzle of why the different mingw versions don't behave the same. One minor puzzle: I don't see anything in the pow code that would let the compiler / library take advantage of hardware support for pow, which I think has been available in the x86 architecture for quite some time now. > >> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >> to by IEEE754). > > If they should give exactly the same result always, then the > first thing that needs to be fixed is any standard that allows > them to differ. Well, it's a question of engineering compromise whether they should always agree. In an ideal world, it would be nice if pow were "exact" (in the sense of closet floatingpoint number). Then if pow were "exact," pow and sqrt should agree. But there is expense in getting pow to be "exact," so it's arguable that the language / IEEE standards should give implementers some flexibility. (And, as you point out below, if pow isn't exact, then forcing pow and sqrt to agree could actually cause other problems.) > The language standards are so much looser than > the numerical standard, even now, that I'd rather see C and C++ > catch up to IEEE7541985 before asking them to go beyond. > >> Browsing some gcc lists, I did see some comments that suggest that >> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x). This would >> make sqrt and pow consistent (whether or not technically correct). > > Not so many years ago (perhaps when msvcrt was written), it was > thought that correctlyrounded transcendental functions weren't > feasible in practice, so library authors had a looser attitude > toward exactness. If you have an ideal implementation of sqrt() > but a pow() that's only approximate, should you specialcase > the power of exactly 0.5? If you do, then pow() and sqrt() will > be consistent, but then the powers > x^(0.5*(1epsilon)) > x^(0.5 ) > x^(0.5*(1+epsilon)) > might not be monotone. Some would call that a bad tradeoff. I actually agree with your concern. The transformation of pow (x, 0.5) to sqrt (x) that it sounded like gcc / linux was doing sounds like a bad idea if it's being done this simply and naively. Any speculation on why Arthur gets different results with different mingw's? Thanks for your perspective. K. Frank 