[SOLVED] Rounding issues (only with 64-bit, not with 32-bit)

[Hartmut]

Because system.round() uses "bankers rounding", long time ago I created a "normal" rounding function:

Code: Pascal  [Select][+][-]

1. function round_normal(x: double): int64;
2.    var k: double;
3.    begin
4.    if x < 0 then k:=-0.5 else k:=0.5;
5.    exit(round(int(x+k)));
6.    end;

This worked all the time with 32-bit. Now with 64-bit I found strange rounding issues for larger values around 4.5E+15:
 - only for these larger values and
 - only for odd values and
 - only with 64-bit
then the result of my function is too big by 1.
And function math.SimpleRoundTo() - which is based on the same formula - has exactly the same behaviour.

Examples:
 - for x=4503,599627,370495 both my round_normal() and math.SimpleRoundTo() return the same as 'x' (which is correct)
 - same for x=4503,599627,370496
 - but for x=4503,599627,370497 both my round_normal() and math.SimpleRoundTo() return 4503,599627,370498, which is too big by 1.

I created a demo (attached as compilable project) to test this:

Code: Pascal  [Select][+][-]

1. {$mode objfpc} {$H+}
2. {$OPTIMIZATION OFF} {same results with Optimitation Level 1}
3.  
4. uses math;
5.  
6. function komma6Str(s: string): string;
7.    {inserts Kommas ',' for each 6 digits; Usable for integers, floats, empty strings and negative numbers}
8.    var i,m: integer;
9.    begin
10.    i:=pos('.',s); if i > 0 then dec(i) else i:=length(s);
11.    m:=6; if s[1]='-' then m:=7;  {auch ok wenn s=''}
12.    while i > m do  begin insert(',', s,i-5); dec(i,6); end;
13.    exit(s);
14.    end;
15.  
16. function i64Str(i: int64): string;
17.    {returns 'i' as a string with Kommas ',' for each 6 digits}
18.    var s: string;
19.    begin
20.    str(i,s);
21.    exit(komma6Str(s)); {inserts Kommas ',' for each 6 digits}
22.    end;
23.  
24. function delTrailingDotZeroes(s: shortstring): shortstring;
25.    {cuts all trailing zeros and decimal dot '.' if last char}
26.    begin
27.    while (length(s) > 0) and (s[length(s)] = '0') do  dec(s[0]);
28.    if (length(s) > 0) and (s[length(s)] = '.') then dec(s[0]);
29.    if s='' then s:='0';
30.    exit(s);
31.    end;
32.  
33. function floatStr(x: extended): string;
34.    {returns 'x' as a string with Kommas ',' for each 6 digits}
35.    var s: string;
36.    begin
37.    str(x:0:30, s);
38.    if pos('E',s) > 0 then exit(s); {if exponential format}
39.  
40.    s:=delTrailingDotZeroes(s); {cuts trailing zeros and decimal dot}
41.    exit(komma6Str(s));         {inserts Kommas ',' for each 6 digits}
42.    end;
43.  
44. procedure show_double_interna(d: double);
45.    {shows internal infos about double 'd'}
46.    const K = extended(4503599627370496.0);
47.    var R: TDoubleRec;
48.        ps: string;
49.        mu,ef,v,e2,p: extended;
50.        m: qword;
51.        e: integer;
52.    begin
53.    R:=TDoubleRec(d);    {access a 'double' as a record}
54.    m:=R.Mantissa(true); {get Mantissa including Hidden Bit}
55.    mu:=m / K;           {converted Mantissa value}
56.    e:=R.Exponent;       {decimal Exponent [-1022..+1023] related to 'mu'}
57.    ef:=IntPower(2,e);   {factor resulting from Exponent}
58.    if R.Sign then v:=-1 else v:=+1; {Sign}
59.    e2:=mu * ef * v;     {must be again the same as 'd'}
60.    p:=ef / K;           {current possible precision (resolution) for the range of 'd'}
61.    if p < 0.000001 then str(p,ps) {use exponential format}
62.       else ps:=floatStr(p);       {use normal format}
63.  
64.    write('Mant=$', hexStr(m,14), {Mantissa including Hidden Bit in Hex}
65.          ' Exp=2^', e,           {decimal Exponent [-1022..+1023]}
66.          ' prec=', ps);          {current possible precision for the range of 'd'}
67. // write(' => ', floatStr(e2));  {must be again the same as 'd'}
68.    end;
69.  
70. procedure test_round(i: int64);
71.    {tests rounding issue for value 'i'}
72.    var x,k,h: double;
73.        i1,i2,i3: int64;
74.    begin
75.    x:=i; {make double}
76.    if x < 0 then k:=-0.5 else k:=0.5;
77.  
78.    i1:=round(int(x+k));                {rounding method #1}
79.    i2:=round(math.SimpleRoundTo(x,0)); {rounding method #2}
80.    h:=x+k; i3:=round(int(h));          {rounding method #3 with stored interim result}
81.  
82.    write('i=', i64Str(i), ' => i1=', i64Str(i1), ' i2=', i64Str(i2), ' i3=', i64Str(i3));
83.    if i1 <> i then write(' i1=BAD');
84.    if i2 <> i then write(' i2=BAD');
85.    if i3 <> i then write(' i3=BAD');
86.    writeln;
87.    write('':3); show_double_interna(x); writeln;
88.    end;
89.  
90. procedure Test_rounding_issue;
91.    {tests rounding issues for a range of values.
92.     'MaxDoublePrec1' is the highest 'double' value up to which you can store
93.     EACH integer value in a 'double' variable; For higher values only each 2nd
94.     integer value can be stored and so on; For that this loop must be stopped
95.     there}
96.    const MaxDoublePrec1 = 9007199254740992; {2^53 = 9E15}
97.    var i,min,max: int64;
98.    begin
99.    writeln('FPC-Version ', {$I %FPCVERSION%}, ', ',
100.    {$IFDEF CPU32} '32-bit' {$ELSE} '64-bit' {$ENDIF} );
101.  
102.    i:=round(IntPower(2,52)); {4503,599627,370496}
103. // i:=round(IntPower(2,53)); {9007,199254,740992}
104.    min:=i-2;
105.    max:=i+4;
106.  
107.    i:=min;
108.    repeat test_round(i);
109.           inc(i);
110.    until  (i > max) or (i > MaxDoublePrec1);
111.    end; {Test_rounding_issue}
112.  
113. begin {main}
114. Test_rounding_issue;
115. end.      

The output is:
FPC-Version 3.2.2, 64-bit
i=4503,599627,370494 => i1=4503,599627,370494 i2=4503,599627,370494 i3=4503,599627,370494
   Mant=$1FFFFFFFFFFFFC Exp=2^51 prec=0.5
i=4503,599627,370495 => i1=4503,599627,370495 i2=4503,599627,370495 i3=4503,599627,370495
   Mant=$1FFFFFFFFFFFFE Exp=2^51 prec=0.5
i=4503,599627,370496 => i1=4503,599627,370496 i2=4503,599627,370496 i3=4503,599627,370496
   Mant=$10000000000000 Exp=2^52 prec=1
i=4503,599627,370497 => i1=4503,599627,370498 i2=4503,599627,370498 i3=4503,599627,370498 i1=BAD i2=BAD i3=BAD
   Mant=$10000000000001 Exp=2^52 prec=1
i=4503,599627,370498 => i1=4503,599627,370498 i2=4503,599627,370498 i3=4503,599627,370498
   Mant=$10000000000002 Exp=2^52 prec=1
i=4503,599627,370499 => i1=4503,599627,370500 i2=4503,599627,370500 i3=4503,599627,370500 i1=BAD i2=BAD i3=BAD
   Mant=$10000000000003 Exp=2^52 prec=1
i=4503,599627,370500 => i1=4503,599627,370500 i2=4503,599627,370500 i3=4503,599627,370500
   Mant=$10000000000004 Exp=2^52 prec=1
 
FPC-Version 3.2.2, 32-bit
i=4503,599627,370494 => i1=4503,599627,370494 i2=4503,599627,370494 i3=4503,599627,370494
   Mant=$1FFFFFFFFFFFFC Exp=2^51 prec=0.5
i=4503,599627,370495 => i1=4503,599627,370495 i2=4503,599627,370495 i3=4503,599627,370495
   Mant=$1FFFFFFFFFFFFE Exp=2^51 prec=0.5
i=4503,599627,370496 => i1=4503,599627,370496 i2=4503,599627,370496 i3=4503,599627,370496
   Mant=$10000000000000 Exp=2^52 prec=1
i=4503,599627,370497 => i1=4503,599627,370497 i2=4503,599627,370497 i3=4503,599627,370498 i3=BAD
   Mant=$10000000000001 Exp=2^52 prec=1
i=4503,599627,370498 => i1=4503,599627,370498 i2=4503,599627,370498 i3=4503,599627,370498
   Mant=$10000000000002 Exp=2^52 prec=1
i=4503,599627,370499 => i1=4503,599627,370499 i2=4503,599627,370499 i3=4503,599627,370500 i3=BAD
   Mant=$10000000000003 Exp=2^52 prec=1
i=4503,599627,370500 => i1=4503,599627,370500 i2=4503,599627,370500 i3=4503,599627,370500
   Mant=$10000000000004 Exp=2^52 prec=1

You see,
 - that both i1=my round_normal() and i2=math.SimpleRoundTo() fail only on 64-bit and only for odd numbers
 - that 'i3', which stores 'x+k' as an interim result, fails always (64- and 32-bit).
Same results with Optimitation Level 1 and Optimitation=off.

I assume I have an explanation, what causes this issue, but I have no explanation, why this only occurs with 64-bit and not with 32-bit.
I tested with FPC 3.2.2 on Linux 64-bit, where this issue occurs (same with FPC 3.2.0). I have no Windows 64-bit.
I tested with FPC 3.2.2 on Linux 32-bit and Windows 32-bit, where this issue not occurs (same with FPC 3.2.0 and 3.0.4).

My explanation, why this issue should always (also with 32-bit) occur, is:
 - a 'double' variable has a precision (mantissa) of 52 bits
 - as long as a value is < 2^52 = 4503,599627,370496 there are enough bits in the mantissa, that a resolution of 0.5 is storable (displayed in above output as 'prec=...')
 - but for values >= 4503,599627,370496 you can store only integers, because there are not enough bits to store a resolution less than 1
 - because of that, 4503,599627,370497 + k = 4503,599627,370497.5 is not storable and must be rounded up to 4503,599627,370498 which is too big by 1.

Question1: why does this occur only with 64-bit and not with 32-bit?

Question2: can please someone, who has Windows 64-bit, run the attached project and report the result (including your FPC version please)?

Question3: has somebody an idea to repair this issue (that all values up to const 'MaxDoublePrec1' = 9007,199254,740992 work correctly)?

Thanks in advance

[Thaddy]

Rounding is configurable.
https://www.freepascal.org/docs-html/rtl/math/setroundmode.html

Freepascal chooses the mode that applies to the ABI.
It is only a bug if it doesn't.

[marcov]

Because system.round() uses "bankers rounding", long time ago I created a "normal" rounding function:

Code: Pascal  [Select][+][-]

1. function round_normal(x: double): int64;
2.    var k: double;
3.    begin
4.    if x < 0 then k:=-0.5 else k:=0.5;
5.    exit(round(int(x+k)));
6.    end;

This worked all the time with 32-bit. Now with 64-bit I found strange rounding issues for larger values around 4.5E+15:
 - only for these larger values and
 - only for odd values and
 - only with 64-bit

Probably because the 64-bit mantissa can no longer record an 0.5 at that (15 digit) point, as it really is a double, and not a 80-bit extended under the hood.

[Thaddy]

As I already suspected, you need to set the round mode.
Your example shows rounding anomalies around values ≈ 4.5×10¹⁵, where the mantissa of a double can no longer represent increments of 0.5.
This is exactly the point where:
2⁵²= 4503599627370496
Above this threshold, double precision cannot represent every integer, only every second integer, then every fourth, etc.

Code: Pascal  [Select][+][-]

1. // do not use this
2. if x < 0 then k := -0.5 else k := 0.5;
3. round(int(x + k));

Run your - very nice code, btw -test like this:

Code: Pascal  [Select][+][-]

1. begin {main}
2.   SetRoundMode(rmDown); // add this, corrects for the above code.
3.   Test_rounding_issue;
4. end.

     
This gives nothing BAD on Win11-x86_64. The ABI specifies rmNearest. rmDown corrects this:

Code: Bash  [Select][+][-]

1. FPC-Version 3.3.1, 64-bit
2. i=4503,599627,370494 => i1=4503,599627,370494 i2=4503,599627,370494 i3=4503,599627,370494
3.    Mant=$1FFFFFFFFFFFFC Exp=2^51 prec=0.5
4. i=4503,599627,370495 => i1=4503,599627,370495 i2=4503,599627,370495 i3=4503,599627,370495
5.    Mant=$1FFFFFFFFFFFFE Exp=2^51 prec=0.5
6. i=4503,599627,370496 => i1=4503,599627,370496 i2=4503,599627,370496 i3=4503,599627,370496
7.    Mant=$10000000000000 Exp=2^52 prec=1
8. i=4503,599627,370497 => i1=4503,599627,370497 i2=4503,599627,370497 i3=4503,599627,370497
9.    Mant=$10000000000001 Exp=2^52 prec=1
10. i=4503,599627,370498 => i1=4503,599627,370498 i2=4503,599627,370498 i3=4503,599627,370498
11.    Mant=$10000000000002 Exp=2^52 prec=1
12. i=4503,599627,370499 => i1=4503,599627,370499 i2=4503,599627,370499 i3=4503,599627,370499
13.    Mant=$10000000000003 Exp=2^52 prec=1
14. i=4503,599627,370500 => i1=4503,599627,370500 i2=4503,599627,370500 i3=4503,599627,370500
15.    Mant=$10000000000004 Exp=2^52 prec=1

rmTruncate also gives correct result in the test, but rmDown is mathematically more correct and is what you expect. Note that for financial applications rmNearest is correct. Just not for - some - mathematics.
SetRoundMode can be set on a per routine basis or use math.RoundTo.
Use rmDown for monotonic rounding (deterministic)
Use rmNearest for financial, statistical or probabilistic computations

Note rmNearest is specified for all x86_64 ABI's, not just Windows, as per IEEE754

* rmDown equals Floor for positive numbers and more negative for negative. School's rounding in many countries.
* rmNearest is bankers rounding, accumulates less error in long calculations.
* rmUp equals Ceil for positive numbers else round towards zero. Never underestimates, think current, voltage, capacitance, resistance in electronics design.
* rmTruncate equals int(value).

(Text slightly editted later with the help of CoPilot, fact checks and suggestions. Original text was accurate.)

[avk]

...
Question1: why does this occur only with 64-bit and not with 32-bit?
...

The key point here seems to be that in the expression

Code: Pascal  [Select][+][-]

1.   round(int(x+k))

the result of adding x + k must be rounded to a representable value according to the current rounding mode.
In 32-bit RTL, by default, calculations are performed in pmExtended(80 bit) precision mode, so all values ​​up to 9007199254740992 + 0.5 are representable.
In Win64, calculations are performed with 64-bit precision, 4503599627370495 + 0.5 is still representable as a Double, but 4503599627370496 + 0.5 is no longer representable (for 4503599627370496, the next representable value is 4503599627370497) and is rounded accordingly.
If you set the precision mode to pmDouble for 32-bit RTL, results should be the same.

...
Question3: has somebody an idea to repair this issue (that all values up to const 'MaxDoublePrec1' = 9007,199254,740992 work correctly)?
...

Maybe

Code: Pascal  [Select][+][-]

1. function round_normal(d: Double): Int64;
2. begin
3.   if Frac(d) = 0 then exit(Trunc(d));
4.   if d < 0 then
5.     Result := Trunc(d - 0.5)
6.   else
7.     Result := Trunc(d + 0.5);
8. end;
9.  

[Thaddy]

@AVK
That is the same as leaving the .5 comparison in and SetRoundMode(rmDown)?

I thoroughly tested and verified my comment.
To many rounding is a rat's nest anyway and I tried to explain what and where you should apply different roundings. And what I wrote is 100% correct.

[Hartmut]

Probably because the 64-bit mantissa can no longer record an 0.5 at that (15 digit) point, as it really is a double, and not a 80-bit extended under the hood.

Hmm...
 - the issue is with 'double' (52 bit mantissa) and not 'extended' (64 bit mantissa)
 - a 64-bit mantissa can store a 0.5 resolution for values like 4503,599627,370497, only a 52 bit mantissa can not
 - on Linux 64-bit 'extended' really exists (I checked: sizeof(extended)=10), so it must not be substituted by a 'double'

Run your - very nice code, btw -test like this:

Code: Pascal  [Select][+][-]

1. begin {main}
2.   SetRoundMode(rmDown); // add this, corrects for the above code.
3.   Test_rounding_issue;
4. end.

'rmDown' is in https://lazarus-ccr.sourceforge.io/docs/rtl/math/setroundmode.html defined as "Round to biggest integer *smaller* than value". The results for 2.75 with 'rmDown' are:
 - system.Round(2.75) => 2
 - math.RoundTo(2.75) => 2
 - math.SimpleRoundTo(2.75) => 3
So obviously 'rmDown' *changes* the behaviour of system.Round() and math.RoundTo(), but obviously does *not* change the behaviour of math.SimpleRoundTo().
So obviously I *cannot* use system.Round() and math.RoundTo() with 'rmDown', because I want "normal" rounding.

You can test this with this code:

Code: Pascal  [Select][+][-]

1. procedure test_RoundMode;
2.    {shows 3 round methods with current RoundMode}
3.  
4.    procedure test(d: double);
5.       begin
6.       writeln('d=', d:0:2, '  Round=', system.Round(d),
7.       '  RoundTo=', math.RoundTo(d,0):0:2,
8.       '  SimpleRoundTo=', math.SimpleRoundTo(d,0):0:2);
9.       end;
10.  
11.    begin
12.    writeln('GetRoundMode=', GetRoundMode);
13.    test(2.5);
14.    test(2.75);
15.    test(4503599627370494.5);
16.    test(4503599627370495.5); {highest value with resolution of 0.5}
17.    test(4503599627370496.0);
18.    test(4503599627370497.0);
19.    end;
20.  
21. begin
22. writeln('FPC-Version ', {$I %FPCVERSION%}, ', ',
23.         {$IFDEF CPU32} '32-bit' {$ELSE} '64-bit' {$ENDIF} );
24. test_RoundMode; {use default RoundMode}
25. writeln;
26. SetRoundMode(rmDown);
27. test_RoundMode;
28. end.

The output is:
FPC-Version 3.2.2, 64-bit
GetRoundMode=rmNearest
d=2.50  Round=2  RoundTo=2.00  SimpleRoundTo=3.00
d=2.75  Round=3  RoundTo=3.00  SimpleRoundTo=3.00
d=4503599627370494.50  Round=4503599627370494  RoundTo=4503599627370494.00  SimpleRoundTo=4503599627370495.00
d=4503599627370495.50  Round=4503599627370496  RoundTo=4503599627370496.00  SimpleRoundTo=4503599627370496.00
d=4503599627370496.00  Round=4503599627370496  RoundTo=4503599627370496.00  SimpleRoundTo=4503599627370496.00
d=4503599627370497.00  Round=4503599627370497  RoundTo=4503599627370497.00  SimpleRoundTo=4503599627370498.00

GetRoundMode=rmDown
d=2.50  Round=2  RoundTo=2.00  SimpleRoundTo=3.00
d=2.75  Round=2  RoundTo=2.00  SimpleRoundTo=3.00
d=4503599627370494.50  Round=4503599627370494  RoundTo=4503599627370494.00  SimpleRoundTo=4503599627370495.00
d=4503599627370495.50  Round=4503599627370495  RoundTo=4503599627370495.00  SimpleRoundTo=4503599627370496.00
d=4503599627370496.00  Round=4503599627370496  RoundTo=4503599627370496.00  SimpleRoundTo=4503599627370496.00
d=4503599627370497.00  Round=4503599627370497  RoundTo=4503599627370497.00  SimpleRoundTo=4503599627370497.00

But I can approve your results: with 'rmDown' all results for the "large values" in my 1st demo on Linux 64-bit are correct. But why ???

Code: Pascal  [Select][+][-]

1. procedure test_round(i: int64);
2. ...
3.    i1:=round(int(x+k));                {rounding method #1}
4.    i2:=round(math.SimpleRoundTo(x,0)); {rounding method #2}
5.    h:=x+k; i3:=round(int(h));          {rounding method #3 with interim result}
6. ...

 - 'i1' should always fail, because 'x+k' = 4503,599627,370497 + 0.5 = 4503,599627,370497.5 is not storable and must be rounded up to 4503,599627,370498 => function int() of that value should always return 4503,599627,370498 (which is wrong)
 - 'i2' is based on the same formula => should always behave the same => should always fail
 - 'i3' stores 'x+k' in an interim variable => should always behave like 'i1' and 'i2' => should always fail

The code for math.SimpleRoundTo() in FPC 3.2.2 is:

Code: Pascal  [Select][+][-]

1. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange = -2): Double;
2. var
3.   RV : Double;
4. begin
5.   RV := IntPower(10, -Digits);
6.   if AValue < 0 then
7.     Result := Int((AValue*RV) - 0.5)/RV
8.   else
9.     Result := Int((AValue*RV) + 0.5)/RV;
10. end;

You see, it's based on the same formula as 'i1', so it *should fail* the same way.

Because math.SimpleRoundTo() works correctly with 'rmDown' for "large values" on Linux 64-bit, but *does not* work with default 'rmNearest', although it's rounding behaviour is *not* changed (as shown above in procedure test_RoundMode), there must be a "hidden" side affect. The code for math.SetRoundMode() in FPC 3.2.2 is:

Code: Pascal  [Select][+][-]

1. function SetRoundMode(const RoundMode: TFPURoundingMode): TFPURoundingMode;
2. var
3.   CtlWord: Word;
4.   SSECSR: dword;
5. begin
6.   CtlWord:=Get8087CW;
7.   SSECSR:=GetMXCSR;
8.   Set8087CW((CtlWord and $F3FF) or (Ord(RoundMode) shl 10));
9.   SetMXCSR((SSECSR and $ffff9fff) or (dword(RoundMode) shl 13));
10. {$ifdef FPC_HAS_TYPE_EXTENDED}
11.   Result:=TFPURoundingMode((CtlWord shr 10) and 3);
12. {$else}
13.   Result:=TFPURoundingMode((SSECSR shr 13) and 3);
14. {$endif FPC_HAS_TYPE_EXTENDED}
15. end;

I understand nothing from that code. How can SetRoundMode(rmDown) "repair" that rounding bug in math.SimpleRoundTo(), although it's rounding behaviour is not changed?? And: math.SimpleRoundTo() does not use any round() function!

So still open questions are:
 - why does SetRoundMode(rmDown) "repair" the rounding issue for 'i1' and 'i2' and 'i3'? All three should always fail!
 - why does the rounding issue occur only with 64-bit and never with 32-bit?

In 32-bit RTL, by default, calculations are performed in pmExtended(80 bit) precision mode, so all values ​​up to 9007199254740992 + 0.5 are representable.
In Win64, calculations are performed with 64-bit precision...

This sounds very interesting, especially your explanation for 32-bits and your repair suggestion. Thank you. But I saw your post just now and need some time to check it more detailed. Will report later.

If you set the precision mode to pmDouble for 32-bit RTL, results should be the same.

Is this possible? How please?

[Thaddy]

- on Linux 64-bit 'extended' really exists (I checked: sizeof(extended)=10), so it must not be substituted by a 'double'

That is only because Freepascal implemented it (8087). It is not part if the Linux64 ABI, nor is it for 64 bit Mac nor for Windows64.
This is due to the fact that these ABI's use SSEx/AVXx instructions for math operations and these are strictly 8 wide.
So beware your mileage may vary...depending on optimization settings. That part I did not explain, (I did just now, though) otherwise follow my instructions. They are 100% correct.

You want a solution? I gave you the correct solution + explanation. Don't even think of using the word but... It is not applicable. I am not even commenting later. With such questions it is useless.

[avk]

...
Is this possible? How please?

Easy

Code: Pascal  [Select][+][-]

1. uses
2.   ..., Math;
3. ...
4.   SetPrecisionMode(pmDouble);
5. ...
6.  

@AVK
That is the same as leaving the .5 comparison in and SetRoundMode(rmDown)?
...

Looks similar.

[Hartmut]

Question1: why does this occur only with 64-bit and not with 32-bit?

The key point here seems to be that in the expression

Code: Pascal  [Select][+][-]

1.   round(int(x+k))

the result of adding x + k must be rounded to a representable value according to the current rounding mode.
In 32-bit RTL, by default, calculations are performed in pmExtended(80 bit) precision mode, so all values up to 9007199254740992 + 0.5 are representable.
In Win64, calculations are performed with 64-bit precision, 4503599627370495 + 0.5 is still representable as a Double, but 4503599627370496 + 0.5 is no longer representable (for 4503599627370496, the next representable value is 4503599627370497) and is rounded accordingly.
If you set the precision mode to pmDouble for 32-bit RTL, results should be the same.

You are right. I tested with my very 1st demo with math.SetPrecisionMode(pmDouble) and then
 - in WIN 32bit and Linux 32bit for 4503,599627,370497 and 4503,599627,370499 additinally "i1=BAD i2=BAD" is shown (so the result is completely the same as in Linux 64bit)
 - in Linux 64bit SetPrecisionMode(pmDouble) makes no difference. Although math.GetPrecisionMode() returns 'pmExtended', it seems to really work with 'pmDouble'.

Question:
Does this mean: In WIN 32bit and Linux 32bit not only the calculation of 'x+k' is done with extended (80bit) precision - then the result is passed to the int() function also with 80bit precision? - without converting the result of 'x+k' between to double?
If yes, I would understand and would have learned a lot.

Question3: has somebody an idea to repair this issue (that all values up to const 'MaxDoublePrec1' = 9007,199254,740992 work correctly)?

Maybe

Code: Pascal  [Select][+][-]

1. function round_normal(d: Double): Int64;
2. begin
3.   if Frac(d) = 0 then exit(Trunc(d));
4.   if d < 0 then
5.     Result := Trunc(d - 0.5)
6.   else
7.     Result := Trunc(d + 0.5);
8. end;

As far as I could test, this function works perfectly. Thank you very much! I prefer it much more than Thaddys solution, because
 - yours is easier, shorter and faster than math.SimpleRoundTo() - see it's source in reply #6
 - and SetRoundMode(rmDown) is not neccessary, which had to be read before and restored after, because otherwise all calls to system.Round() and math.RoundTo() in the rest of the whole project (including foreign libraries) would return wrong values.

- on Linux 64-bit 'extended' really exists (I checked: sizeof(extended)=10), so it must not be substituted by a 'double'

That is only because Freepascal implemented it (8087). It is not part if the Linux64 ABI, nor is it for 64 bit Mac nor for Windows64.
This is due to the fact that these ABI's use SSEx/AVXx instructions for math operations and these are strictly 8 wide.

I'm not familiar with ABI, SSEx and AVXx stuff. So please what are the consequences of what you write?
 - does FPC on Linux 64bit have variables of type 'extended'?
 - can FPC on Linux 64bit store values of type 'extended' in those variables correctly?
 - does FPC on Linux 64bit all calculations (at least where the input comes from 'extended' variables and the result is stored in 'extended' variables) also with 'extended' (80bit) precision?
 - or does FPC on Linux 64bit all those calculations only with 'double' precision?

And the last 2 open questions are:
 - why does the rounding issue not occur on Thaddys WIN 64bit (see reply #3)? I learned that on WIN 64 'extended' is not available and is replaced by 'double'?
 - why does SetRoundMode(rmDown) "repair" the rounding issue in math.SimpleRoundTo() on Linux 64bit, although it's rounding behaviour is not changed by that? Does SetRoundMode(rmDown) have a hidden "side affect"?

[Thaddy]

A standard 64 bit build:(ppcx64 -if on linux64)
Supported FPU instruction sets:
  NONE,SSE64,X86-64-V1,SSE3,SSSE3,SSE41,SSE42,X86-64-V2,AVX,FMA,AVX2,
  X86-64-V3,AVX512F,X86-64-V4
Does not contain 8087! None of these are 10 wide.
A standard 32 bit build (ppc386:
Supported FPU instruction sets:
  NONE,X87,SSE,SSE2,SSE3,SSSE3,SSE41,SSE42,AVX,FMA,AVX2,AVX512F
But you can manually build a FPC64/PPCX64 to add that x87 support because Florian himself added that option.

[edit] Although x87 is not listed, this small program indicates that in Linux64 the extended type is 10.

Code: Pascal  [Select][+][-]

1. begin writeln(SizeOf(extended));end.// writes 8 on win64, 10 on linux64

Some additional research and a quirk you should be aware of: this is NOT an ABI violation, like I first thought: FPC64 does not pass extended in registers, like it does on i386, but passes extended by reference on Linux64 and that is allowed in the ABI.
It is only not allowed to pass extended values in (x87) registers, like is done on i386. Internally it is free to use 8087 instructions
Learned something new here, I really wasn't aware of that and the underlying mechanism.

FreePascal’s 80‑bit extended on Linux/x86‑64 is not an ABI violation.
The SysV ABI only forbids passing 80‑bit values in x87 registers.
FPC64 passes extended by reference, which is allowed, and uses x87 internally for extended expressions.
So the compiler is fully ABI‑compliant while still supporting true 80‑bit precision.

This is technically precise and matches GCC’s behavior.

Quote by CoPilot. Explains quite a lot about the precision.

[Hartmut]

A standard 64 bit build:(ppcx64 -if on linux64)
Supported FPU instruction sets:
  NONE,SSE64,X86-64-V1,SSE3,SSSE3,SSE41,SSE42,X86-64-V2,AVX,FMA,AVX2,
  X86-64-V3,AVX512F,X86-64-V4
Does not contain 8087! None of these are 10 wide.
A standard 32 bit build (ppc386:
Supported FPU instruction sets:
  NONE,X87,SSE,SSE2,SSE3,SSSE3,SSE41,SSE42,AVX,FMA,AVX2,AVX512F

As I wrote, I'm not familiar with this ABI, FPU, SSEx, AVXx etc. stuff.
I understand nothing from that what Thaddy writes.
Obviously is Thaddy not able to answer my questions with a clear understandable yes or no.
It makes no sense to deduce the answers (on the basis of those difficult details) by my self, because my answers most probably would be wrong.

Can please someone else help and answer these questions with a clear understandable yes or no?
 - does FPC on Linux 64bit have variables of type 'extended'?
 - can FPC on Linux 64bit store values of type 'extended' in those variables correctly?
 - does FPC on Linux 64bit compute all calculations (at least where the input comes from 'extended' variables and the result is stored in 'extended' variables) also with 'extended' (80bit) precision?
 - or does FPC on Linux 64bit compute all those calculations only with 'double' precision?

[LeP]

We're always talking about 80-bit extended values. This format was introduced by Intel with its x87 coprocessor (and also by Motorola, if I recall correctly) for the sole purpose of ensuring correct intermediate processing between calculations.

INTEL has always discouraged and deprecated the direct use of the extended format.

Indeed, with 64-bit technology, this format was no longer officially supported (for other reasons as well).

Continuing to use and support this format, except for legacy applications due to historical reasons, is a high-risk activity.

Support for x87 and its very "presence" could cease at any time in the near future.

New applications or portings should not use the 80-bit extended format. If you want to use a mathematical format for complex calculations, where precision is crucial, you should use specific libraries and not rely on deprecated hardware.

[Hartmut]

Hello LeP, I do not want to use 'extended'.
I only want to understand why some rounding issues with type 'double' occur differently in 64bit and 32bit.
Therefore some questions came up, which precision FPC currently uses to computes (internal) calculations on Linux 64bit.

It would be great, if someone except Thaddy (because he is not able to answer with yes or no) could answer my questions with a clear understandable yes or no:
 - does FPC on Linux 64bit have variables of type 'extended'?
 - can FPC on Linux 64bit store values of type 'extended' in those variables correctly?
 - does FPC on Linux 64bit compute all calculations (at least where the input comes from 'extended' variables and the result is stored in 'extended' variables) also with 'extended' (80bit) precision?
 - or does FPC on Linux 64bit compute all those calculations only with 'double' precision?

[Thaddy]

Can please someone else help and answer these questions with a clear understandable yes or no?
 - does FPC on Linux 64bit have variables of type 'extended'?
 - can FPC on Linux 64bit store values of type 'extended' in those variables correctly?
 - does FPC on Linux 64bit compute all calculations (at least where the input comes from 'extended' variables and the result is stored in 'extended' variables) also with 'extended' (80bit) precision?
 - or does FPC on Linux 64bit compute all those calculations only with 'double' precision?

Sorry but you did not read me correctly, so I answer (correctly):
1. Yes full 80 bit width internally, passed by reference, not via fpu registers as in i386
2. Yes
3. Yes but as reference, not register, the compiler references the result it calculates internally. No fpu registers involved in the return values.
4. No, but Windows64 does compile only in double. On Windows64, extended is an alias to double.

I suggest to compile with -al and examine the .s file to see how it works.
Given:

Code: Pascal  [Select][+][-]

1. // only for linux intel/amd 64 bit
2. {$mode objfpc}
3. var
4.   a,b:extended;
5. begin
6.   a := 0.3;
7.   b := 10.0;
8.   a := b/ a;
9.   writeln(a:1:17);
10. end.

Outputs for Linux64 Intel/AMD.

Code: ASM  [Select][+][-]

1. PASCALMAIN: # this is just the essential fragment
2. # [8] a := b/ a;
3.         fldt    U_$P$TESTEXTENDED_$$_A  # load by reference into fpu 8087 space
4.         fldt    U_$P$TESTEXTENDED_$$_B  # using 8087 instructions
5.         fdivp   %st,%st(1)              # perform 8087 division, 80 bits
6.         fstpt   U_$P$TESTEXTENDED_$$_A  # store in non - fpu 8087 space
7. # [9] writeln(a:1:17);
8.         call    fpc_get_output

If that is not clear enough, then someone else should answer it. You could have done that yourself.
I clearly commented what happens here. This is slower than the 32 bit i386 compiler.