Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: K. Frank <kfrank29.c@gm...>  20110429 19:56:00

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

My comments don't seem to be making it into the thread this week, possibly because of a switch between my replyto 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 reducedrange mantissa in a modified Chebychev approximation. Arithmetic coprocessors do sqrt(), ln(), and exp() internally with 80bit floating point, then store in 32, 64, or 80 bit floating point format. So, the coprocessors can get full mantissa accuracy with 32bit and 64bit sqrt() and pow(x,0.5). So, the problem is nonexistent 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:gchicares@...] Sent: Friday, April 29, 2011 7:38 PM To: MinGW Users List Cc: mingw64 Subject: Re: [Mingww64public] [Mingwusers] Math library discrepancies that surprised me. On 20110429 19:55Z, K. Frank wrote: > > By the way, could someone confirm that mingw does use msvcrt for > sqrt and pow? http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot= src Looks to me like sqrt(double) calls into msvcrt. But pow() is a different story. Danny implemented it in libmingwex because of some qualityofimplementation problem with msvcrt; then people complained that the libmingwex version was slower, and I can't remember where it wound up, but it's all in the archives. > sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required > to by IEEE754). If they should give exactly the same result always, then the first thing that needs to be fixed is any standard that allows them to differ. The language standards are so much looser than the numerical standard, even now, that I'd rather see C and C++ catch up to IEEE7541985 before asking them to go beyond. > Browsing some gcc lists, I did see some comments that suggest that > gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x). This would > make sqrt and pow consistent (whether or not technically correct). Not so many years ago (perhaps when msvcrt was written), it was thought that correctlyrounded transcendental functions weren't feasible in practice, so library authors had a looser attitude toward exactness. If you have an ideal implementation of sqrt() but a pow() that's only approximate, should you specialcase the power of exactly 0.5? If you do, then pow() and sqrt() will be consistent, but then the powers x^(0.5*(1epsilon)) x^(0.5 ) x^(0.5*(1+epsilon)) might not be monotone. Some would call that a bad tradeoff.   WhatsUp Gold  Download Free Network Management Software The most intuitive, comprehensive, and costeffective network management toolset available today. Delivers lowest initial acquisition cost and overall TCO of any competing solution. http://p.sf.net/sfu/whatsupgoldsd _______________________________________________ Mingww64public mailing list Mingww64public@... https://lists.sourceforge.net/lists/listinfo/mingww64public 
From: Kai Tietz <ktietz70@go...>  20110430 15:34:00

