Sign in Join

Live chat for Free Pascal, Lazarus, Delphi and every Pascal programmer.

Mailing lists / fpc-pascal / 2012-07

Floating point discrepancy in x64 mode

3 messages · 2 people ·

Reply by email Subscribe to this list Your subscription

Read-only archive. Use j and k to move between messages.

OBones

Hello all,

Let's consider the following test program:

program test;
function DirectCompute: Double;
var
   k: Integer;
begin
   k := 0;
   Result := (k+0.5)/12;
end;
function IntermediateCompute: Double;
var
   k: Integer;
begin
   k := 0;
   Result := k+0.5;
   Result := Result / 12;
end;
begin
   WriteLn(DirectCompute);
   WriteLn(IntermediateCompute);
end.

When I compile it with ppc386 and run it, I get the following output:

4.16666666666667E-002
4.16666666666667E-002

which is totally expected.
However, when I compile it with ppcrossx64 and run it, I get the following output:

4.16666679084301E-002
4.16666666666667E-002

The error is very significant here as in my real world program, it gets propagated through successive multiplications.
Looking at the bit values, I get this:

$3FA5555560000000 for the "bogus value"
$3FA5555555555555 for the expected value

I tried setting the optimization options ( O-, O1, O2, O3) but to no avail. All this is with FPC 2.6.0 (2011/12/25) as it is the latest stable release.

Many thanks in advance for your help

Jonas Maebe
to OBones
OBones wrote 10 lines
When I compile it with ppc386 and run it, I get the following output: 4.16666666666667E-002 4.16666666666667E-002 which is totally expected. However, when I compile it with ppcrossx64 and run it, I get the following output: 4.16666679084301E-002 4.16666666666667E-002

The reason is that Win64 (and only Win64, not other x86-64 platforms) has deprecated the use of the 80 bit 80x87 fpu. You now have to use the SSE unit for floating point math, which supports up to 64 bit double precision like most FPU's on other platforms.

Now, apart from supporting up to 80 bit precision, there's another peculiarity about the 80x87: it does not have any opcodes to perform operations at less than 80 bits precision. So a division always uses 80 bits precision, and afterwards you can truncate the result back down to whatever you need. The SSE unit on the other hand has separate opcodes for both 32 bit and 64 bit floating point math.

Now, your code has this expression:

 Result := (k+0.5)/12;

with k an integer. 0.5 and 12 can both be represented exactly using single precision, so the compiler evaluates that entire expression using single precision math. However, as mentioned on the 80x87 the addition and division will nevertheless still be evaluated using 80 bits precision due to the 80x87 limitations. On all other platforms, the expression will be evaluated entirely using 32 bit floating point math.

Possible solutions:
a) typecast one of the operands of the division to double (not that this is Delphi-incompatible) b) Compile with -CF64 to force all floating point constants to become at least 64 bit in precision, even if they can be represented exactly using 32 bit

Jonas

OBones
to Jonas Maebe
Jonas Maebe wrote
Possible solutions: a) typecast one of the operands of the division to double (not that this is Delphi-incompatible) b) Compile with -CF64 to force all floating point constants to become at least 64 bit in precision, even if they can be represented exactly using 32 bit

I went with the "-CF64" option and it works just fine or us

Many thanks for the fast answer.