[racket-dev] Floating-Point Compliance Testing

Back on list.

A lot of things point to general sloppiness in either the FPU or C 
libraries, but I'd like more information just in case. Can you reply 
with the values of the following expressions on the Athlon?

   (flexpt -1001.0 -1.3407807929942596e+154)
   (flexpt -1001.0 1.3407807929942596e+154)
   (flexpt -0.1 -1.3407807929942596e+154)
   (flexpt -0.1 1.3407807929942596e+154)
   (flexpt -744.4400719213812 -1.3407807929942596e+154)
   (flexpt -744.4400719213812 1.3407807929942596e+154)
   (flexpt -1.0 -1.3407807929942596e+154)
   (flexpt -1.0 1.3407807929942596e+154)

You should get these values:

   0.0 +inf.0 +inf.0 0.0 0.0 +inf.0 1.0 1.0

I think these are cases Racket handles specially instead of handing off 
to C's `pow' on platforms where we know `pow' handles them wrongly. We 
might need to ask Matthew to expand that set of platforms.

Small rounding errors like this:

   ((fl*/error -6.87181640380727e-156 2.3341566035927725e-153)
    0.7079979032237692)

which are only 1 bit off, are probably the cause of these errors:

   ((fl2log 1.5124390004859715e-308 0.0) 4294967296.220986)
   ((fl2log1p 3.799205740343036e+246 1.4492752004468653e+230)
    549755813887.9473))

`fl2log' and `fl2log1p' are 103-ish-bit logarithm implementations used 
in certain tricky subdomains of a few special functions. They assume 
arithmetic is always correct and are very sensitive to rounding errors. 
You're getting about 65-bit output precision for certain inputs. We can 
almost certainly blame the FPU because IEEE 754 requires arithmetic to 
be implemented and correctly rounded.

IEEE 754 only *recommends* typical irrational functions. When they're 
not implemented in hardware, C libraries compute them in software. So I 
don't know whether these and others are the FPU's or C library's fault:

   ((flsin 2.5489254492488616e+52) 22637349860729424.0)
   ((flcos 3.91520018229574e+49) 6369061509154398.0)
   ((fltan 1.6614020610450763e+21) 9158003261155244.0)
   ((flexp 16.938769136983012) 7.0)
   ((flexp 282.52374429474406) 102.0)
   ((flexp -10.0) 4.0)
   ((flexp -708.3964185322641) 124.0)

These errors come from not doing argument reductions carefully enough. 
The trigonometric functions probably don't compute x % 2*pi using a 
high-precision 2*pi when x is large. They seem to be correct enough when 
x is small, though.

The exponential function wasn't tested well at the boundaries of its domain:

   ((flexp 709.782712893384) +inf.0)

709.782712893384 is (fllog +max.0), and (flexp (fllog +max.0)) should be 
near +max.0, not (as I suspect) +inf.0.

I don't know whether Racket should do anything about these errors. I 
don't think it would be too hard to do something like Java's StrictMath, 
but it would take time.

In the meantime, don't use the Athlon for serious numerical computation.

Neil ⊥

On 02/08/2013 12:13 AM, Laurent wrote:
> Hi Neil,
>
> Interaction in a terminal is attached, using Racket 5.3.2.3.
>
> Some details about my machine:
> Linux 3.2.0-37-generic-pae #58-Ubuntu SMP Thu Jan 24 15:51:02 UTC 2013
> i686 athlon i386 GNU/Linux
>
> In particular, I use a 32bits Ubuntu 12.04.2 on a 686 processor, if
> that's of any interest.
>
> Cheers,
> Laurent
>
>
>
>
> On Fri, Feb 8, 2013 at 1:15 AM, Neil Toronto <neil.toronto at gmail.com
> <mailto:neil.toronto at gmail.com>> wrote:
>
>     On 02/07/2013 12:09 PM, Laurent wrote:
>
>         On Thu, Feb 7, 2013 at 5:50 PM, Neil Toronto
>         <neil.toronto at gmail.com <mailto:neil.toronto at gmail.com>
>         <mailto:neil.toronto at gmail.com
>         <mailto:neil.toronto at gmail.com>__>> wrote:
>
>              Today is not that day, but thanks for asking about this
>         anyway. :)
>
>
>         On one machine with Ubuntu 12.10, I get no error, but on another
>         machine
>         with Ubuntu 12.04, I get more than 14000 errors, many of them being
>         +inf.0 and other numbers with big exponents (is my machine
>         really that
>         bad?).
>         Is this exactly the kind of reply you want to avoid for now or
>         are you
>         interested in a report?
>
>
>     Alrighty, you've piqued my interest. Better send it off-list, though. :)
>
>     Neil ⊥