2011/4/30 K. Frank <kfrank29.c@...>: > Hello Kai! > > On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote: >> 2011/4/30 Kai Tietz <ktietz70@...>: >>> 2011/4/30 Greg Chicares <gchicares@...>: >>>> On 20110429 19:55Z, K. Frank wrote: >>>>> >>>>> By the way, could someone confirm that mingw does use msvcrt for >>>>> sqrt and pow? >>>> >>>> http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src >>>> >>>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a >>>> different story. Danny implemented it in libmingwex because of >>>> some qualityofimplementation problem with msvcrt; then people >>>> complained that the libmingwex version was slower, and I can't >>>> remember where it wound up, but it's all in the archives. >>>> >>>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >>>>> to by IEEE754). >>>> ... >>>> Not so many years ago (perhaps when msvcrt was written), it was >>>> thought that correctlyrounded transcendental functions weren't >>>> feasible in practice, so library authors had a looser attitude >>>> toward exactness. If you have an ideal implementation of sqrt() >>>> but a pow() that's only approximate, should you specialcase >>>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>>> be consistent, but then the powers >>>> x^(0.5*(1epsilon)) >>>> x^(0.5 ) >>>> x^(0.5*(1+epsilon)) >>>> might not be monotone. Some would call that a bad tradeoff. >>> >>> Well, I spent some efforts into log, exp, and pow implementation. Also >>> some other basemath functions I improved in mingww64's math library >>> in two terms. a) Make those functions behave as ISOC99 specification >>> tells in cornercases. and b) Improve accuracy in terms of gcc's >>> internal used gmp. As we noticed here result differences. >>> >>> The following example illustrates the issues pretty well: >>> >>> #include <stdio.h> >>> #include <math.h> >>> >>> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; >>> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; >>> int main() >>> { >>> double r[3]; >>> long double rl[3]; >>> int i; >>> for (i = 0; i < 3; i++) >>> __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], >>> (r[i] = pow (abc[i][0], abc[i][1]))); >>> __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, >>> 3.0), pow (4.0, 0.5)); >>> r[0] = pow (2.2, 3.1); >>> r[1] = pow (2.0, 3.0); >>> r[2] = pow (4.0, 0.5); >>> __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); >>> >>> for (i = 0; i < 3; i++) >>> __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], >>> (rl[i] = powl (ab[i][0], ab[i][1]))); >>> __mingw_printf ("%.20LF %.20LF %.20LF\n", >>> powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl >>> (4.0L, 0.5L)); >>> rl[0] = powl (2.2L, 3.1L); >>> rl[1] = powl (2.0L, 3.0L); >>> rl[2] = powl (4.0L, 0.5L); >>> __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); >>> return 1; >>> } > > Do I understand your point correctly that you are comparing gcc's > compiletime and runtime values for pow? > > For the record, when I run this on a mingw64 build: > > gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902  mingww64/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 runtime evaluation only for powl (2.2, 3.1), i.e., the > longdouble case. > >>> The interesting issue is here that gcc uses for constant >>> mathcalculations gmplibrary to get result. For more complex >>> variants, it uses crt's math. By this it is important for a gcc based >>> runtime to be in the IEEE 754 floating point compatible to gmp's >>> results. > > Yes, that's a good point. I hadn't considered the issue of compiletime > evaluation, and, ideally, it should agree with runtime evaluation, adding > yet another constraint to our wish list. > > Furthermore, if we take the idea of a crosscompiler seriously, we > would also like the compiletime evaluations of a crosscompiler and > native compiler to agree. > > Hmm... > > So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like > compiletime and runtime evaluations to agree. We would like > crosscompilers 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 agreedupon standard, > such as IEEE754). > > Just to make the reasoning explicit, if I perform my runtime evaluation > on platform A, e.g., using platform A's FPU, then a B>A crosscompiler > running on platform B won't have platform A's FPU available to perform > consistent compiletime evaluations. > > So you either give up some of the abovewishedfor 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 specialcase 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 tradeoff. (Ideally, both pow and sqrt would > be "exact" and therefore agree with one another automatically.) > > By the way, I also see a problem with runtime powl for powl (2.2L, 0.5L). > It disagrees slightly from sqrtl (2.2L) and from the compiletime 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  mingww64/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 mingww64's pow(fl) 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 32bit variant and 64bit 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 shortpath calculation in an ISOC variant (as gcc expects internally such a math). Also complexmath part is completely rewritten on mingww64'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 importsymbols 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 Tietz <ktietz70@go...>  20110430 15:48:06

