Hi Fabian, This was a bit bizarre particularly as when I tried it on Windows there didn't seem to be a problem. It is a problem with the precision and when rounding occurs. Poly/ML uses the old x87 floating point unit on both X86/32 and X86/64. As far as I'm aware all x86/64 machines support the SSE2 instructions so other compilers use SSE2 instructions on 64-bit but, by default use the x87 instructions on 32-bit.
The x87 instructions use 80-bit registers internally and only round to 64-bits when the registers are stored to double-precision values. The SSE2 instructions are 64-bit so all the results are the same whether values are stored to memory or not. The x87 has bits in the control word that can allow it to compute intermediate results to 80-bit, 64-bit or 32-bit precision. The question is whether this should be set to 80-bit or 64-bits?
Currently Poly/ML sets it to 80-bit precision at the start. As far as I've been able to ascertain that is what other x86/32 compilers do. The idea is, I suppose, to compute intermediate results to the best precision available. It does, though, produce the strange results that you've encountered. I think the reason for the difference in results is whether the function can be computed at compile-time or not.
I'm coming round to the view that it should instead set the FPU to 64-bits so that intermediate results are the same whether they remain in the register or are written out to memory and reloaded. It looks as though the libraries in Visual C set that anyway which is why I only saw the problem when I tested on Linux. This would also mean that results would not change when, as I plan, I change to using the SSE2 instructions on x86/64.
Best regards, David
On 07/10/2015 14:07, Fabian Immler wrote:
Hi David,
Perhaps I should rephrase my enquiry to: "Do I get any guarantees from the implementation of type real in PolyML"? Naively, I would have expected, that every subexpression is evaluated according to the IEEE-754 standard (in double precision).
I skimmed the implementation in libpolyml/reals.cpp, it looks like all intermediate results are stored in and retrieved from memory (real_arg, real_result). Is that sufficient to comply to IEEE-754? I.e., can there not occur "double rounding" (as in [1, section 3.3.1])? Or isn't there the need to use something like FLT_EVAL_METHOD=1 [2,3]?
[1] http://www.springer.com/us/book/9780817647049 [2] https://en.wikipedia.org/wiki/C99#IEEE.C2.A0754_floating_point_support [3] http://www.cplusplus.com/reference/cfloat/
Best regards, Fabian
Am 07.10.2015 um 12:26 schrieb Fabian Immler <immler at in.tum.de>:
Hi David,
I noticed some strange behavior when working with type real in PolyML: In the following example, where I have "testres = twosum a", depending on whether I call "testres" or "twosum a" I get different results.
If I take the definitions out of the structure Test, the results are the same. They are also the same if e.g. print-statements are inserted in between the assignments to b, s, sb and r. The results are also the same when evaluating inside Isabelle with the debugger turned on.
Do you have an idea what is going on there? Is it possible that some intermediate results are kept in higher precision in the FPU? Does the compiler perform some "optimizations" on expressions of type real?
Best regards, Fabian
structure Test = struct
fun one n = (if n = 0 then 100.0 else one (n - 1));
val a : real = one 1;
fun twosum a = let val b = 1.0/1243313.0; val s = a + b val sb = s - b val r = s - sb in r end;
val testres = twosum a; end;
open Test; val a = twosum a val b = testres