On Apr 29, 2007, at 1:01 PM, Tim Prince wrote:
[EMAIL PROTECTED] 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
