Re: Effects of newly introduced -mpcX 80387 precision flag

[Bradley Lucier <lucier at math dot purdue dot edu>]

On Apr 29, 2007, at 1:01 PM, Tim Prince wrote:

> lucier@math.purdue.edu wrote:
>
> > I just (re-)discovered these tables giving maximum known errors in  
> > some libm functions when extended precision is enabled:
> > http://people.inf.ethz.ch/gonnet/FPAccuracy/linux/summary.html
> > and when the precision of the mantissa is set to 53 bits (double  
> > precision):
> > http://people.inf.ethz.ch/gonnet/FPAccuracy/linux64/summary.html
> > This is from 2002, and indeed, some of the errors in double- 
> > precision results are hundreds or thousands of times bigger when  
> > the precision is set to 53 bits.
> This isn't very helpful.  I can't find an indication of whose libm  
> is being tested, it appears to be an unspecified non-standard  
> version of gcc, and a lot of digging would be needed to find out  
> what the tests are.

The C code for the tests is in a subdirectory of:

http://people.inf.ethz.ch/gonnet/FPAccuracy/

There is also a file "all.tar.Z" there, and a few html files with  
commentary.

> It makes no sense at all for sqrt() to break down with change in  
> precision mode.

If you do an extended-precision (80-bit) sqrt and then round the  
result again to a double (64-bit) then those two roundings will  
increase the error, sometimes to > 1/2 ulp.

To give current results on a machine I have access to, I ran the  
tests there on

vendor_id       : AuthenticAMD
cpu family      : 15
model           : 33
model name      : Dual Core AMD Opteron(tm) Processor 875

using

euler-59% gcc -v
Using built-in specs.
Target: x86_64-unknown-linux-gnu
Configured with: ../configure --prefix=/pkgs/gcc-4.1.2
Thread model: posix
gcc version 4.1.2

on an up-to-date RHEL 4.0 server (so whatever libm is offered there),  
and, indeed, the only differences that it found were in 1/x, sqrt(x),  
and Pi*x because of double rounding.  In other words, the code that  
went through libm gave identical answers whether running on sse, x87  
(extended precision), or x87 (double precision).

I don't know whether there are still math libraries for which  
Gonnet's 2002 results prevail.

Brad

The only change I made to the code was

euler-63% rcsdiff header.h
===================================================================
RCS file: RCS/header.h,v
retrieving revision 1.1
diff -r1.1 header.h
66c66,68
< double scalb( double x, double n );
---
> /* double scalb( double x, double n ); */
> #include <math.h>
> #include <stdio.h>

and the scripts I used to run the code were

euler-60% cat generate-results
#! /bin/tcsh
set files = `ls *.c | sed 's/\.c//'`
foreach file ( $files )
  gcc -O3 -mfpmath=sse -o $file $file.c -lm
  ./$file >! $file.-mfpmath=sse
  gcc -O3 -mfpmath=387 -o $file $file.c -lm
  ./$file >! $file.-mfpmath=387
  gcc -O3 -mfpmath=387 -DFORCE_FPU_DOUBLE -o $file $file.c -lm
  ./$file >! $file.-mfpmath=387-DFORCE_FPU_DOUBLE
  rm $file
end
euler-61% cat compare
#! /bin/tcsh
set files = `ls *.c | sed 's/\.c//'`
foreach file ( $files )
    echo comparing $file.-mfpmath=387 $file.-mfpmath=387- 
DFORCE_FPU_DOUBLE
    diff $file.-mfpmath=387 $file.-mfpmath=387-DFORCE_FPU_DOUBLE
    echo comparing $file.-mfpmath=387 $file.-mfpmath=sse
    diff $file.-mfpmath=387 $file.-mfpmath=sse
end