On Aug 11, 6:03 pm, Jonathan Bober <[EMAIL PROTECTED]> wrote: > I have just noticed that using the C type long double from within sage > doesn't work the way that I've expected it to. > > The issue is a little complicated, and other people on this list > probably know more about it than I do, but, briefly, the problem stems > from the fact that x86 fpu registers support extended precision (64 bit > mantissa), while a C double is just double precision (53 bit mantissa). > This causes unexpected results sometimes because floating point > arithmetic will be done on the processor at extended precision, and then > stored in memory at double precision, which means that that value of a > double could theoretically change at any time. This is very bad for > packages like quaddouble. > > It is possible to set the fpu control word so that the fpu only uses > double precision, which is one way to solve this problem. Apparently, > this is done by default on Windows, but not in Linux (and I have no idea > about OS X). One problem with this, however, at least on Linux, is that > long doubles no longer 64 bits of precision. > > It seems that, perhaps accidentally, sage is setting the fpu to use > double precision. ... > There are a few possible ways to fix this, I think. > > (1)One way is to rewrite the quaddouble wrapper. Then the fpu control > word would need to be saved, set, and reset every time arithmetic was > done on a quad double (not a good thing, probably). > > (2)Another possibility is to decide that it is a good idea to set the > fpu to just use double precision, and to make sure that no sage > libraries ever expect extended precision. This might actually be the > only really cross platform solution that isn't lots of work -- I'm not > sure how long doubles will work on Windows or on OS X.
Note that the transcendental functions (trigonometric, etc.) in libm on my Debian Linux box are much more accurate if the fpu is in extended precision mode than if the fpu is in double precision mode. (If you're using RDF floating-point numbers in SAGE, then SAGE will use the gsl_ routines instead of the libm routines. I don't know what the effect of extended precision is on the gsl_ routines. If you're using RealField floating-point numbers, then SAGE uses the mpfr routines, and the fpu is not used at all.) > (3a)What I think might be the best idea, at least on Linux, is to change > the compilation settings for quad double so that the fpu fix is not > needed. There are two ways to do this: If a processor supports sse2, > then passing gcc -march=whatever -msse2 -mfpmath=sse (maybe the -march > isn't needed) will cause gcc to use sse registers and instructions for > doubles, and these have the proper precision. In fact, gcc already does > this by default for x86-64 cpus, so the quad double package doesn't even > need the fpu fix on those architectures. Also, this has the added > benefit of being faster. Yes, this certainly sounds like a good idea. In fact, most of SAGE should be built this way whenever possible. (This might mean having two different binary downloads for x86 Linux, though.) Do people use pre-SSE2 x86 processors to run SAGE? (According to http://en.wikipedia.org/wiki/SSE2, Intel added SSE2 with the Pentium 4 in 2001, and AMD added SSE2 with their 64-bit processors in 2003.) > (3b)For processors that don't support sse2, gcc can be passed the > -ffloat-store option, which fixes the problem by storing doubles in > memory after every operation, ensuring that they are always correctly > rounded to double precision. This slows things down a little bit, but > would probably be much simpler than option (1). Unfortunately this doesn't work. The float-store option rounds all numbers to double precision, avoiding the problem where a number invisibly changes in the middle of a computation; but it does not correctly round to double precision, so it is not usable for quad double. > Personally, I think that I like options 3a and 3b. These would probably > require rewriting the configure script for quad double. I don't know how > to do that, but it probably isn't that hard. > > Option 2 might be better, though, depending on how Windows (I know sage > doesn't work on Windows now, but it would be best not to make it any > harder to port to Windows in the future -- also, I'm not whether this > effects vmware) and OS X work. Hopefully other people on this list know > more about floating point arithmetic issues than I do. I've got my own code that needs to use double precision floating point on x86 in SAGE (my real root isolation routines, which I am still planning to contribute to SAGE soon); so I hope whatever solution we come up with isn't specialized to only work with the quad double package. Carl --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~----------~----~----~----~------~----~------~--~---
