On 4/11/2018 7:03 PM, Philip McGrath wrote:
From one following along who knows fairly little about floating-point
math (the Toronto/McCarthy paper looks very informative!):
On Tue, Apr 10, 2018 at 12:13 AM, George Neuner <gneun...@comcast.net
<mailto:gneun...@comcast.net>> wrote:
As Philip McGrath mentioned already, you can specify single
precision at least for input values ... but I think inexact (i.e.
floating point) math in Racket all is done at maximum hardware
precision and the result coerced back to the input precision if
possible. For purposes of numerical analysis, that is not the
same as performing all operations at the input precision.
Is this confirmed to be what Racket does?
Yes and no. I did not intend to attribute the issue to Racket - it's
really a hardware issue.
Racket follows the modern convention of using the CPU's SIMD hardware to
do FP math. SIMD hardware (typically) does not widen results, so you
get back the same precision you put in ... but internally there was an
intermediate double-wide mantissa computed that maybe didn't fit into
the output.
You see:
single X single -> single
double X double -> double
The hardware does:
single X single -> double -> single
double X double -> quad -> double
If the result overflows or underflows, you get INF or zero rather than a
correct but wider answer.
If you are trying to check your FPU using software emulations: i.e.
bigint-rationals or bigfloats, then you need to know exactly how the
computations are being done. More on that below.
Modern x86 have both SIMD and x87 co-processor. The SIMD unit was
introduced on the Pentium III with single precision and was expanded to
double precision on the Pentium IV and later. The modern convention to
use the SIMD unit exclusively leaves the x87 idle all the time. There's
a good reason for it - the SIMD unit is faster because it is more
tightly integrated ... the x87 originally was a separate chip and even
though it now is integrated, there still is a performance penalty for
having to use the co-processor interface protocol. And the SIMD unit
produces more consistent results [see below].
The x87 and the SIMD unit do not work the same way, and they may produce
different results from the same calculation. The SIMD unit offers
single or double precision. The x87 offers a 3rd choice: extended
precision - an 80 bit format that has 64 bits of mantissa [vs 53 bits in
an IEEE double]. Extended precision is the default mode of the x87.
The x87 always computes a full (64-bit mantissa) extended precision
result before rounding AT the selected precision and then truncating TO
the selected precision [which clears any trailing bits in the
register]. In the SIMD unit, precision is specified by the
instruction: the registers are (virtually) subdivided into slots which
are sized appropriately for the specified precision [no trailing bits],
and all slots are computed in parallel TO the specified precision.
Very different. And potentially confusing if you are expecting to get
bit-identical results.
GPU or DSP calculations can be even more challenging to check with a CPU
because their FPU hardware often is very different - not just not IEEE
compliant, but sometimes not even close.
Aside for Dale:
If you're on x86, and you find that rationals or bigfloats are too slow,
the x87 may offer a faster way to your more precise calculations.
see:
https://docs.racket-lang.org/reference/extflonums.html?q=math%2Fextflonum#%28tech._extflonum%29
If so, are the observable differences vs. the operations working at
the input precision because of error and error propagation
Yes. Consider what happens at the limits of the selected precision.
Switching to C for a moment, consider:
float X = FLT_MAX;
float Y = FLT_MAX - 4;
float Z = (X + Y) / 2;
float W = (float) (((double) X + Y) / 2);
Z gets compiled something like:
float t1 = X + Y
float t2 = t1 / 2
Z = t2
W get compiled something like:
double t1 = (double) X
double t2 = (double) Y
double t3 = t1 + t2
double t4 = t3 / 2
W = (float) t4
In the case of Z the temporaries are single precision, so (X + Y)
overflows to INF and that carries forward through the rest of the
calculation. In the case of W, the temporaries are double precision, so
the calculation of (X + Y) succeeds and W gets the right answer.
Now you don't see these temporaries [except maybe by compiling in debug
mode and specifying "consistent" floating point], but they exist in the
FPU registers. The actual answers you would get depend on the FPU
hardware and how the code was compiled.
E.g., on x87, keeping the temporaries in registers, both Z and W would
get the same (right) answer if the FPU were in double or extended
precision mode even though the compiler thinks the temporaries are
single precision values. But if the code were compiled for
"consistent" floating point - which forces temporaries to be truncated
to *storage* precision as if saved to memory and reloaded - then after
t1 computed some actual double precision value, that value would be
truncated to single precision regardless of FPU mode. The result would
be wrong, but it would not be INF as it would be with the SIMD unit or
if the x87 had been in single precision mode from the beginning.
Computing Z on the SIMD unit, t1 will always overflow to INF because the
internal double precision intermediate can't be stored back into the
single precision register. So the SIMD unit produces more consistent
results in keeping with programmer expectations. [another reason it is
preferred now]
(i.e., if starting from exact values, the implementation using
double-precision under the hood would sometimes produce results with
less error), or something else?
Double precision calculation MAY have less error than single precision
calculation ... it depends on whether temporaries over or under flow.
-Philip
Hope this helps,
George
--
You received this message because you are subscribed to the Google Groups "Racket
Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to racket-users+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