2011/4/30 Kai Tietz <ktietz70@...>: > 2011/4/30 K. Frank <kfrank29.c@...>: >> Hello Kai! >> >> On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote: >>> 2011/4/30 Kai Tietz <ktietz70@...>: >>>> 2011/4/30 Greg Chicares <gchicares@...>: >>>>> On 20110429 19:55Z, K. Frank wrote: >>>>>> >>>>>> By the way, could someone confirm that mingw does use msvcrt for >>>>>> sqrt and pow? >>>>> >>>>> http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src >>>>> >>>>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a >>>>> different story. Danny implemented it in libmingwex because of >>>>> some qualityofimplementation problem with msvcrt; then people >>>>> complained that the libmingwex version was slower, and I can't >>>>> remember where it wound up, but it's all in the archives. >>>>> >>>>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >>>>>> to by IEEE754). >>>>> ... >>>>> Not so many years ago (perhaps when msvcrt was written), it was >>>>> thought that correctlyrounded transcendental functions weren't >>>>> feasible in practice, so library authors had a looser attitude >>>>> toward exactness. If you have an ideal implementation of sqrt() >>>>> but a pow() that's only approximate, should you specialcase >>>>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>>>> be consistent, but then the powers >>>>> x^(0.5*(1epsilon)) >>>>> x^(0.5 ) >>>>> x^(0.5*(1+epsilon)) >>>>> might not be monotone. Some would call that a bad tradeoff. >>>> >>>> Well, I spent some efforts into log, exp, and pow implementation. Also >>>> some other basemath functions I improved in mingww64's math library >>>> in two terms. a) Make those functions behave as ISOC99 specification >>>> tells in cornercases. and b) Improve accuracy in terms of gcc's >>>> internal used gmp. As we noticed here result differences. >>>> >>>> The following example illustrates the issues pretty well: >>>> >>>> #include <stdio.h> >>>> #include <math.h> >>>> >>>> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; >>>> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; >>>> int main() >>>> { >>>> double r[3]; >>>> long double rl[3]; >>>> int i; >>>> for (i = 0; i < 3; i++) >>>> __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], >>>> (r[i] = pow (abc[i][0], abc[i][1]))); >>>> __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, >>>> 3.0), pow (4.0, 0.5)); >>>> r[0] = pow (2.2, 3.1); >>>> r[1] = pow (2.0, 3.0); >>>> r[2] = pow (4.0, 0.5); >>>> __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); >>>> >>>> for (i = 0; i < 3; i++) >>>> __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], >>>> (rl[i] = powl (ab[i][0], ab[i][1]))); >>>> __mingw_printf ("%.20LF %.20LF %.20LF\n", >>>> powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl >>>> (4.0L, 0.5L)); >>>> rl[0] = powl (2.2L, 3.1L); >>>> rl[1] = powl (2.0L, 3.0L); >>>> rl[2] = powl (4.0L, 0.5L); >>>> __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); >>>> return 1; >>>> } >> >> Do I understand your point correctly that you are comparing gcc's >> compiletime and runtime values for pow? >> >> For the record, when I run this on a mingw64 build: >> >> gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902  mingww64/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 runtime evaluation only for powl (2.2, 3.1), i.e., the >> longdouble case. >> >>>> The interesting issue is here that gcc uses for constant >>>> mathcalculations gmplibrary to get result. For more complex >>>> variants, it uses crt's math. By this it is important for a gcc based >>>> runtime to be in the IEEE 754 floating point compatible to gmp's >>>> results. >> >> Yes, that's a good point. I hadn't considered the issue of compiletime >> evaluation, and, ideally, it should agree with runtime evaluation, adding >> yet another constraint to our wish list. >> >> Furthermore, if we take the idea of a crosscompiler seriously, we >> would also like the compiletime evaluations of a crosscompiler and >> native compiler to agree. >> >> Hmm... >> >> So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like >> compiletime and runtime evaluations to agree. We would like >> crosscompilers 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 agreedupon standard, >> such as IEEE754). >> >> Just to make the reasoning explicit, if I perform my runtime evaluation >> on platform A, e.g., using platform A's FPU, then a B>A crosscompiler >> running on platform B won't have platform A's FPU available to perform >> consistent compiletime evaluations. >> >> So you either give up some of the abovewishedfor 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 specialcase 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 tradeoff. (Ideally, both pow and sqrt would >> be "exact" and therefore agree with one another automatically.) >> >> By the way, I also see a problem with runtime powl for powl (2.2L, 0.5L). >> It disagrees slightly from sqrtl (2.2L) and from the compiletime 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  mingww64/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 mingww64's pow(fl) 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 32bit variant and 64bit 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 > shortpath calculation in an ISOC variant (as gcc expects internally > such a math). Also complexmath part is completely rewritten on > mingww64'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 importsymbols 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 32bit aren't caused by inaccuarcy of used math directly. Instead it is caused by the fact that 32bit use FPU registers (which are 80bit 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 32bit a delta of zero. Interesting this behavior is on 64bit x86_64 different, as it uses always precise accuracy to store results in registers. By this you have really just a 64bit result for double type, and not indirecty 80bit (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 <jon_y@us...>  20110501 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 32bit and 64bit 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 90s200Xs era) don't have 80bit FPUs, but they do have 64bit 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.08420217248550443401e19L accuracy absolutely require an 80bit 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. Frank <kfrank29.c@gm...>  20110501 15:28:53

Hello Jon! On Sun, May 1, 2011 at 10:17 AM, JonY <jon_y@...> 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 90s200Xs era) don't have 80bit FPUs, but they > do have 64bit 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.08420217248550443401e19L accuracy > absolutely require an 80bit FPU? If your willing to do your floatingpoint math in software, then no. But if you want to get 19decimaldigit 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 80bit floatingpoint format has a 64bit 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 80bit format, and it would be hard to justify making a nonstandard 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 floatingpoint you support in software, and use floatingpoint hardware when available. If you want to do heavy floatingpoint computation, it's up to you to make sure you have the necessary floatingpoint hardware or eat the (significant) performance hit. Nowadays I would think that floatingpoint 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 Chicares <gchicares@sb...>  20110501 15:11:35

