From: Kai Tietz <ktietz70@go...>  20110430 08:42:15

2011/4/30 Greg Chicares <gchicares@...>: > 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. > >> 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. 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. Well, I spent some efforts into log, exp, and pow implementation. Also some other basemath functions I improved in mingww64's math library in two terms. a) Make those functions behave as ISOC99 specification tells in cornercases. and b) Improve accuracy in terms of gcc's internal used gmp. As we noticed here result differences. The following example illustrates the issues pretty well: #include <stdio.h> #include <math.h> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; int main() { double r[3]; long double rl[3]; int i; for (i = 0; i < 3; i++) __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], (r[i] = pow (abc[i][0], abc[i][1]))); __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, 3.0), pow (4.0, 0.5)); r[0] = pow (2.2, 3.1); r[1] = pow (2.0, 3.0); r[2] = pow (4.0, 0.5); __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); for (i = 0; i < 3; i++) __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], (rl[i] = powl (ab[i][0], ab[i][1]))); __mingw_printf ("%.20LF %.20LF %.20LF\n", powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl (4.0L, 0.5L)); rl[0] = powl (2.2L, 3.1L); rl[1] = powl (2.0L, 3.0L); rl[2] = powl (4.0L, 0.5L); __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); return 1; } The interesting issue is here that gcc uses for constant mathcalculations gmplibrary to get result. For more complex variants, it uses crt's math. By this it is important for a gcc based runtime to be in the IEEE 754 floating point compatible to gmp's results. Btw new pow, log, and exp variant are slightly faster then variant in msvcrt, but this is more or less just a sideeffect. I admit that this might be not that obvious to users, why results here are different, but for a gcc based toolchain we need to play nice with gcc's internal assumptions. Regards, Kai 