computing exp(1.0) on SPARC.
in [Solaris]
|
From: John D Groenveld on 31 Dec 2009 12:27 In article <a9121688-3815-4fef-be3f-b7be1f633959(a)l30g2000yqb.googlegroups.com>, David Kirkby <drkirkby(a)gmail.com> wrote: >Sun Studio on SPARC is unique among all the different compilers tried >(gcc, HP, Sun Studio) on different platforms (x86, SPARC, PA-RISC and The Sun Studio developers have some web forums: <URL:http://developers.sun.com/sunstudio/community/index.jsp> John groenveld(a)acm.org
From: Chris Ridd on 31 Dec 2009 12:40 On 2009-12-31 15:33:40 +0000, David Kirkby said: > Thank you Stuart. > > But the issue here is not simply different compilers on different > architectures, but different compilers on the *same* architecture. > > Sun Studio on SPARC is unique among all the different compilers tried > (gcc, HP, Sun Studio) on different platforms (x86, SPARC, PA-RISC and > some unknown to me processor on which AIX 6.1 is running). I can only > check with gcc on AIX - there is no valid IBM compiler. > > 1) The correct answer is > 2.71828182845904523536028747135266249775724709369995957496696... goes > on forever. > > 2) Rounded to within the limits specific by DBL_EPSILON in float.h > (2.2204460492503130808473 x 10^-16), E should be 2.7182818284590451 (a > better answer, using no more decimal digits would be > 2.7182818284590452) > > 3) Sun Studio on SPARC gives: 2.7182818284590455, which is not a > correctly rounded representation of E. cc -xlibmopt (Studio 12) on SPARC gives 2.7182818284590451 BTW. The C User's Guide has a number of options that affect floating point, but -xlibmopt was the only one I found that changed the result to end in 51. -- Chris
From: Raymond Toy on 31 Dec 2009 13:02 David Kirkby wrote: > On Dec 31, 12:59 pm, Stuart Biggar <sbig...(a)email.arizona.edu> wrote: >> David, >> >> I'm no expert in floating point but do use it. I wouldn't >> expect different compilers on different architectures to get >> consistent results in double precision to better than what >> is defined in float.h (from a SPARC running Solaris 10U4): >> >> DBL_DIG is defined as 15 and >> DBL_EPSILON is defined as a bit over 2.2e-16 >> >> So I wouldn't expect computations in double precision to >> have more than 15 significant digits. Some (few) values may >> be exact but most will not be. >> >> If you need better precision, try 80-bit on x86 or quad >> precision on SPARC (LDBL). That should result in 18 or 33 >> significant digits. In any case, floating point is only >> a (probably inaccurate to some degree) representation of >> the real number you are trying to calculate. There are >> all sorts of traps for the unwary in floating point. I've >> found that you have to be very careful how you compute >> things - mathematically equal ways of representing a >> quantity can compute quite differently with precision >> VERY much worse than you might expect. Extended precision >> can help some but isn't the real fix (and it may slow >> things down considerably). >> >> Stuart > > Thank you Stuart. > > But the issue here is not simply different compilers on different > architectures, but different compilers on the *same* architecture. > > Sun Studio on SPARC is unique among all the different compilers tried > (gcc, HP, Sun Studio) on different platforms (x86, SPARC, PA-RISC and > some unknown to me processor on which AIX 6.1 is running). I can only > check with gcc on AIX - there is no valid IBM compiler. > > 1) The correct answer is > 2.71828182845904523536028747135266249775724709369995957496696... goes > on forever. > > 2) Rounded to within the limits specific by DBL_EPSILON in float.h > (2.2204460492503130808473 x 10^-16), E should be 2.7182818284590451 (a > better answer, using no more decimal digits would be > 2.7182818284590452) > > 3) Sun Studio on SPARC gives: 2.7182818284590455, which is not a > correctly rounded representation of E. > > 4) Every other compiler / OS gives the the correctly rounded number. > > * GCC on SPARC gives 2.7182818284590451 > * Sun Studio on Solaris x86 gives 2.7182818284590451 > * GCC on Solaris x86 gives 2.7182818284590451 > * GCC on Linux x86 gives 2.7182818284590451 > * HP's compiler on HP-UX gives 2.7182818284590451 > * GCC on HP-UX gives 2.7182818284590451 > * GCC on AIX gives 2.7182818284590451 Why are these correctly rounded when the true result is 2.718...459045235? Why should they not end in 452, not 451? Perhaps gcc on sparc uses libm but sunstudio provides its own, different libm? The source code for Sun's libm is available and can be found in glibc. That version might be old though, but it could tell you why the answer is slightly different. And the printed output depends on libc so that can also affect the output. Glibc used to have some issues what that but perhaps that has been fixed too. Ray
From: David Kirkby on 31 Dec 2009 20:42 On Dec 31 2009, 6:02 pm, Raymond Toy <toy.raym...(a)gmail.com> wrote: > David Kirkby wrote: > > * GCC on SPARC gives 2.7182818284590451 > > * Sun Studio on Solaris x86 gives 2.7182818284590451 > > * GCC on Solaris x86 gives 2.7182818284590451 > > * GCC on Linux x86 gives 2.7182818284590451 > > * HP's compiler on HP-UX gives 2.7182818284590451 > > * GCC on HP-UX gives 2.7182818284590451 > > * GCC on AIX gives 2.7182818284590451 > > Why are these correctly rounded when the true result is > 2.718...459045235? Why should they not end in 452, not 451? Hi Ray, It's impossible to 2.7182818284590452 in the IEEE format. All double precision real numbers in IEEE 754 format are of the form: (S)^-1 * C * 2^q S=Sign bit C=Significand q=Exponent The following two numbers 6121026514868073 * 2^(-51) = 2.71828182845904509079559829842764884233474731445312... and 6121026514868074 * 2^(-51) = 2.71828182845904523536028747135266249775724709369995957 ... can be represented, but nothing in between them. The closest of them to E is 2.71828182845904509079559829842764884233474731445312. > Perhaps gcc on sparc uses libm but sunstudio provides its own, different > libm? It would appear from some work done by someone else, that the Sun maths library computes this to be 2.7182818284590455.... In the trivial program I wrote, gcc inlines E which gives a more accurate result. But in a more complex program, gcc will call the Sun maths library, which will return 2.7182818284590455. The program does not have to be very complex. This example, written by Gonzalo Tornaria (see the sage-devel mailing list for more), shows how gcc will also give 2.7182818284590455 on SPARC. #include <math.h> #include <stdio.h> #include <stdlib.h> int main(int argc, char **argv) { double x = 1.0; if (argc>1) x = atof(argv[1]); printf("%.16lf\n",exp(x)); } From what I gather, this does not necessarily mean the maths library is broken. The libraries do not guarantee to return the most accurate result, so there will be some degree of tradeoff with accuracy and speed. I'm guessing that since the SPARC processor uses 64-bit for the computations, but the Intel processors uses 80 bits internally, there is justification for choosing a slightly less accurate but faster implementation on SPARC. I'm no expert on this, but have a better understanding of it then I did a day or two ago. Dave
From: Stuart Biggar on 1 Jan 2010 09:06
On 12/31/09 6:42 PM, David Kirkby wrote: > I'm guessing that since the SPARC processor uses 64-bit for the > computations, but the Intel processors uses 80 bits internally, there > is justification for choosing a slightly less accurate but faster > implementation on SPARC. > > I'm no expert on this, but have a better understanding of it then I > did a day or two ago. > > Dave Dave, On various Intel platforms, double precision may use SSE where the math is done at 64-bit, not the extended precision of 80-bit used in the 80N87 (math co-processor). Also at least one Windows C compiler silently converts long doubles to doubles to allow use of the faster SSE stuff - gives "interesting" results when the code specifically tries to use 80-bits for accuracy. If you need an accurate value for e, it is defined in one of the include files along with various other constants like pi and fractions of pi (from math.h on SPARC S10U4): #define M_E 2.7182818284590452354 #define M_LOG2E 1.4426950408889634074 #define M_LOG10E 0.43429448190325182765 #define M_LN2 0.69314718055994530942 #define M_LN10 2.30258509299404568402 #define M_PI 3.14159265358979323846 #define M_PI_2 1.57079632679489661923 #define M_PI_4 0.78539816339744830962 #define M_1_PI 0.31830988618379067154 #define M_2_PI 0.63661977236758134308 #define M_2_SQRTPI 1.12837916709551257390 #define M_SQRT2 1.41421356237309504880 #define M_SQRT1_2 0.70710678118654752440 With Studio you can also use optimized math libraries. I have no idea if the values of things would change by linking with them (typically done automatically if you optimize using -fast and maybe other options or manually during the link step). Stuart |