On 20110501 13:38Z, K. Frank wrote: [...] > I'm speaking from memory, but  assuming that gcc long doubles are > x87 80bit floatingpoint numbers  yes, they do not use an implicit bit. Yes, MinGW gcc's 'long double' is intel's 'extended precision' 80bit floating point. BTW, IIRC, sizeof(long double) is 12for reasons of alignment, I guessbut only ten bytes are used. Indeed 80bit 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. Frank <kfrank29.c@gm...>  20110501 15:29:54

Hello Greg! Thanks for clearing this up. On Sun, May 1, 2011 at 11:11 AM, Greg Chicares <gchicares@...> wrote: > On 20110501 13:38Z, K. Frank wrote: > [...] >> I'm speaking from memory, but  assuming that gcc long doubles are >> x87 80bit floatingpoint numbers  yes, they do not use an implicit bit. > > Yes, MinGW gcc's 'long double' is intel's 'extended precision' 80bit > floating point. BTW, IIRC, sizeof(long double) is 12for reasons of > alignment, I guessbut only ten bytes are used. > > Indeed 80bit 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 Rockett <p.rockett@sh...>  20110502 19:42:09

<!DOCTYPE HTML PUBLIC "//W3C//DTD HTML 4.01 Transitional//EN"> <html> <head> <meta content="text/html; charset=ISO88591" httpequiv="ContentType"> </head> <body text="#000000" bgcolor="#ffffff"> Forgive the toppost 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 leastsignificant mantissa bit are indeed zeroes. Machine numbers have a important role in numerical algebra since they are zeroerror 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+QqnLAG_FGPm0J8t6Zp3gcw@..." 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 floatingpoint 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 floatingpoint 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 floatingpoint 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 [valueulps/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 floatingpoint 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 floatingpoint 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 floatingpoint numbers. But, to reiterate my point, there are use cases where floatingpoint numbers are used to exactly represent a special subset of the real numbers. These floatingpoint 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 floatingpoint numbers and FPU's in situations where the correctness of the calculations depended upon the floatingpoint 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+QqnLAG_FGPm0J8t6Zp3gcw@..." 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 floatingpoint 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. Frank <kfrank29.c@gm...>  20110502 21:13:38

