The standard process for computation of sqrt(x) has been the use of Newton's method. This means that an approximation to SQRT(x) is found by whatever means and the iteration result = (result_old^2+x)/(2*result_old) Repeated until result = result_old. A check isn't necessary because a maximum iteration count can be obtained from the worst-case accuracy of the initial approximation. Thus there isn't any reason that the librarys should vary in their results.
A bit or so may result from the non-algebraically-reduced form of the iteration, result = result_old + (x - result_old^2)/(2*result_old) which some of little faith use because the second term can be tested for zero - or near zero. Don't go there. In any case this non-reduced form is less accurate. All this works for complex numbers, too. Watch out for the sign convention for complex numbers, though. James K Beard -----Original Message----- From: K. Frank [mailto:kfrank2...@gmail.com] Sent: Friday, April 29, 2011 3:56 PM To: mingw64; MinGW Users List Subject: Re: [Mingw-w64-public] Math library discrepancies that surprised me. 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 quality-of-implementation 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 64-bit tested) no output > Mac (snow leopard) no output > gcc (cygwin) no output > gcc-3 -mno-cygwin ) > i686-w64-mingw32-gcc ) all different from each other > x86_64-w64-mingw32-gcc ) 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 run-time 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, IEEE-754 compliance requires sqrt to be calculated "correctly", i.e., to return the floating-point number closest to the exact mathematical value. (The c++ standard, I believe, does not, however, require IEEE-754 compliance.) I think, however, the IEEE-754 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 linux-based gcc, and see whether the problem is in sqrt (for which IEEE-754 specifies the result) or in pow (for which I think IEEE-754 permits some leeway). In my opinion: mingw / msvcrt should be IEEE-754 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 IEEE-754). 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@panamint temp > $ gcc expt.c -o expt > acn1@panamint temp > $ ./expt > acn1@panamint temp > $ gcc-3 -mno-cygwin expt.c -o expt > acn1@panamint temp > $ ./expt > 1.0241210164238141 1.011988644414459 1.0119886444144588 > acn1@panamint temp > $ i686-w64-mingw32-gcc expt.c -o expt > acn1@panamint temp > $ ./expt > 1.0004603659307831 1.000230156479389 1.0002301564793892 > acn1@panamint temp > $ x86_64-w64-mingw32-gcc expt.c -o expt > acn1@panamint temp > $ ./expt > 1.0000179101601867 1.0000089550399969 1.0000089550399971 > acn1@panamint temp > ============================================================== ---------------------------------------------------------------------------- -- 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-w64-public mailing list Mingw-w64-public@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/mingw-w64-public ------------------------------------------------------------------------------ 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-w64-public mailing list Mingw-w64-public@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
