From: K. Frank <kfrank29.c@gm...>  20110429 19:56:00

Hello Arthur! (I've take the liberty of copying this to the mingw list, as well.) On Fri, Apr 29, 2011 at 9:54 AM, Arthur Norman wrote: > The following short program compares the value of pow(x, 0.5) to that > returned by sqrt(x) for a range of values of x, and prints out any first > case where the two do not agree. Well, this is unfortunate. Ideally they should, and as a qualityofimplementation issue, I think we should be able to expect them to agree. > With the mingw family of compilers I had > expected the Microsoft C library to be in use. > > What I appear to see is > linux (ubuntu 11.04 64bit tested) no output > Mac (snow leopard) no output > gcc (cygwin) no output > gcc3 mnocygwin ) > i686w64mingw32gcc ) all different from each other > x86_64w64mingw32gcc ) but all find discrepancies. Very interesting. > It is a totally fair cop if the Microsoft C library gives different > results in the last bit for floating point elementary functions to gnu > libraries, Well, from the mingw perspective, if the microsoft runtime library is flawed in this way, then it's a legitimate cop for mingw not to fix it. By the way, could someone confirm that mingw does use msvcrt for sqrt and pow? > but is there an obvious reason why the three mingw varients > all give different answers. If I understand your point correctly, if the problem is with msvcrt, then the three mingw variants should all give the same flawed results. > This showed up in regression testing something > much larger across platforms. > > Arthur > ... As I understand it, IEEE754 compliance requires sqrt to be calculated "correctly", i.e., to return the floatingpoint number closest to the exact mathematical value. (The c++ standard, I believe, does not, however, require IEEE754 compliance.) I think, however, the IEEE754 does not require all mathematical functions, and, in particular, pow, to be "correct," in the above sense sense. It would be interesting to print out in binary (or hex) the results for a few of the values that lead to discrepancies with both mingw and a linuxbased gcc, and see whether the problem is in sqrt (for which IEEE754 specifies the result) or in pow (for which I think IEEE754 permits some leeway). In my opinion: mingw / msvcrt should be IEEE754 compliant. Of course, if msvcrt isn't, then it's not mingw's fault. sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required to by IEEE754). 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). Maybe gcc / linux pow and msvcrt pow are both not quite exact, but this is hidden in your test, because gcc / linux replaces pow (x, 0.5) with sqrt (x), and msvcrt doesn't. (Just speculation.) It is odd, however, that the various mingw's don't agree, and this speaks to a possible bug in mingw. Thanks. K. Frank > ... > ==================== > #include <stdio.h> > #include <math.h> > > int main(int argc, char *argv[]) > { > double a, b, c; > int i; > for (i=0, a=1.0; i<10000000; i++) > { a *= 1.00000001; > b = sqrt(a); c = pow(a, 0.5); > if (b == c) continue; > printf("%.17g %.17g %.17g\n", a, b, c); > return 0; > } > return 0; > } > ==================== > > acn1@... temp > $ gcc expt.c o expt > acn1@... temp > $ ./expt > acn1@... temp > $ gcc3 mnocygwin expt.c o expt > acn1@... temp > $ ./expt > 1.0241210164238141 1.011988644414459 1.0119886444144588 > acn1@... temp > $ i686w64mingw32gcc expt.c o expt > acn1@... temp > $ ./expt > 1.0004603659307831 1.000230156479389 1.0002301564793892 > acn1@... temp > $ x86_64w64mingw32gcc expt.c o expt > acn1@... temp > $ ./expt > 1.0000179101601867 1.0000089550399969 1.0000089550399971 > acn1@... temp > ============================================================== 