From: K. F. <kfr...@gm...> - 2011-04-29 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 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 > ============================================================== |
From: James K B. <bea...@ve...> - 2011-04-30 13:36:33
|
My comments don't seem to be making it into the thread this week, possibly because of a switch between my reply-to email address and another. So, if it gets lost, perhaps nightstrike can intervene, or echo the information. Perhaps not. The sqrt() function can be implemented with excellent speed to within a bit of full mantissa accuracy using iterations of Newton's method, as is classically done in libraries. Not so with ln() or exp(), which are classically done by extracting the exponent and using the reduced-range mantissa in a modified Chebychev approximation. Arithmetic co-processors do sqrt(), ln(), and exp() internally with 80-bit floating point, then store in 32, 64, or 80 bit floating point format. So, the co-processors can get full mantissa accuracy with 32-bit and 64-bit sqrt() and pow(x,0.5). So, the problem is non-existent for Intel/PPC and most other common workstation processors. Yet we are looking at software transcendental function libraries. James K Beard -----Original Message----- From: Greg Chicares [mailto:gch...@sb...] Sent: Friday, April 29, 2011 7:38 PM To: MinGW Users List Cc: mingw64 Subject: Re: [Mingw-w64-public] [Mingw-users] Math library discrepancies that surprised me. On 2011-04-29 19:55Z, K. Frank wrote: > > By the way, could someone confirm that mingw does use msvcrt for > sqrt and pow? http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). 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 IEEE-754-1985 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 correctly-rounded 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 special-case the power of exactly 0.5? If you do, then pow() and sqrt() will be consistent, but then the powers x^(0.5*(1-epsilon)) x^(0.5 ) x^(0.5*(1+epsilon)) might not be monotone. Some would call that a bad tradeoff. ---------------------------------------------------------------------------- -- 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 Min...@li... https://lists.sourceforge.net/lists/listinfo/mingw-w64-public |
From: Kai T. <kti...@go...> - 2011-04-30 15:34:00
|
2011/4/30 K. Frank <kfr...@gm...>: > Hello Kai! > > On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <kti...@go...> wrote: >> 2011/4/30 Kai Tietz <kti...@go...>: >>> 2011/4/30 Greg Chicares <gch...@sb...>: >>>> On 2011-04-29 19:55Z, K. Frank wrote: >>>>> >>>>> By the way, could someone confirm that mingw does use msvcrt for >>>>> sqrt and pow? >>>> >>>> http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). >>>> ... >>>> Not so many years ago (perhaps when msvcrt was written), it was >>>> thought that correctly-rounded 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 special-case >>>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>>> be consistent, but then the powers >>>> x^(0.5*(1-epsilon)) >>>> 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 base-math functions I improved in mingw-w64's math library >>> in two terms. a) Make those functions behave as ISO-C99 specification >>> tells in corner-cases. 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; >>> } > > Do I understand your point correctly that you are comparing gcc's > compile-time and run-time values for pow? > > For the record, when I run this on a mingw64 build: > > gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] > > I get: > > C:\>gcc -o pow_test pow_test.c > > C:\>pow_test > 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 > 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 > 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 > 11.52153412678571697825 8.00000000000000000000 2.00000000000000000000 > 0.00000000000000177636 0.00000000000000000000 0.00000000000000000000 > 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 > 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 > 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 > 11.52153412678571414718 8.00000000000000000000 2.00000000000000000000 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > > As I understand it, this means I am seeing a small difference between > compile and run-time evaluation only for powl (2.2, 3.1), i.e., the > long-double case. > >>> The interesting issue is here that gcc uses for constant >>> math-calculations gmp-library 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. > > Yes, that's a good point. I hadn't considered the issue of compile-time > evaluation, and, ideally, it should agree with run-time evaluation, adding > yet another constraint to our wish list. > > Furthermore, if we take the idea of a cross-compiler seriously, we > would also like the compile-time evaluations of a cross-compiler and > native compiler to agree. > > Hmm... > > So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like > compile-time and run-time evaluations to agree. We would like > cross-compilers and native compilers to agree. > > Taken at face value, this would seem to preclude using the platform > library and platform hardware to compute things like sqrt and pow > (unless we know they comply with some agreed-upon standard, > such as IEEE-754). > > Just to make the reasoning explicit, if I perform my run-time evaluation > on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler > running on platform B won't have platform A's FPU available to perform > consistent compile-time evaluations. > > So you either give up some of the above-wished-for consistency, you > forgo the potential performance benefits of using the platform's library > and / or hardware, or you get everybody to use a common standard. > >>> ... >>> 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 >> >> On closer look for long double, we can have rouding issues for long >> double precission. I see that powl (0.01L, 0.5L) we have a rounding >> difference to sqrtl (0.01L). Interestingly not for float or double >> types. So we we might need to special-case in pow the case for y == >> 0.5 to solve this. > > Be, however, cognizant of Greg's concern. Forcing pow to agree with > sqrt just for y = 0.5 has the potential to make pow, in a sense, internally > inconsistent, It's a tricky trade-off. (Ideally, both pow and sqrt would > be "exact" and therefore agree with one another automatically.) > > By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L). > It disagrees slightly from sqrtl (2.2L) and from the compile-time value > for powl (2.2L, 0.5L). > > I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl > with powl). For the record: > > C:\>gcc --version > gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] > > C:\>gcc -o sqrt_pow sqrt_pow.c > > C:\>sqrt_pow > 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000 > 1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000 > 1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000 > 1.04880884817015163080 1.04880884817015163080 > 1.48323969741913264109 1.48323969741913264109 > 1.81659021245849494619 1.81659021245849494619 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 > 1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011 > 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 > 1.04880884817015154699 1.04880884817015154699 > 1.48323969741913258980 1.48323969741913258980 > 1.81659021245849499921 1.81659021245849499921 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000 > > where sqrt_pow.c is: > > #include <stdio.h> > #include <math.h> > > double abc[3] = { 1.1, 2.2, 3.3 }; > long double abcl[3] = { 1.1L, 2.2L, 3.3L }; > int main() > { > double r1[3], r2[3], r3[3]; > long double rl1[3], rl2[3], rl3[3]; > int i; > for (i = 0; i < 3; i++) > __mingw_printf ("%.20F, %.20F, %.20F\n", > (r1[i] = sqrt (abc[i])), > (r2[i] = pow (abc[i], 0.5)), > (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i]))); > __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5)); > __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5)); > __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5)); > r1[0] -= sqrt (1.1); > r1[1] -= sqrt (2.2); > r1[2] -= sqrt (3.3); > __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]); > r2[0] -= pow (1.1, 0.5); > r2[1] -= pow (2.2, 0.5); > r2[2] -= pow (3.3, 0.5); > __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]); > > for (i = 0; i < 3; i++) > __mingw_printf ("%.20LF, %.20LF, %.20LF\n", > (rl1[i] = sqrtl (abcl[i])), > (rl2[i] = powl (abcl[i], 0.5L)), > (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i]))); > __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L)); > __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L)); > __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L)); > rl1[0] -= sqrtl (1.1L); > rl1[1] -= sqrtl (2.2L); > rl1[2] -= sqrtl (3.3L); > __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]); > rl2[0] -= powl (1.1L, 0.5L); > rl2[1] -= powl (2.2L, 0.5L); > rl2[2] -= powl (3.3L, 0.5L); > __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]); > > return 1; > } > >> >> Kai > > Thanks, Kai, for your insight and explanations. > > > K. Frank I add a special case to mingw-w64's pow(f|l) function about y=0.5 and get now for your adjusted testcase the following results, which are looking much better IMHO. $ ./tst.exe 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000004142 1.48323969741913264109, 1.48323969741913264109, -0.00000000000000000857 1.81659021245849494619, 1.81659021245849494619, -0.00000000000000000412 1.04880884817015163080 1.04880884817015163080 1.48323969741913264109 1.48323969741913264109 1.81659021245849494619 1.81659021245849494619 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 1.04880884817015154699 1.04880884817015154699 1.48323969741913258980 1.48323969741913258980 1.81659021245849499921 1.81659021245849499921 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 I used for this a 32-bit variant and 64-bit variant, which made no differences in output. The issue here is that 1.0 branch still uses the inaccurate variant from Danny, but nowaday's trunk uses for pow, sqrt, log, and exp functions new implementations, which are treating rounding and short-path calculation in an ISO-C variant (as gcc expects internally such a math). Also complex-math part is completely rewritten on mingw-w64's trunk version. IMHO it might be a good idea to replace step by step the classical runtime functions, so that they are not using here msvcrt in first place. Of course the import-symbols of msvcrt (for compatibility with VC generated objects need to remain in .def exports, but the actual functions should come directly from mingwex. |
From: Kai T. <kti...@go...> - 2011-04-30 15:48:06
|
2011/4/30 Kai Tietz <kti...@go...>: > 2011/4/30 K. Frank <kfr...@gm...>: >> Hello Kai! >> >> On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <kti...@go...> wrote: >>> 2011/4/30 Kai Tietz <kti...@go...>: >>>> 2011/4/30 Greg Chicares <gch...@sb...>: >>>>> On 2011-04-29 19:55Z, K. Frank wrote: >>>>>> >>>>>> By the way, could someone confirm that mingw does use msvcrt for >>>>>> sqrt and pow? >>>>> >>>>> http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). >>>>> ... >>>>> Not so many years ago (perhaps when msvcrt was written), it was >>>>> thought that correctly-rounded 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 special-case >>>>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>>>> be consistent, but then the powers >>>>> x^(0.5*(1-epsilon)) >>>>> 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 base-math functions I improved in mingw-w64's math library >>>> in two terms. a) Make those functions behave as ISO-C99 specification >>>> tells in corner-cases. 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; >>>> } >> >> Do I understand your point correctly that you are comparing gcc's >> compile-time and run-time values for pow? >> >> For the record, when I run this on a mingw64 build: >> >> gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] >> >> I get: >> >> C:\>gcc -o pow_test pow_test.c >> >> C:\>pow_test >> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 >> 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 >> 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 >> 11.52153412678571697825 8.00000000000000000000 2.00000000000000000000 >> 0.00000000000000177636 0.00000000000000000000 0.00000000000000000000 >> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 >> 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 >> 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 >> 11.52153412678571414718 8.00000000000000000000 2.00000000000000000000 >> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 >> >> As I understand it, this means I am seeing a small difference between >> compile and run-time evaluation only for powl (2.2, 3.1), i.e., the >> long-double case. >> >>>> The interesting issue is here that gcc uses for constant >>>> math-calculations gmp-library 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. >> >> Yes, that's a good point. I hadn't considered the issue of compile-time >> evaluation, and, ideally, it should agree with run-time evaluation, adding >> yet another constraint to our wish list. >> >> Furthermore, if we take the idea of a cross-compiler seriously, we >> would also like the compile-time evaluations of a cross-compiler and >> native compiler to agree. >> >> Hmm... >> >> So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like >> compile-time and run-time evaluations to agree. We would like >> cross-compilers and native compilers to agree. >> >> Taken at face value, this would seem to preclude using the platform >> library and platform hardware to compute things like sqrt and pow >> (unless we know they comply with some agreed-upon standard, >> such as IEEE-754). >> >> Just to make the reasoning explicit, if I perform my run-time evaluation >> on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler >> running on platform B won't have platform A's FPU available to perform >> consistent compile-time evaluations. >> >> So you either give up some of the above-wished-for consistency, you >> forgo the potential performance benefits of using the platform's library >> and / or hardware, or you get everybody to use a common standard. >> >>>> ... >>>> 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 >>> >>> On closer look for long double, we can have rouding issues for long >>> double precission. I see that powl (0.01L, 0.5L) we have a rounding >>> difference to sqrtl (0.01L). Interestingly not for float or double >>> types. So we we might need to special-case in pow the case for y == >>> 0.5 to solve this. >> >> Be, however, cognizant of Greg's concern. Forcing pow to agree with >> sqrt just for y = 0.5 has the potential to make pow, in a sense, internally >> inconsistent, It's a tricky trade-off. (Ideally, both pow and sqrt would >> be "exact" and therefore agree with one another automatically.) >> >> By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L). >> It disagrees slightly from sqrtl (2.2L) and from the compile-time value >> for powl (2.2L, 0.5L). >> >> I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl >> with powl). For the record: >> >> C:\>gcc --version >> gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] >> >> C:\>gcc -o sqrt_pow sqrt_pow.c >> >> C:\>sqrt_pow >> 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000 >> 1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000 >> 1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000 >> 1.04880884817015163080 1.04880884817015163080 >> 1.48323969741913264109 1.48323969741913264109 >> 1.81659021245849494619 1.81659021245849494619 >> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 >> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 >> 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 >> 1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011 >> 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 >> 1.04880884817015154699 1.04880884817015154699 >> 1.48323969741913258980 1.48323969741913258980 >> 1.81659021245849499921 1.81659021245849499921 >> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 >> 0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000 >> >> where sqrt_pow.c is: >> >> #include <stdio.h> >> #include <math.h> >> >> double abc[3] = { 1.1, 2.2, 3.3 }; >> long double abcl[3] = { 1.1L, 2.2L, 3.3L }; >> int main() >> { >> double r1[3], r2[3], r3[3]; >> long double rl1[3], rl2[3], rl3[3]; >> int i; >> for (i = 0; i < 3; i++) >> __mingw_printf ("%.20F, %.20F, %.20F\n", >> (r1[i] = sqrt (abc[i])), >> (r2[i] = pow (abc[i], 0.5)), >> (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i]))); >> __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5)); >> __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5)); >> __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5)); >> r1[0] -= sqrt (1.1); >> r1[1] -= sqrt (2.2); >> r1[2] -= sqrt (3.3); >> __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]); >> r2[0] -= pow (1.1, 0.5); >> r2[1] -= pow (2.2, 0.5); >> r2[2] -= pow (3.3, 0.5); >> __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]); >> >> for (i = 0; i < 3; i++) >> __mingw_printf ("%.20LF, %.20LF, %.20LF\n", >> (rl1[i] = sqrtl (abcl[i])), >> (rl2[i] = powl (abcl[i], 0.5L)), >> (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i]))); >> __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L)); >> __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L)); >> __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L)); >> rl1[0] -= sqrtl (1.1L); >> rl1[1] -= sqrtl (2.2L); >> rl1[2] -= sqrtl (3.3L); >> __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]); >> rl2[0] -= powl (1.1L, 0.5L); >> rl2[1] -= powl (2.2L, 0.5L); >> rl2[2] -= powl (3.3L, 0.5L); >> __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]); >> >> return 1; >> } >> >>> >>> Kai >> >> Thanks, Kai, for your insight and explanations. >> >> >> K. Frank > > I add a special case to mingw-w64's pow(f|l) function about y=0.5 and > get now for your adjusted testcase the following results, which are > looking much better IMHO. > > $ ./tst.exe > 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000004142 > 1.48323969741913264109, 1.48323969741913264109, -0.00000000000000000857 > 1.81659021245849494619, 1.81659021245849494619, -0.00000000000000000412 > 1.04880884817015163080 1.04880884817015163080 > 1.48323969741913264109 1.48323969741913264109 > 1.81659021245849494619 1.81659021245849494619 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 > 1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000 > 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 > 1.04880884817015154699 1.04880884817015154699 > 1.48323969741913258980 1.48323969741913258980 > 1.81659021245849499921 1.81659021245849499921 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 > > I used for this a 32-bit variant and 64-bit variant, which made no > differences in output. > The issue here is that 1.0 branch still uses the inaccurate variant > from Danny, but nowaday's trunk uses for pow, sqrt, log, and exp > functions new implementations, which are treating rounding and > short-path calculation in an ISO-C variant (as gcc expects internally > such a math). Also complex-math part is completely rewritten on > mingw-w64's trunk version. > > IMHO it might be a good idea to replace step by step the classical > runtime functions, so that they are not using here msvcrt in first > place. Of course the import-symbols of msvcrt (for compatibility with > VC generated objects need to remain in .def exports, but the actual > functions should come directly from mingwex. > PS: The differences shown in results for 32-bit aren't caused by inaccuarcy of used math directly. Instead it is caused by the fact that 32-bit use FPU registers (which are 80-bit always). By this the inline difference without using store, have for double small differences. By first storing values and then printing those stored values you see also for 32-bit a delta of zero. Interesting this behavior is on 64-bit x86_64 different, as it uses always precise accuracy to store results in registers. By this you have really just a 64-bit result for double type, and not indirecty 80-bit (as shown in first case). I altered your testcase a bit in following lines ... for (i = 0; i < 3; i++) { r1[i] = sqrt (abc[i]); r2[i] = pow (abc[i], 0.5); r3[i] = pow (abc[i], 0.5) - sqrt (abc[i]); __mingw_printf ("%.20F, %.20F, %.20F\n", r1[i], r2[i], r3[i]); } ... and get now $ ./tst.exe 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000 1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000 1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000 1.04880884817015163080 1.04880884817015163080 1.48323969741913264109 1.48323969741913264109 1.81659021245849494619 1.81659021245849494619 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 1.04880884817015154699 1.04880884817015154699 1.48323969741913258980 1.48323969741913258980 1.81659021245849499921 1.81659021245849499921 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 as result. it is an interesting fact that this happens. Regards, Kai |
From: JonY <jo...@us...> - 2011-05-01 14:38:07
Attachments:
signature.asc
0xED74C077.asc
|
On 5/1/2011 22:03, James K Beard wrote: > Now, about that math library problem -- it seems to me that asking full > mantissa accuracy and agreement for 32-bit and 64-bit floating point is > reasonable. Differences between FPUs makes full mantissa accuracy and > agreement at 80 bits a bit of a question, and I wouldn't demand it so that > the FPU can be used - unless you are using your own algorithms for quad > precision or other situation where the FPU doesn't provide the result. > > AFAIK, all currently produced FPUs are IEEE compliant. Are there some out > there that are still bucking the inevitable? -- Putting mingw lists back -- I heard older AMDs (late 90s-200Xs era) don't have 80-bit FPUs, but they do have 64-bit FPUs, so it wasn't strange that results do differ slightly between CPUs. Maybe that has changed for newer CPUs. Does the long double epsilon 1.08420217248550443401e-19L accuracy absolutely require an 80-bit FPU? IMHO, we don't need to produce results more accurate than that, but its a bonus if it doesn't impact performance too much. |
From: K. F. <kfr...@gm...> - 2011-05-01 15:28:53
|
Hello Jon! On Sun, May 1, 2011 at 10:17 AM, JonY <jo...@us...> wrote: > On 5/1/2011 22:03, James K Beard wrote: >> Now, about that math library problem... >> ... >> AFAIK, all currently produced FPUs are IEEE compliant. Are there some out >> there that are still bucking the inevitable? > > -- Putting mingw lists back -- > > I heard older AMDs (late 90s-200Xs era) don't have 80-bit FPUs, but they > do have 64-bit FPUs, so it wasn't strange that results do differ > slightly between CPUs. Maybe that has changed for newer CPUs. > > Does the long double epsilon 1.08420217248550443401e-19L accuracy > absolutely require an 80-bit FPU? If your willing to do your floating-point math in software, then no. But if you want to get 19-decimal-digit precision in your FPU, it pretty much needs to be 80 bits. (You can always trade off exponent bits for mantissa bits, but it's much better stick with IEEE choice.) The 80-bit floating-point format has a 64-bit mantissa (and, I believe, no implicit bit): log_base_10 (2^64) = 19.27, so 64 bits gets you a little over 19 decimal digits of precision. So to get 19 decimal digits, you need 64 bits plus exponent plus sign, and that gives you something pretty similar to the 80-bit format, and it would be hard to justify making a non-standard choice. > IMHO, we don't need to produce results more accurate than that, but its > a bonus if it doesn't impact performance too much. Historically the approach (within a compiler / platform family) has been to support whatever floating-point you support in software, and use floating-point hardware when available. If you want to do heavy floating-point computation, it's up to you to make sure you have the necessary floating-point hardware or eat the (significant) performance hit. Nowadays I would think that floating-point hardware is pretty much a given on desktops and servers. Maybe this isn't true for some embedded systems (that gcc targets) -- I just don't know. Best. K. Frank |
From: Greg C. <gch...@sb...> - 2011-05-01 15:11:35
|
On 2011-05-01 13:38Z, K. Frank wrote: [...] > I'm speaking from memory, but -- assuming that gcc long doubles are > x87 80-bit floating-point numbers -- yes, they do not use an implicit bit. Yes, MinGW gcc's 'long double' is intel's 'extended precision' 80-bit floating point. BTW, IIRC, sizeof(long double) is 12--for reasons of alignment, I guess--but only ten bytes are used. Indeed 80-bit extended precision has no implicit (hidden) bit. Here's an authoritative explanation: http://www.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF | 680x0 / ix87 Extendeds have an explicit Nth bit for historical reasons; | it allowed the Intel 8087 to suppress the normalization of subnormals | advantageously for certain scalar products in matrix computations, but | this and other features of the 8087 were later deemed too arcane to | include in IEEE 754, and have atrophied. |
From: K. F. <kfr...@gm...> - 2011-05-01 15:29:54
|
Hello Greg! Thanks for clearing this up. On Sun, May 1, 2011 at 11:11 AM, Greg Chicares <gch...@sb...> wrote: > On 2011-05-01 13:38Z, K. Frank wrote: > [...] >> I'm speaking from memory, but -- assuming that gcc long doubles are >> x87 80-bit floating-point numbers -- yes, they do not use an implicit bit. > > Yes, MinGW gcc's 'long double' is intel's 'extended precision' 80-bit > floating point. BTW, IIRC, sizeof(long double) is 12--for reasons of > alignment, I guess--but only ten bytes are used. > > Indeed 80-bit extended precision has no implicit (hidden) bit. Here's > an authoritative explanation: > > http://www.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF > | 680x0 / ix87 Extendeds have an explicit Nth bit for historical reasons; > | it allowed the Intel 8087 to suppress the normalization of subnormals > | advantageously for certain scalar products in matrix computations, but > | this and other features of the 8087 were later deemed too arcane to > | include in IEEE 754, and have atrophied. Thanks for this reference. Glad to have the definitive answer. Best. K. Frank |
From: Peter R. <p.r...@sh...> - 2011-05-02 19:42:09
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> <html> <head> <meta content="text/html; charset=ISO-8859-1" http-equiv="Content-Type"> </head> <body text="#000000" bgcolor="#ffffff"> Forgive the top-post but this whole disagreement seems concerned with <i>machine numbers</i>, that small set of floats that can be represented exactly: the trailing zeroes beyond the least-significant mantissa bit are indeed zeroes. Machine numbers have a important role in numerical algebra since they are zero-error numbers and hence probe the intrinsic accuracy of a routine rather than routine + error of input quantity convolved together in some way. But as far as I am aware, machine numbers are only ever used for testing numerical routines - since you cannot say that a general float is a machine number or just the nearest approximation to something nearby on the real number line, I am puzzled by the points annotated at the bottom of this post. Or do these refer only to testing?<br> <br> P.<br> <br> <br> On 02/05/2011 18:17, K. Frank wrote: <blockquote cite="mid:BANLkTi=aMb...@ma..." type="cite"> <pre wrap="">Hello Charles and Keith! On Mon, May 2, 2011 at 11:30 AM, Charles Wilson wrote: </pre> <blockquote type="cite"> <pre wrap="">On 5/2/2011 2:34 AM, Keith Marshall wrote: </pre> <blockquote type="cite"> <pre wrap="">On 01/05/11 14:36, K. Frank wrote: </pre> <blockquote type="cite"> <pre wrap="">[snip] </pre> </blockquote> <pre wrap="">My complaint is that, long after the last bit of precision has been interpreted, gdtoa (which is at the heart of printf()'s floating point output formatting) will continue to spew out extra decimal digits, based solely on the residual remainder from the preceding digit conversion, with an arbitrary number of extra zero bits appended. (Thus, gdtoa makes the unjustified and technically invalid assumption that the known bit precision may be arbitrarily extended to ANY LENGTH AT ALL, simply by appending zero bits in place of the less significant unknowns). </pre> <blockquote type="cite"> <pre wrap="">I don't agree with this. </pre> </blockquote> <pre wrap=""> Well, you are entitled to your own opinion; we may agree to disagree. </pre> <blockquote type="cite"> <pre wrap="">In most cases, it is not helpful to print out a long double to more than twenty decimal place, but sometimes it is. The point is that it is not the case that floating-point numbers represent all real numbers inexactly; rather, they represent only a subset of real numbers exactly. </pre> </blockquote> </blockquote> <pre wrap=""> But the problem is, if I send you a floating point number that represents the specific real number which I have in mind, exactly, YOU don't know that. </pre> </blockquote> <pre wrap=""> True. But this does not apply to all use cases. If I have generated a floating-point number, and I happen to know that it was generated in a way that it exactly represents a specific real number, I am fully allowed to make use of this information. (If you send me a floating-point number, and don't provide me the guarantee that it is being used to represent a specific real number exactly -- and this indeed is the most common, but not the only use case -- then your comments are correct.) </pre> <blockquote type="cite"> <pre wrap="">All you have is a particular floating point number that represents the range [value-ulps/2, value+ulps/w). You have no idea that I actually INTENDED to communicate EXACTLY "value" to you. </pre> </blockquote> <pre wrap=""> Yes, in your use case you did not intend to communicate an exact value to me. Therefore I should not imagine the floating-point number you sent me to be exact. And I shouldn't use gdtoa to print it out to fifty decimal places. </pre> <blockquote type="cite"> <pre wrap="">... </pre> <blockquote type="cite"> <blockquote type="cite"> <pre wrap="">If I happen to be representing a real number exactly with a long double, I might wish to print it out with lots (more than twenty) decimal digits. Such a use case is admittedly rare, but not illegitimate. </pre> </blockquote> </blockquote> <pre wrap=""> No, that's always illegitimate (i.e. misleading). </pre> </blockquote> <pre wrap=""> It's the "always" I disagree with. Your comments are correct for the use cases you describe, but you are incorrectly implying other use cases don't exist. </pre> <blockquote type="cite"> <pre wrap="">... </pre> <blockquote type="cite"> <pre wrap="">... </pre> <blockquote type="cite"> <pre wrap="">Let's say that I have a floating-point number with ten binary digits, so it gives about three decimal digits of precision (2^10 == 1024 ~= 10^3). I can use such a number to represent 1 + 2^-10 exactly. </pre> </blockquote> <pre wrap=""> Well, yes, you can if we allow you an implied eleventh bit, as most significant, normalised to 1; thus your mantissa bit pattern becomes: 10000000001B </pre> <blockquote type="cite"> <pre wrap="">I can print this number out exactly in decimal using ten digits after the decimal point: 1.0009765625. That's legitimate, and potentially a good thing. </pre> </blockquote> </blockquote> <pre wrap=""> But 10000000001B does NOT mean "1 + 2^-10". </pre> </blockquote> <pre wrap=""> In some specialized use cases, it means precisely this. In my program (assuming it's written correctly, etc.) the value means precisely what my program deems it to mean. </pre> <blockquote type="cite"> <pre wrap="">It means "with the limited precision I have, I can't represent the actual value of real number R more accurately with any other bit pattern than this one". </pre> </blockquote> <pre wrap=""> In the common use cases, this is a very good way to understand floating-point numbers. But, to reiterate my point, there are use cases where floating-point numbers are used to exactly represent a special subset of the real numbers. These floating-point numbers mean something exact, and therefore mean something different than your phrase, "with the limited precision I have." </pre> <blockquote type="cite"> <pre wrap="">... </pre> <blockquote type="cite"> <blockquote type="cite"> <pre wrap="">Sure, this is not a common use case, but I would prefer that the software let me do this, and leave it up to me to know what I'm doing. </pre> </blockquote> <pre wrap=""> I would prefer that software didn't try to claim the indefensible. ... </pre> </blockquote> </blockquote> <pre wrap=""> I have use floating-point numbers and FPU's in situations where the correctness of the calculations depended upon the floating-point numbers representing specific real numbers exactly. I had to be very careful doing this to get it right, but I was careful and the calculations were correct. What more can I say? Such use cases exist. I have found it convenient (although hardly essential) to sometimes print those numbers out in decimal (because decimal representations are more familiar to me), and it was therefore helpful that my formatting routine didn't prevent me from doing this exactly. This is exactly my argument for why the behavior of gdtoa that you object to is good in certain specialized instances. This may be a use case that you have never encountered, but nevertheless, there it is. </pre> <blockquote type="cite"> <blockquote type="cite"> <pre wrap="">... 1000000000100000000000000000000000000B Since you have only 11 bits of guaranteed binary precision available, you are making a sweeping assumption about those extra 26 bits; (they must ALL be zero). If you know for certain that this is so, then okay, </pre> </blockquote> </blockquote> <pre wrap=""> That's the point. There are specialized cases where I know with mathematically provable certainty that those extra bits are all zero. Not because they're stored in memory -- they're not -- but because my calculation was specifically structured so that they had to be. </pre> <blockquote type="cite"> <blockquote type="cite"> <pre wrap="">but since you don't have those bits available, you have no technically defensible basis, in the general case, for making such an assumption; your argument is flawed, and IMO technically invalid. </pre> </blockquote> </blockquote> <pre wrap=""> On the contrary, even though the bits are not available (i.e. are not stored in memory), it is still possible in my specific case (not your general case) to know (not assume) that these bits are zero. My technically defensible basis for knowing this is that my calculation is structured so that these bits being zero is an invariant of the calculation. </pre> </blockquote> <br> This puzzles me! Do you mean testing of numerical routines? Can't see how general calculations can make use of machine numbers...<br> <br> <blockquote cite="mid:BANLkTi=aMb...@ma..." type="cite"> <pre wrap=""> Look, I've never argued or implied that this sort of use case is in any way common, and I stated explicitly at the beginning of the discussion that this sort of use case is atypical. Your comments are quite correct for the way floating-point numbers are used the vast majority of the time. I'm just pointing out that there are some unusual use cases, and that I think it's good that gdtoa supports them. </pre> <blockquote type="cite"> <pre wrap="">... Chuck </pre> </blockquote> <pre wrap=""> Best regards. K. Frank </pre> </blockquote> </body> </html> |
From: K. F. <kfr...@gm...> - 2011-05-02 21:13:38
|
Hi Peter! On Mon, May 2, 2011 at 3:42 PM, Peter Rockett <p.r...@sh...> wrote: > Forgive the top-post but this whole disagreement seems concerned with > machine numbers, I like your term "machine numbers." I hadn't heard it before. I believe this is what I've been talking about. > that small set of floats that can be represented exactly: > the trailing zeroes beyond the least-significant mantissa bit are indeed > zeroes Yes, I'm talking about unusual use cases where the calculations work with that special set of real numbers that can be represented exactly by floating-point numbers. (Unless I misunderstand, these are your "machine numbers.") > Machine numbers have a important role in numerical algebra since > they are zero-error numbers and hence probe the intrinsic accuracy of a > routine rather than routine + error of input quantity convolved together in > some way. But as far as I am aware, machine numbers are only ever used for > testing numerical routines - since you cannot say that a general float is a > machine number or just the nearest approximation to something nearby on the > real number line, I am puzzled by the points annotated at the bottom of this > post. Or do these refer only to testing? Testing of numerical routines -- that's a good point. This is a use case I wasn't aware of before, but it certainly makes sense. > P. > > On 02/05/2011 18:17, K. Frank wrote: > > Hello Charles and Keith! > On Mon, May 2, 2011 at 11:30 AM, Charles Wilson wrote: > On 5/2/2011 2:34 AM, Keith Marshall wrote: > On 01/05/11 14:36, K. Frank wrote: > [snip] > ... >>> but since you don't have those bits available, you have no technically >>> defensible basis, in the general case, for making such an assumption; >>> your argument is flawed, and IMO technically invalid. >> >> On the contrary, even though the bits are not available (i.e. are not stored >> in memory), it is still possible in my specific case (not your general case) >> to know (not assume) that these bits are zero. My technically defensible >> basis for knowing this is that my calculation is structured so that these >> bits being zero is an invariant of the calculation. > > This puzzles me! Do you mean testing of numerical routines? Can't see how > general calculations can make use of machine numbers... I agree that general calculations don't make use of machine numbers. My point is that there are certain specialized calculations that do. You mentioned a use case new to me -- testing of numerical routines. Some others are: Interval arithmetic. in which a floating-point number is an exact representation of a real number that is a strict upper (or lower) bound to a quantity being studied. As an optimization in arbitrary-precision arithmetic: There was a symbolic algebra package (I believe that it was Mathematica, but it may have been one of its predecessor) that used floating-point numbers to exactly represent real numbers, and then cut over to a full-blown arbitrary precision package when the result of a calculation required it. (As I recall, this optimization was only used in the first few version of the package. I believe that it was never bug free, presumably because it was too hard to detect the necessary cut-over point correctly while still achieving useful gains in performance.) The use case I dealt with directly was that of implementing a random-number generator in the stack of an 80287 FPU. (I don't remember for sure, but I believe it was a linear-congruential generator.) The mathematics that gives you good quality random numbers requires that the calculations be carried out exactly, not just to some finite precision, that is, that they be carried out solely within the set of machine numbers. This is all I've been saying -- that there are relatively rare, but legitimate use cases where you use floating-point numbers as machine numbers, that is to represent specific real numbers exactly. It's hardly a big deal, but in these cases it's nice to be able to print these numbers out exactly in decimal. Best. K. Frank |
From: Kai T. <kti...@go...> - 2011-05-03 10:58:19
|
2011/5/3 Peter Rockett <p.r...@sh...>: > On 03/05/2011 10:11, Kai Tietz wrote: >> 2011/5/3 JonY<jo...@us...>: >>> On 5/3/2011 04:05, Kai Tietz wrote: >>>> 2011/5/2 Peter Rockett<p.r...@sh...>: >>>>> On 02/05/2011 17:00, Keith Marshall wrote: >>>>>> On 02/05/11 16:30, Charles Wilson wrote: >>>>>>>> So, we would like sqrt (x) and pow (x, 0.5) to agree. >>>>> The point below about off-topic diversions into precision is valid but >>>>> can I reiterate, there is no good reason why x^(0.5) calculated two >>>>> different ways should agree since the two different routes will >>>>> accumulate different rounding errors. Adding a special case to force >>>>> agreement for x^n when n = 0.5 will introduce a discontinuity into the >>>>> pow() function and is a very bad idea. The value of returned by pow(x, >>>>> 0.5) might be less accurate compared to sqrt(x) but that is part of the >>>>> fundamental issue with using floating-point numbers: they are only ever >>>>> approximations. You have to take that fact on board. >>>> Correct, the ways of calculation of pow (x, 0.5) and sqrt (x) are different. >>>> Nevertheless is it an issue in gcc. As if you are writing pow (1.1, >>>> 0.5) or sqrt (1.1) - which means gcc optimize it via gmp to a >>>> precalculated result - it has identical values. But by using >>>> math-library results in a - expectable - difference in result. >>>> So IMHO pow should special case here to provide same result as pow >>>> does for y = 0.5. >>> I would agree with the others that its a horrible hack to treat "0.5" >>> special, I don't remember reading anything about requiring pow() results >>> to agree with sqrt(). >>> >>> Can we at least follow a standard spec like: >>> <http://pubs.opengroup.org/onlinepubs/009695399/functions/pow.html> >>> >>> We can remove the special case, right? >> Well, as some people aren't able to use reply-all here, but have wise >> suggestion about implementation details they don't satisfy themself, I >> am a bit annoyed. >> Anyway, I see that this special-casing is neither explicitly demanded, >> but also it isn't explicitly forbidden. So I see here not much reason >> for removing this. > > 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. > > Sorry, but that's a really bad call IMO for the reason stated above. You > are ignoring the consequences to try to (mistakenly) enforce > mathematical correctness. Even basic FP ops like addition and > subtraction deviate from pure mathematical notions of equality. For > example, 1 + a, where a is less than the machine epsilon will return > precisely unity! Floating point numbers ain't reals and you can't force > them to behave exactly like reals... > > P. 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. >> Well, it might be a good idea on mingw.org's side to recheck the exp2, >> exp2f, exp2l, and expl implementation, as here you can easily prove >> that the calculated results have a rounding issue, which leads for >> combined operations - like exp2 (y * log2 (fabs (x)) - to >> mis-calculations over two eps. But well, this might be treated as >> typical inaccuracy ... >> >> Regards, >> Kai Regards, Kai |
From: Peter R. <p.r...@sh...> - 2011-05-03 12:10:52
|
On 03/05/2011 11:58, Kai Tietz wrote: > 2011/5/3 Peter Rockett<p.r...@sh...>: >> On 03/05/2011 10:11, Kai Tietz wrote: >>> 2011/5/3 JonY<jo...@us...>: >>>> On 5/3/2011 04:05, Kai Tietz wrote: >>>>> 2011/5/2 Peter Rockett<p.r...@sh...>: >>>>>> On 02/05/2011 17:00, Keith Marshall wrote: >>>>>>> On 02/05/11 16:30, Charles Wilson wrote: >>>>>>>>> So, we would like sqrt (x) and pow (x, 0.5) to agree. >>>>>> The point below about off-topic diversions into precision is valid but >>>>>> can I reiterate, there is no good reason why x^(0.5) calculated two >>>>>> different ways should agree since the two different routes will >>>>>> accumulate different rounding errors. Adding a special case to force >>>>>> agreement for x^n when n = 0.5 will introduce a discontinuity into the >>>>>> pow() function and is a very bad idea. The value of returned by pow(x, >>>>>> 0.5) might be less accurate compared to sqrt(x) but that is part of the >>>>>> fundamental issue with using floating-point numbers: they are only ever >>>>>> approximations. You have to take that fact on board. >>>>> Correct, the ways of calculation of pow (x, 0.5) and sqrt (x) are different. >>>>> Nevertheless is it an issue in gcc. As if you are writing pow (1.1, >>>>> 0.5) or sqrt (1.1) - which means gcc optimize it via gmp to a >>>>> precalculated result - it has identical values. But by using >>>>> math-library results in a - expectable - difference in result. >>>>> So IMHO pow should special case here to provide same result as pow >>>>> does for y = 0.5. >>>> I would agree with the others that its a horrible hack to treat "0.5" >>>> special, I don't remember reading anything about requiring pow() results >>>> to agree with sqrt(). >>>> >>>> Can we at least follow a standard spec like: >>>> <http://pubs.opengroup.org/onlinepubs/009695399/functions/pow.html> >>>> >>>> We can remove the special case, right? >>> Well, as some people aren't able to use reply-all here, but have wise >>> suggestion about implementation details they don't satisfy themself, I >>> am a bit annoyed. >>> Anyway, I see that this special-casing is neither explicitly demanded, >>> but also it isn't explicitly forbidden. So I see here not much reason >>> for removing this. >> 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. >> Sorry, but that's a really bad call IMO for the reason stated above. You >> are ignoring the consequences to try to (mistakenly) enforce >> mathematical correctness. Even basic FP ops like addition and >> subtraction deviate from pure mathematical notions of equality. For >> example, 1 + a, where a is less than the machine epsilon will return >> precisely unity! Floating point numbers ain't reals and you can't force >> them to behave exactly like reals... >> >> P. > 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 can see no mention of thus in the standard. Is this part of the mingw64 implementation? >>> Well, it might be a good idea on mingw.org's side to recheck the exp2, >>> exp2f, exp2l, and expl implementation, as here you can easily prove >>> that the calculated results have a rounding issue, which leads for >>> combined operations - like exp2 (y * log2 (fabs (x)) - to >>> mis-calculations over two eps. But well, this might be treated as >>> typical inaccuracy ... >>> >>> Regards, >>> Kai > Regards, > Kai > > ------------------------------------------------------------------------------ > 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-users mailing list > Min...@li... > > This list observes the Etiquette found at > http://www.mingw.org/Mailing_Lists. > We ask that you be polite and do the same. Disregard for the list etiquette may cause your account to be moderated. > > _______________________________________________ > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingw-users > Also: mailto:min...@li...?subject=unsubscribe |
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 |
From: James K B. <bea...@ve...> - 2011-05-03 16:37:55
|
There have been special cases in scientific libraries for quite a long time. The special case for atan(1.0) is used by many programmers to extract the system stored value for pi/4, knowing that that the math library (or co-processor) wouldnt compute it instead of simply retrieving the known result. The simple fact that a function like pow(*,*) will have an accuracy requirement that allows a bit or two of numerical error while the nature of the special case sqrt(*) allows full mantissa accuracy to within a possible toggle of the last bit, is not the same thing. The arc tangent is computed by libraries by expansion about the argument points 0 and 1, and identities like arctan(|x|) = pi/2 - arctan(1/|x|) when |x|>1 and such. But this doesn't apply to pow(*,*) because, other than extraction of the exponents of the two arguments to limit the range of the logarithm and exponential in the identity a^b = exp(b*ln(a)) there is no breaking up of things with a natural boundary at b=0.5, although special cases for b=0 and 1 are often taken, and sometimes checks are made to determine whether or not b is an integer. The usual process is breaking up a into exponent and mantissa, a = (2^ea)*ma then computation of the logarithm of the result, ln(a^b) = lr = b*ln(a) = b*(eb*ln(2)+ln(ma)) = ilr + mlr, |mlr|<1 and finally exponentiation of the log result, a^b = (e^ilr)*(e^mlr). There is no obvious advantage for checking for any particular value of b other than 0 and 1 and possibly integers. Thus it seems to me, the obvious thing to do is to accept the specified accuracy of pow(*,*) and not do special cases other than 0, 1, and possibly integers. One might stretch things a bit by evaluating what we would call pow(a^2, b/2) but that would not likely survive an analysis of speed and accuracy as compared to pow(a,b). James K Beard -----Original Message----- From: K. Frank [mailto:kfr...@gm...] Sent: Tuesday, May 03, 2011 9:32 AM To: mingw64; MinGW Users List Subject: Re: [Mingw-users] [Mingw-w64-public] Math library discrepancies that surprised me. 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 ---------------------------------------------------------------------------- -- 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-users mailing list Min...@li... This list observes the Etiquette found at http://www.mingw.org/Mailing_Lists. We ask that you be polite and do the same. Disregard for the list etiquette may cause your account to be moderated. _______________________________________________ You may change your MinGW Account Options or unsubscribe at: https://lists.sourceforge.net/lists/listinfo/mingw-users Also: mailto:min...@li...?subject=unsubscribe |
From: Greg C. <gch...@sb...> - 2011-04-29 23:38:10
|
On 2011-04-29 19:55Z, K. Frank wrote: > > By the way, could someone confirm that mingw does use msvcrt for > sqrt and pow? http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). 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 IEEE-754-1985 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 correctly-rounded 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 special-case the power of exactly 0.5? If you do, then pow() and sqrt() will be consistent, but then the powers x^(0.5*(1-epsilon)) x^(0.5 ) x^(0.5*(1+epsilon)) might not be monotone. Some would call that a bad tradeoff. |
From: K. F. <kfr...@gm...> - 2011-04-30 00:51:08
|
Hi Greg! On Fri, Apr 29, 2011 at 7:38 PM, Greg Chicares <gch...@sb...> wrote: > On 2011-04-29 19:55Z, K. Frank wrote: >> >> By the way, could someone confirm that mingw does use msvcrt for >> sqrt and pow? > > http://cygwin.com/cgi-bin/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 quality-of-implementation 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/cgi-bin/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.2e-16 7.7e-17 33: * DEC -26,26 60000 4.8e-17 9.1e-18 34: * 1/26 < x < 26, with log(x) uniformly distributed. 35: * -26 < y < 26, y uniformly distributed. 36: * IEEE 0,8700 30000 1.5e-14 2.1e-15 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 IEEE-754). > > 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 floating-point 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 IEEE-754-1985 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 correctly-rounded 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 special-case > the power of exactly 0.5? If you do, then pow() and sqrt() will > be consistent, but then the powers > x^(0.5*(1-epsilon)) > 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 |
From: Kai T. <kti...@go...> - 2011-04-30 08:42:16
|
2011/4/30 Greg Chicares <gch...@sb...>: > On 2011-04-29 19:55Z, K. Frank wrote: >> >> By the way, could someone confirm that mingw does use msvcrt for >> sqrt and pow? > > http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). > > 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 IEEE-754-1985 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 correctly-rounded 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 special-case > the power of exactly 0.5? If you do, then pow() and sqrt() will > be consistent, but then the powers > x^(0.5*(1-epsilon)) > 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 base-math functions I improved in mingw-w64's math library in two terms. a) Make those functions behave as ISO-C99 specification tells in corner-cases. 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 math-calculations gmp-library 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 side-effect. 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 |
From: Kai T. <kti...@go...> - 2011-04-30 09:51:40
|
2011/4/30 Kai Tietz <kti...@go...>: > 2011/4/30 Greg Chicares <gch...@sb...>: >> On 2011-04-29 19:55Z, K. Frank wrote: >>> >>> By the way, could someone confirm that mingw does use msvcrt for >>> sqrt and pow? >> >> http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). >> >> 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 IEEE-754-1985 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 correctly-rounded 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 special-case >> the power of exactly 0.5? If you do, then pow() and sqrt() will >> be consistent, but then the powers >> x^(0.5*(1-epsilon)) >> 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 base-math functions I improved in mingw-w64's math library > in two terms. a) Make those functions behave as ISO-C99 specification > tells in corner-cases. 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 > math-calculations gmp-library 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 side-effect. > > 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 > On closer look for long double, we can have rouding issues for long double precission. I see that powl (0.01L, 0.5L) we have a rounding difference to sqrtl (0.01L). Interestingly not for float or double types. So we we might need to special-case in pow the case for y == 0.5 to solve this. Kai |
From: K. F. <kfr...@gm...> - 2011-04-30 15:15:19
|
Hello Kai! On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <kti...@go...> wrote: > 2011/4/30 Kai Tietz <kti...@go...>: >> 2011/4/30 Greg Chicares <gch...@sb...>: >>> On 2011-04-29 19:55Z, K. Frank wrote: >>>> >>>> By the way, could someone confirm that mingw does use msvcrt for >>>> sqrt and pow? >>> >>> http://cygwin.com/cgi-bin/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 quality-of-implementation 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 IEEE-754). >>> ... >>> Not so many years ago (perhaps when msvcrt was written), it was >>> thought that correctly-rounded 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 special-case >>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>> be consistent, but then the powers >>> x^(0.5*(1-epsilon)) >>> 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 base-math functions I improved in mingw-w64's math library >> in two terms. a) Make those functions behave as ISO-C99 specification >> tells in corner-cases. 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; >> } Do I understand your point correctly that you are comparing gcc's compile-time and run-time values for pow? For the record, when I run this on a mingw64 build: gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] I get: C:\>gcc -o pow_test pow_test.c C:\>pow_test 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 11.52153412678571697825 8.00000000000000000000 2.00000000000000000000 0.00000000000000177636 0.00000000000000000000 0.00000000000000000000 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000 4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000 11.52153412678571414718 8.00000000000000000000 2.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 As I understand it, this means I am seeing a small difference between compile and run-time evaluation only for powl (2.2, 3.1), i.e., the long-double case. >> The interesting issue is here that gcc uses for constant >> math-calculations gmp-library 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. Yes, that's a good point. I hadn't considered the issue of compile-time evaluation, and, ideally, it should agree with run-time evaluation, adding yet another constraint to our wish list. Furthermore, if we take the idea of a cross-compiler seriously, we would also like the compile-time evaluations of a cross-compiler and native compiler to agree. Hmm... So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like compile-time and run-time evaluations to agree. We would like cross-compilers and native compilers to agree. Taken at face value, this would seem to preclude using the platform library and platform hardware to compute things like sqrt and pow (unless we know they comply with some agreed-upon standard, such as IEEE-754). Just to make the reasoning explicit, if I perform my run-time evaluation on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler running on platform B won't have platform A's FPU available to perform consistent compile-time evaluations. So you either give up some of the above-wished-for consistency, you forgo the potential performance benefits of using the platform's library and / or hardware, or you get everybody to use a common standard. >> ... >> 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 > > On closer look for long double, we can have rouding issues for long > double precission. I see that powl (0.01L, 0.5L) we have a rounding > difference to sqrtl (0.01L). Interestingly not for float or double > types. So we we might need to special-case in pow the case for y == > 0.5 to solve this. Be, however, cognizant of Greg's concern. Forcing pow to agree with sqrt just for y = 0.5 has the potential to make pow, in a sense, internally inconsistent, It's a tricky trade-off. (Ideally, both pow and sqrt would be "exact" and therefore agree with one another automatically.) By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L). It disagrees slightly from sqrtl (2.2L) and from the compile-time value for powl (2.2L, 0.5L). I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl with powl). For the record: C:\>gcc --version gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz] C:\>gcc -o sqrt_pow sqrt_pow.c C:\>sqrt_pow 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000 1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000 1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000 1.04880884817015163080 1.04880884817015163080 1.48323969741913264109 1.48323969741913264109 1.81659021245849494619 1.81659021245849494619 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000 1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000 1.04880884817015154699 1.04880884817015154699 1.48323969741913258980 1.48323969741913258980 1.81659021245849499921 1.81659021245849499921 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000 where sqrt_pow.c is: #include <stdio.h> #include <math.h> double abc[3] = { 1.1, 2.2, 3.3 }; long double abcl[3] = { 1.1L, 2.2L, 3.3L }; int main() { double r1[3], r2[3], r3[3]; long double rl1[3], rl2[3], rl3[3]; int i; for (i = 0; i < 3; i++) __mingw_printf ("%.20F, %.20F, %.20F\n", (r1[i] = sqrt (abc[i])), (r2[i] = pow (abc[i], 0.5)), (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i]))); __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5)); __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5)); __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5)); r1[0] -= sqrt (1.1); r1[1] -= sqrt (2.2); r1[2] -= sqrt (3.3); __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]); r2[0] -= pow (1.1, 0.5); r2[1] -= pow (2.2, 0.5); r2[2] -= pow (3.3, 0.5); __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]); for (i = 0; i < 3; i++) __mingw_printf ("%.20LF, %.20LF, %.20LF\n", (rl1[i] = sqrtl (abcl[i])), (rl2[i] = powl (abcl[i], 0.5L)), (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i]))); __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L)); __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L)); __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L)); rl1[0] -= sqrtl (1.1L); rl1[1] -= sqrtl (2.2L); rl1[2] -= sqrtl (3.3L); __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]); rl2[0] -= powl (1.1L, 0.5L); rl2[1] -= powl (2.2L, 0.5L); rl2[2] -= powl (3.3L, 0.5L); __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]); return 1; } > > Kai Thanks, Kai, for your insight and explanations. K. Frank |
From: Bob D. <del...@ea...> - 2011-04-30 21:21:26
|
On Apr 30, 2011, at 10:15 AM, K. Frank wrote: > I get: > > C:\>gcc -o pow_test pow_test.c > > C:\>pow_test > 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 Victor Shoup's NTL and my own fp, to 22 digit accuracy, give: 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832 so your result is accurate to only 16 digits. Your 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 should be: 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732 so your result is accurate to about 20 digits. Bob |
From: Kai T. <kti...@go...> - 2011-04-30 21:42:03
|
2011/4/30 Bob Delaney <del...@ea...>: > > On Apr 30, 2011, at 10:15 AM, K. Frank wrote: > >> I get: >> >> C:\>gcc -o pow_test pow_test.c >> >> C:\>pow_test >> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 > > > Victor Shoup's NTL and my own fp, to 22 digit accuracy, give: > > 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832 > > so your result is accurate to only 16 digits. Your > > 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 > > should be: > > 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732 > > so your result is accurate to about 20 digits. > > Bob Well, what is to be expected as we used here in test application an %.20F. By specifying here a %.40F its calculated result is 11.5215341267857141471786519559827866032720 and it is has binary identity for what gmp calculates for 80-bit IEEE 754 floating point. Which means for our purposes in gcc, it has right accuracy. Kai |
From: K. F. <kfr...@gm...> - 2011-04-30 22:13:15
|
Hello Bob and Kai! On Sat, Apr 30, 2011 at 5:41 PM, Kai Tietz <kti...@go...> wrote: > 2011/4/30 Bob Delaney <del...@ea...>: >> >> On Apr 30, 2011, at 10:15 AM, K. Frank wrote: >> >>> I get: >>> >>> C:\>gcc -o pow_test pow_test.c >>> >>> C:\>pow_test >>> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 >> >> >> Victor Shoup's NTL and my own fp, to 22 digit accuracy, give: >> >> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832 >> >> so your result is accurate to only 16 digits. Your This is exactly what I would expect (short of looking at the actual binary representation). An IEEE double is 64 bits, with a 52-bit mantissa, which, with the "implicit" bit, gives 53 bits of precision. This gives just under 16 decimal digits of precision, which is what you see. >> >> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 >> >> should be: >> >> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732 >> >> so your result is accurate to about 20 digits. As I understand it, gcc defines long double as the x87 80-bit floating point number. This has a 64-bit mantissa, and I believe that it does not use an implicit bit, so this gives slightly more than 19 decimal digits of precision, again consistent with your result. >> >> Bob > > Well, what is to be expected as we used here in test application an > %.20F. By specifying here a %.40F its calculated result is > 11.5215341267857141471786519559827866032720 and it is has binary > identity for what gmp calculates for 80-bit IEEE 754 floating point. > Which means for our purposes in gcc, it has right accuracy. (Kai's comment seems to confirm that gcc's long double is, indeed, the 80-bit floating-point number.) > Kai Thanks to all for the insight. K. Frank |
From: Kai T. <kti...@go...> - 2011-04-30 22:14:14
|
2011/4/30 Kai Tietz <kti...@go...>: > 2011/4/30 Bob Delaney <del...@ea...>: >> >> On Apr 30, 2011, at 10:15 AM, K. Frank wrote: >> >>> I get: >>> >>> C:\>gcc -o pow_test pow_test.c >>> >>> C:\>pow_test >>> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460 >> >> >> Victor Shoup's NTL and my own fp, to 22 digit accuracy, give: >> >> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832 >> >> so your result is accurate to only 16 digits. Your >> >> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718 >> >> should be: >> >> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732 >> >> so your result is accurate to about 20 digits. >> >> Bob > > Well, what is to be expected as we used here in test application an > %.20F. By specifying here a %.40F its calculated result is > 11.5215341267857141471786519559827866032720 and it is has binary > identity for what gmp calculates for 80-bit IEEE 754 floating point. > Which means for our purposes in gcc, it has right accuracy. > > Kai > To be more precise here IEEE 754 floats for Intel/Amd architecture x86/x86_64 has following definitions: float (32-bit): digits for mantissa:24 digit:6, min exp:-125 (-37 10th) max exp: 128 (38 10th) Epsilon:1.19209e-007 double (64-bit): digits for mantissa:53 digit:15, min exp:-1021 (-307 10th) max exp: 1024 (308 10th) Epsilon:2.22045e-016 long double (80-bit): digits for mantissa:64 digit:18, min exp:-16381 (-4931 10th) max exp: 16384 (4932 10th) Epsilon:1.0842e-019 Which means that indeed the accuracy of 20 digits after decimal point are pretty precise here. Regards, Kai |
From: Keith M. <kei...@us...> - 2011-05-01 07:43:46
|
On 30/04/11 23:14, Kai Tietz wrote: > long double (80-bit): > digits for mantissa:64 > digit:18, min exp:-16381 (-4931 10th) > max exp: 16384 (4932 10th) > Epsilon:1.0842e-019 > > Which means that indeed the accuracy of 20 digits after decimal point > are pretty precise here. You cannot get 20 decimal digits precision after the decimal point, with 80-bit floating point. You cannot even get 20 significant digit precision -- there is a subtle but extremely important distinction. The best you can hope for is 19 significant digits. Alas, this sort of mathematical ineptitude seems to be all too common amongst the programming fraternity. It isn't helped by a flaw in the gdtoa implementation common on BSD and GNU/Linux, (and also adopted by Danny into mingwrt); it fails to implement the Steele and White stop condition correctly, continuing to spew out garbage beyond the last bit of available precision, creating an illusion of better accuracy than is really achievable. -- Regards, Keith. |
From: Peter R. <p.r...@sh...> - 2011-05-01 11:01:13
|
On 01/05/11 08:43, Keith Marshall wrote: > On 30/04/11 23:14, Kai Tietz wrote: >> long double (80-bit): >> digits for mantissa:64 >> digit:18, min exp:-16381 (-4931 10th) >> max exp: 16384 (4932 10th) >> Epsilon:1.0842e-019 >> >> Which means that indeed the accuracy of 20 digits after decimal point >> are pretty precise here. I've followed with some unease at the way "accuracy" and "precision" have been used in this thread. To nuance Keith's point below, you can write down 20 digits after the decimal point but that does not mean you have 20-digit accuracy. They are very different things. > You cannot get 20 decimal digits precision after the decimal point, with > 80-bit floating point. You cannot even get 20 significant digit > precision -- there is a subtle but extremely important distinction. The > best you can hope for is 19 significant digits. So long doubles don't do normalisation, and thus have no extra, implicit mantissa bit? (Have come across something that suggests this as an historical quirk but it was unclear.) That would set long doubles apart from singles and doubles in the IEEE754 scheme of things... Interesting. To add my two pennyworth on pow/sqrt: Although changing the behaviour of pow() to give the same answer as sqrt() for the special case of 0.5 may seem like a good idea, this is potentially introducing a discontinuity into pow() which could be disastrous. Fundamentally, why should you expect two transcendental quantities calculated in different ways with finite-precision arithmetic to give exactly the same result? Floats are not real numbers! In fact some of the comparisons in this thread have involved differences of a few ulps (units of least precision) - which looks really good to me. > Alas, this sort of mathematical ineptitude seems to be all too common > amongst the programming fraternity. It isn't helped by a flaw in the > gdtoa implementation common on BSD and GNU/Linux, (and also adopted by > Danny into mingwrt); it fails to implement the Steele and White stop > condition correctly, continuing to spew out garbage beyond the last bit > of available precision, creating an illusion of better accuracy than is > really achievable. Can't answer for the gdtoa implementation but surely this is only a problem if you're naive enough to actually believe all the digits! If your algorithm relies on getting accuracy down to the ulp level for transcendentals then you really should be googling "ill-conditioning". P. |