Hi Peter! On Mon, May 2, 2011 at 3:42 PM, Peter Rockett <p.rockett@...> wrote: > Forgive the toppost 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 leastsignificant 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 floatingpoint numbers. (Unless I misunderstand, these are your "machine numbers.") > Machine numbers have a important role in numerical algebra since > they are zeroerror 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 floatingpoint 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 arbitraryprecision arithmetic: There was a symbolic algebra package (I believe that it was Mathematica, but it may have been one of its predecessor) that used floatingpoint numbers to exactly represent real numbers, and then cut over to a fullblown 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 cutover point correctly while still achieving useful gains in performance.) The use case I dealt with directly was that of implementing a randomnumber generator in the stack of an 80287 FPU. (I don't remember for sure, but I believe it was a linearcongruential 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 floatingpoint 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 Tietz <ktietz70@go...>  20110503 10:58:19

2011/5/3 Peter Rockett <p.rockett@...>: > On 03/05/2011 10:11, Kai Tietz wrote: >> 2011/5/3 JonY<jon_y@...>: >>> On 5/3/2011 04:05, Kai Tietz wrote: >>>> 2011/5/2 Peter Rockett<p.rockett@...>: >>>>> 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 offtopic 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 floatingpoint 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 >>>> mathlibrary 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 replyall here, but have wise >> suggestion about implementation details they don't satisfy themself, I >> am a bit annoyed. >> Anyway, I see that this specialcasing 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 specialcase 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 softfloatingpoint library. >> As I don't see in spec that this special casing shouldn't be done, I >> think we will keep it on mingww64'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 specialcasing 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 >> miscalculations over two eps. But well, this might be treated as >> typical inaccuracy ... >> >> Regards, >> Kai Regards, Kai 
From: Peter Rockett <p.rockett@sh...>  20110503 12:10:52

On 03/05/2011 11:58, Kai Tietz wrote: > 2011/5/3 Peter Rockett<p.rockett@...>: >> On 03/05/2011 10:11, Kai Tietz wrote: >>> 2011/5/3 JonY<jon_y@...>: >>>> On 5/3/2011 04:05, Kai Tietz wrote: >>>>> 2011/5/2 Peter Rockett<p.rockett@...>: >>>>>> 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 offtopic 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 floatingpoint 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 >>>>> mathlibrary 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 replyall here, but have wise >>> suggestion about implementation details they don't satisfy themself, I >>> am a bit annoyed. >>> Anyway, I see that this specialcasing 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 specialcase 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 softfloatingpoint library. >>> As I don't see in spec that this special casing shouldn't be done, I >>> think we will keep it on mingww64'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 specialcasing 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 >>> miscalculations 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 costeffective network > management toolset available today. Delivers lowest initial > acquisition cost and overall TCO of any competing solution. > http://p.sf.net/sfu/whatsupgoldsd > _______________________________________________ > MinGWusers mailing list > MinGWusers@... > > 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/mingwusers > Also: mailto:mingwusersrequest@...?subject=unsubscribe 
From: K. Frank <kfrank29.c@gm...>  20110503 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.rockett@...>: >> 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 specialcase 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 softfloatingpoint library. >>> As I don't see in spec that this special casing shouldn't be done, I >>> think we will keep it on mingww64's side. >> ... > I am aware of this what you are telling me, but there is in pow > function already specialcasing 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 IEEE754 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 nonmonotonicity 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 nonmonotonicity into pow.) Lastly, I have to believe that this is a wellunderstood and "solved" problem in the numericalanalysis community (of which I am not a part). Perhaps it would make sense to find out how the experts who have doing this for fiftyplus 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 Beard <beardjamesk@ve...>  20110503 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 coprocessor) 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:kfrank29.c@...] Sent: Tuesday, May 03, 2011 9:32 AM To: mingw64; MinGW Users List Subject: Re: [Mingwusers] [Mingww64public] 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.rockett@...>: >> 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 specialcase 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 softfloatingpoint library. >>> As I don't see in spec that this special casing shouldn't be done, I >>> think we will keep it on mingww64's side. >> ... > I am aware of this what you are telling me, but there is in pow > function already specialcasing 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 IEEE754 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 nonmonotonicity 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 nonmonotonicity into pow.) Lastly, I have to believe that this is a wellunderstood and "solved" problem in the numericalanalysis community (of which I am not a part). Perhaps it would make sense to find out how the experts who have doing this for fiftyplus 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 costeffective network management toolset available today. Delivers lowest initial acquisition cost and overall TCO of any competing solution. http://p.sf.net/sfu/whatsupgoldsd _______________________________________________ MinGWusers mailing list MinGWusers@... 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/mingwusers Also: mailto:mingwusersrequest@...?subject=unsubscribe 
From: Greg Chicares <gchicares@sb...>  20110429 23:38:10

On 20110429 19:55Z, K. Frank wrote: > > By the way, could someone confirm that mingw does use msvcrt for > sqrt and pow? http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src Looks to me like sqrt(double) calls into msvcrt. But pow() is a different story. Danny implemented it in libmingwex because of some qualityofimplementation problem with msvcrt; then people complained that the libmingwex version was slower, and I can't remember where it wound up, but it's all in the archives. > sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required > to by IEEE754). If they should give exactly the same result always, then the first thing that needs to be fixed is any standard that allows them to differ. The language standards are so much looser than the numerical standard, even now, that I'd rather see C and C++ catch up to IEEE7541985 before asking them to go beyond. > Browsing some gcc lists, I did see some comments that suggest that > gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x). This would > make sqrt and pow consistent (whether or not technically correct). Not so many years ago (perhaps when msvcrt was written), it was thought that correctlyrounded transcendental functions weren't feasible in practice, so library authors had a looser attitude toward exactness. If you have an ideal implementation of sqrt() but a pow() that's only approximate, should you specialcase the power of exactly 0.5? If you do, then pow() and sqrt() will be consistent, but then the powers x^(0.5*(1epsilon)) x^(0.5 ) x^(0.5*(1+epsilon)) might not be monotone. Some would call that a bad tradeoff. 
From: K. Frank <kfrank29.c@gm...>  20110430 00:51:08

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

