RE: [boost] Re: Maths Constants

From: Paul A. Bristow (pbristow_at_[hidden])
Date: 2001-03-26 17:39:26

Sorry for my delay in continuing this work on math constants.

You are correct.  Victor Shoup gives more details in his NTL package,
where the 'solution' is to switch off use of the extra bits,
if you value portability more than accuracy.
I propose only to document this.

The latest version of NTL is available at http://www.shoup.net.

An extract follows:

THE INTEL x86 PROBLEM

Although just about every modern processor implements the IEEE
floating point standard, there is still can be problems
on processors that support IEEE extended double precision.
The only processor I know of that supports this is the x86/Pentium.

While extended double precision may sound like a nice thing,
it is not.  Normal double precision has 53 bits of precision.
Extended has 64.  On x86s, the FP registers have 53 or 64 bits
of precision---this can be set at run-time by modifying
the cpu "control word" (something that can be done
only in assembly code).
However, doubles stored in memory always have only 53 bits.
Compilers may move values between memory and registers
whenever they want, which can effectively change the value
of a floating point number even though at the C/C++ level,
nothing has happened that should have changed the value.
Is that sick, or what?

This is a real headache, and if one is not just a bit careful,
the quad_float code will break.  This breaking is not at all subtle,
and the program QuadTest will catch the problem if it exists.

You should not need to worry about any of this, because NTL automatically
detects and works around these problems as best it can, as described below.
It shouldn't make a mistake, but if it does, you will
catch it in the QuadTest program.
If things don't work quite right, you might try
setting NTL_FIX_X86 or NTL_NO_FIX_X86 flags in ntl_config.h,
but this should not be necessary.

Here are the details about how NTL fixes the problem.

The first and best way is to have the default setting of the control word
be 53 bits.  However, you are at the mercy of your platform
(compiler, OS, run-time libraries).  Windows does this,
and so the problem simply does not arise here, and NTL neither
detects nor fixes the problem.  Linux, however, does not do this,
which really sucks.  Can we talk these Linux people into changing this?

The second way to fix the problem is by having NTL
fiddle with control word itself.  If you compile NTL using a GNU compiler
on an x86, this should happen automatically.
On the one hand, this is not a general, portable solution,
since it will only work if you use a GNU compiler, or at least one that
supports GNU 'asm' syntax.
On the other hand, almost everybody who compiles C++ on x86/Linux
platforms uses GNU compilers (although are some commercial
compilers out there that I don't know too much about).

The third way to fix the problem is to 'force' all intermediate
floating point results into memory.  This is not an 'ideal' fix,
since it is not fully equivalent to 53-bit precision (because of
double rounding), but it works (although to be honest, I've never seen
a full proof of correctness in this case).
NTL's quad_float code does this by storing intermediate results
in local variables declared to be 'volatile'.
This is the solution to the problem that NTL uses if it detects
the problem and can't fix it using the GNU 'asm' hack mentioned above.
This solution should work on any platform that faithfully
implements 'volatile' according to the ANSI C standard.

> -----Original Message-----
> From: k.hagan_at_[hidden] [mailto:k.hagan_at_[hidden]]
> Sent: Tuesday, February 13, 2001 9:52 AM
> To: boost_at_[hidden]
> Subject: [boost] Re: Maths Constants
> 
> 
> <pbristow_at_[hidden]> wrote:
> >
> > So using builtin constants might provide a more accurate result,
> > but one that is not as portable (that may be different from
> > other processors). So builtin pi, e ... is better, but maybe
> > also badder?
> 
> I think differences in accuracy are inevitable, even if processors
> claim to be IEEE. On x86 and 68k, the FPU has 80 bit registers.
> Even if you set the operating precision to 64 or 32, it still isn't
> entirely equivalent to an FPU that has true 32 or 64-bit wide
> registers, such as (amusingly enough) the SSE registers added in
> the Pentium 3 and 4, or the registers in a PowerPC or Alpha.
> Intel's 64-bit chip has FP registers even wider still. They are
> deliberately a couple of bits wider than *any* memory based
> floating point type.
> 
> In all these cases, the accuracy that you get out depends on how
> the compiler balances memory accesses against register pressure.
> The standard allows FP calculations to be carried out at "inflated"
> precision anyway, for these kinds of reason, so float calculations
> may be done at double precision as specified in "classic C".
> 
> In such circumstances, I don't think we have the luxury of calling
> the extra accuracy "badder", so we'll just have to call it better.
> 
> (The main reason why the x86 provides these "load constant"
> instructions is to support rounding modes correctly. The most
> accurate representation differs if the rounding mode is towards
> +/- infinity. In one case we want the closest value below the
> constant, and in the other case we want the closest value above it.
> No compile-time constant can reproduce this behaviour, and working
> it out at run-time is not terribly fast, at least for the "long
> double" case! We probably don't care about this. Jens might, in the
> context of his boost interval library, where he needs a pi for
> argument reduction, if I remember correctly.)
> 
> > This is one reason for providing the constants in crude pre-
> > processor #define form for all possible uses, and then providing
> > a C++ representation which should be as portable as possible for
> > the Standard.
> 
> I still don't like the namespace pollution that #define causes.
> What was the objection to something like...?
> 
>    template<class T> struct math_constants /*==namespace*/
>    {
>       static T e() { return T(2.718...L); }
>    };
> 
> 
> 
> 
> 
> 
>