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.
sage/rings/real_rqdf.pyx contains the following code:
cdef class RealQuadDoubleField_class(Field):
"""
Real Quad Double Field
"""
def __init__(self):
fpu_fix_start(self.cwf)
def __dealloc__(self):
fpu_fix_end(self.cwf)
[etc]
__dealloc__() is never called until sage exits, however, since a global
instance of this class is created at startup. This means that all
computations with long doubles in the libraries that sage uses will not
have the expected precision, at least on Linux. (I noticed this while
working on the partition counting code -- the current version in sage
probably won't produce wrong answers because it already is using too
much precision, but this causes problems on the new version that I have
been working on.)
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.
(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.
(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).
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.
Anyway, I mostly just wanted to point out the issue, since other people
likely know more about it than I do. I don't know how this email got to
be so long.
-Jon
--~--~---------~--~----~------------~-------~--~----~
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/
-~----------~----~----~----~------~----~------~--~---