2011/4/30 Greg Chicares <gchicares@...>: > On 20110429 19:55Z, K. Frank wrote: >> >> By the way, could someone confirm that mingw does use msvcrt for >> sqrt and pow? > > http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src > > Looks to me like sqrt(double) calls into msvcrt. But pow() is a > different story. Danny implemented it in libmingwex because of > some qualityofimplementation problem with msvcrt; then people > complained that the libmingwex version was slower, and I can't > remember where it wound up, but it's all in the archives. > >> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >> to by IEEE754). > > If they should give exactly the same result always, then the > first thing that needs to be fixed is any standard that allows > them to differ. The language standards are so much looser than > the numerical standard, even now, that I'd rather see C and C++ > catch up to IEEE7541985 before asking them to go beyond. > >> Browsing some gcc lists, I did see some comments that suggest that >> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x). This would >> make sqrt and pow consistent (whether or not technically correct). > > Not so many years ago (perhaps when msvcrt was written), it was > thought that correctlyrounded transcendental functions weren't > feasible in practice, so library authors had a looser attitude > toward exactness. If you have an ideal implementation of sqrt() > but a pow() that's only approximate, should you specialcase > the power of exactly 0.5? If you do, then pow() and sqrt() will > be consistent, but then the powers > x^(0.5*(1epsilon)) > x^(0.5 ) > x^(0.5*(1+epsilon)) > might not be monotone. Some would call that a bad tradeoff. Well, I spent some efforts into log, exp, and pow implementation. Also some other basemath functions I improved in mingww64's math library in two terms. a) Make those functions behave as ISOC99 specification tells in cornercases. and b) Improve accuracy in terms of gcc's internal used gmp. As we noticed here result differences. The following example illustrates the issues pretty well: #include <stdio.h> #include <math.h> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; int main() { double r[3]; long double rl[3]; int i; for (i = 0; i < 3; i++) __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], (r[i] = pow (abc[i][0], abc[i][1]))); __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, 3.0), pow (4.0, 0.5)); r[0] = pow (2.2, 3.1); r[1] = pow (2.0, 3.0); r[2] = pow (4.0, 0.5); __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); for (i = 0; i < 3; i++) __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], (rl[i] = powl (ab[i][0], ab[i][1]))); __mingw_printf ("%.20LF %.20LF %.20LF\n", powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl (4.0L, 0.5L)); rl[0] = powl (2.2L, 3.1L); rl[1] = powl (2.0L, 3.0L); rl[2] = powl (4.0L, 0.5L); __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); return 1; } The interesting issue is here that gcc uses for constant mathcalculations gmplibrary to get result. For more complex variants, it uses crt's math. By this it is important for a gcc based runtime to be in the IEEE 754 floating point compatible to gmp's results. Btw new pow, log, and exp variant are slightly faster then variant in msvcrt, but this is more or less just a sideeffect. I admit that this might be not that obvious to users, why results here are different, but for a gcc based toolchain we need to play nice with gcc's internal assumptions. Regards, Kai 
From: Kai Tietz <ktietz70@go...>  20110430 09:51:40

2011/4/30 Kai Tietz <ktietz70@...>: > 2011/4/30 Greg Chicares <gchicares@...>: >> On 20110429 19:55Z, K. Frank wrote: >>> >>> By the way, could someone confirm that mingw does use msvcrt for >>> sqrt and pow? >> >> http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src >> >> Looks to me like sqrt(double) calls into msvcrt. But pow() is a >> different story. Danny implemented it in libmingwex because of >> some qualityofimplementation problem with msvcrt; then people >> complained that the libmingwex version was slower, and I can't >> remember where it wound up, but it's all in the archives. >> >>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >>> to by IEEE754). >> >> If they should give exactly the same result always, then the >> first thing that needs to be fixed is any standard that allows >> them to differ. The language standards are so much looser than >> the numerical standard, even now, that I'd rather see C and C++ >> catch up to IEEE7541985 before asking them to go beyond. >> >>> Browsing some gcc lists, I did see some comments that suggest that >>> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x). This would >>> make sqrt and pow consistent (whether or not technically correct). >> >> Not so many years ago (perhaps when msvcrt was written), it was >> thought that correctlyrounded transcendental functions weren't >> feasible in practice, so library authors had a looser attitude >> toward exactness. If you have an ideal implementation of sqrt() >> but a pow() that's only approximate, should you specialcase >> the power of exactly 0.5? If you do, then pow() and sqrt() will >> be consistent, but then the powers >> x^(0.5*(1epsilon)) >> x^(0.5 ) >> x^(0.5*(1+epsilon)) >> might not be monotone. Some would call that a bad tradeoff. > > Well, I spent some efforts into log, exp, and pow implementation. Also > some other basemath functions I improved in mingww64's math library > in two terms. a) Make those functions behave as ISOC99 specification > tells in cornercases. and b) Improve accuracy in terms of gcc's > internal used gmp. As we noticed here result differences. > > The following example illustrates the issues pretty well: > > #include <stdio.h> > #include <math.h> > > double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; > long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; > int main() > { > double r[3]; > long double rl[3]; > int i; > for (i = 0; i < 3; i++) > __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], > (r[i] = pow (abc[i][0], abc[i][1]))); > __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, > 3.0), pow (4.0, 0.5)); > r[0] = pow (2.2, 3.1); > r[1] = pow (2.0, 3.0); > r[2] = pow (4.0, 0.5); > __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); > > for (i = 0; i < 3; i++) > __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], > (rl[i] = powl (ab[i][0], ab[i][1]))); > __mingw_printf ("%.20LF %.20LF %.20LF\n", > powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl > (4.0L, 0.5L)); > rl[0] = powl (2.2L, 3.1L); > rl[1] = powl (2.0L, 3.0L); > rl[2] = powl (4.0L, 0.5L); > __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); > return 1; > } > > The interesting issue is here that gcc uses for constant > mathcalculations gmplibrary to get result. For more complex > variants, it uses crt's math. By this it is important for a gcc based > runtime to be in the IEEE 754 floating point compatible to gmp's > results. > > Btw new pow, log, and exp variant are slightly faster then variant in > msvcrt, but this is more or less just a sideeffect. > > I admit that this might be not that obvious to users, why results here > are different, but for a gcc based toolchain we need to play nice with > gcc's internal assumptions. > > Regards, > Kai > 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 specialcase in pow the case for y == 0.5 to solve this. Kai 
From: K. Frank <kfrank29.c@gm...>  20110430 15:15:19

Hello Kai! On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote: > 2011/4/30 Kai Tietz <ktietz70@...>: >> 2011/4/30 Greg Chicares <gchicares@...>: >>> On 20110429 19:55Z, K. Frank wrote: >>>> >>>> By the way, could someone confirm that mingw does use msvcrt for >>>> sqrt and pow? >>> >>> http://cygwin.com/cgibin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src >>> >>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a >>> different story. Danny implemented it in libmingwex because of >>> some qualityofimplementation problem with msvcrt; then people >>> complained that the libmingwex version was slower, and I can't >>> remember where it wound up, but it's all in the archives. >>> >>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required >>>> to by IEEE754). >>> ... >>> Not so many years ago (perhaps when msvcrt was written), it was >>> thought that correctlyrounded transcendental functions weren't >>> feasible in practice, so library authors had a looser attitude >>> toward exactness. If you have an ideal implementation of sqrt() >>> but a pow() that's only approximate, should you specialcase >>> the power of exactly 0.5? If you do, then pow() and sqrt() will >>> be consistent, but then the powers >>> x^(0.5*(1epsilon)) >>> x^(0.5 ) >>> x^(0.5*(1+epsilon)) >>> might not be monotone. Some would call that a bad tradeoff. >> >> Well, I spent some efforts into log, exp, and pow implementation. Also >> some other basemath functions I improved in mingww64's math library >> in two terms. a) Make those functions behave as ISOC99 specification >> tells in cornercases. and b) Improve accuracy in terms of gcc's >> internal used gmp. As we noticed here result differences. >> >> The following example illustrates the issues pretty well: >> >> #include <stdio.h> >> #include <math.h> >> >> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } }; >> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } }; >> int main() >> { >> double r[3]; >> long double rl[3]; >> int i; >> for (i = 0; i < 3; i++) >> __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1], >> (r[i] = pow (abc[i][0], abc[i][1]))); >> __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0, >> 3.0), pow (4.0, 0.5)); >> r[0] = pow (2.2, 3.1); >> r[1] = pow (2.0, 3.0); >> r[2] = pow (4.0, 0.5); >> __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]); >> >> for (i = 0; i < 3; i++) >> __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1], >> (rl[i] = powl (ab[i][0], ab[i][1]))); >> __mingw_printf ("%.20LF %.20LF %.20LF\n", >> powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl >> (4.0L, 0.5L)); >> rl[0] = powl (2.2L, 3.1L); >> rl[1] = powl (2.0L, 3.0L); >> rl[2] = powl (4.0L, 0.5L); >> __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]); >> return 1; >> } Do I understand your point correctly that you are comparing gcc's compiletime and runtime values for pow? For the record, when I run this on a mingw64 build: gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902  mingww64/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 runtime evaluation only for powl (2.2, 3.1), i.e., the longdouble case. >> The interesting issue is here that gcc uses for constant >> mathcalculations gmplibrary to get result. For more complex >> variants, it uses crt's math. By this it is important for a gcc based >> runtime to be in the IEEE 754 floating point compatible to gmp's >> results. Yes, that's a good point. I hadn't considered the issue of compiletime evaluation, and, ideally, it should agree with runtime evaluation, adding yet another constraint to our wish list. Furthermore, if we take the idea of a crosscompiler seriously, we would also like the compiletime evaluations of a crosscompiler and native compiler to agree. Hmm... So, we would like sqrt (x) and pow (x, 0.5) to agree. We would like compiletime and runtime evaluations to agree. We would like crosscompilers 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 agreedupon standard, such as IEEE754). Just to make the reasoning explicit, if I perform my runtime evaluation on platform A, e.g., using platform A's FPU, then a B>A crosscompiler running on platform B won't have platform A's FPU available to perform consistent compiletime evaluations. So you either give up some of the abovewishedfor 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 specialcase 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 tradeoff. (Ideally, both pow and sqrt would be "exact" and therefore agree with one another automatically.) By the way, I also see a problem with runtime powl for powl (2.2L, 0.5L). It disagrees slightly from sqrtl (2.2L) and from the compiletime 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  mingww64/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 Delaney <delaneyrm@ea...>  20110430 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 Tietz <ktietz70@go...>  20110430 21:42:03

2011/4/30 Bob Delaney <delaneyrm@...>: > > 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 80bit IEEE 754 floating point. Which means for our purposes in gcc, it has right accuracy. Kai 
From: K. Frank <kfrank29.c@gm...>  20110430 22:13:15

Hello Bob and Kai! On Sat, Apr 30, 2011 at 5:41 PM, Kai Tietz <ktietz70@...> wrote: > 2011/4/30 Bob Delaney <delaneyrm@...>: >> >> 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 52bit 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 80bit floating point number. This has a 64bit 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 80bit 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 80bit floatingpoint number.) > Kai Thanks to all for the insight. K. Frank 
From: Kai Tietz <ktietz70@go...>  20110430 22:14:14

2011/4/30 Kai Tietz <ktietz70@...>: > 2011/4/30 Bob Delaney <delaneyrm@...>: >> >> 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 80bit 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 (32bit): digits for mantissa:24 digit:6, min exp:125 (37 10th) max exp: 128 (38 10th) Epsilon:1.19209e007 double (64bit): digits for mantissa:53 digit:15, min exp:1021 (307 10th) max exp: 1024 (308 10th) Epsilon:2.22045e016 long double (80bit): digits for mantissa:64 digit:18, min exp:16381 (4931 10th) max exp: 16384 (4932 10th) Epsilon:1.0842e019 Which means that indeed the accuracy of 20 digits after decimal point are pretty precise here. Regards, Kai 
From: Keith Marshall <keithmarshall@us...>  20110501 07:43:46

On 30/04/11 23:14, Kai Tietz wrote: > long double (80bit): > digits for mantissa:64 > digit:18, min exp:16381 (4931 10th) > max exp: 16384 (4932 10th) > Epsilon:1.0842e019 > > 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 80bit 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 Rockett <p.rockett@sh...>  20110501 11:01:13

On 01/05/11 08:43, Keith Marshall wrote: > On 30/04/11 23:14, Kai Tietz wrote: >> long double (80bit): >> digits for mantissa:64 >> digit:18, min exp:16381 (4931 10th) >> max exp: 16384 (4932 10th) >> Epsilon:1.0842e019 >> >> 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 20digit accuracy. They are very different things. > You cannot get 20 decimal digits precision after the decimal point, with > 80bit 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 finiteprecision 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 "illconditioning". P. 