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

[H₂SO₄]

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.

If you had taken the time to familiarize yourself even just a bit with how floating-point operations are implemented on x86 instead of bitching about how Thaddy didn't dumb down his answers enough, you'd have understood by now that there's a global switch that influences every floating-point operation on a particular CPU core, and that on different platforms this switch is set to different positions at program startup.

Take a look at the  System  unit on Linux:

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

1. procedure SysEntry(constref info: TEntryInformation);[public,alias:'FPC_SysEntry'];
2. begin
3.   SetupEntryInformation(info);
4. {$ifdef cpui386}
5.   Set8087CW(Default8087CW);
6. {$endif cpui386}
7.   info.PascalMain();
8. end;

That  {$ifdef}  subtly changes how floats work depending on whether the CPU is 32-bit or not. You're stuck with that uncertainty unless you explicitly tell the CPU how to behave (via  SetRoundMode,  SetPrecisionMode  et al.).

[Hartmut]

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.

NOW your answers are clear understandable. Thanks.

Because in the past answers from Thaddy were often (completely or partly) wrong:
can someone else confirm, whether these 4 answers are correct?

If ... you'd have understood by now that there's a global switch that influences every floating-point operation on a particular CPU core, and that on different platforms this switch is set to different positions at program startup.

I checked your Link, which is about a bug in Visual C++ 2008.
How should I know, whether this has any relevance to current Free Pascal on LINUX?
And even if I would have known: I found there nothing, what could answer any of the still open questions.

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

1. procedure SysEntry(constref info: TEntryInformation);[public,alias:'FPC_SysEntry'];
2. begin
3.   SetupEntryInformation(info);
4. {$ifdef cpui386}
5.   Set8087CW(Default8087CW);
6. {$endif cpui386}
7.   info.PascalMain();
8. end;

I already knew, that there is a difference in float handling between 32bit and 64bit.
But this source does not answer any of the still open questions.

Maybe you are an expert for all that stuff - maybe you forgot, that there are people, who's knowledge about that is near zero.
With this kind of post you don't help me.

You're stuck with that uncertainty unless you explicitly tell the CPU how to behave (via  SetRoundMode,  SetPrecisionMode  et al.).

I already examined SetRoundMode and SetPrecisionMode and reported the results above. Didn't you read that?

This are the 2 left open questions:
 - 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"?

[LeP]

@Hartmut, I use Delphi and the result is identical to @Thaddy test: works in WIN32 bit mode or WIN64 mode ad the same wtihout "BAD" signals.

I see that for Math.SimpleRoundTo there is a commente in Delphi side:

{ This variation of the RoundTo function follows the symmetric arithmetic
  rounding algorithm ....... (the same as FPC) ..........
  Note: FPU rounding modes affect the behavior of this function. }

So rounding mode can affect the results for SimpleRoundTo .... in Delphi, but in FPC to I think ....

[Hartmut]

@Hartmut, I use Delphi and the result is identical to @Thaddy test: works in WIN32 bit mode or WIN64 mode ad the same wtihout "BAD" signals.

Thanks LeP for your test.
Has anybody an explanation, why the rounding issue does *not* occur on WIN 64bit, although on WIN 64bit 'extended' is not available and is replaced by 'double', so the issue *should* occur?

{ This variation of the RoundTo function follows the symmetric arithmetic
  rounding algorithm ....... (the same as FPC) ..........
  Note: FPU rounding modes affect the behavior of this function. }

So rounding mode can affect the results for SimpleRoundTo .... in Delphi, but in FPC to I think ....

I showed in reply #6, that SetRoundMode() does change the rounding behaviour *only* of system.Round() and math.RoundTo(), but does *not* change the rounding behaviour of math.SimpleRoundTo().
In reply #6 I showed this only for 'rmNearest' and 'rmDown', but I tested 'rmUp' and 'rmTruncate' too and both do not change the rounding behaviour of math.SimpleRoundTo().

What I do not understand: why does SetRoundMode(rmDown) "repair" the rounding issue in math.SimpleRoundTo() on Linux 64bit?
Because SetRoundMode(rmDown) does not change the "rounding mode".
Might SetRoundMode(rmDown) have a hidden "side affect"?

[Thaddy]

You are very difficult ro follow:
SetRoundMode does not "repair" but change the rounding mode.
SimpleRoundTo? That uses integer rounding based on IntPower so bypasses math processing altogether. Don't ever use that in serious programs.
It has only one way of rounding.

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

1. // don't use this in real programs:
2. {$ifdef FPC_HAS_TYPE_DOUBLE}
3. function SimpleRoundTo(const AValue: Double; const Digits: TRoundToRange = -2): Double;
4. var
5.   RV : Double;
6. begin
7.   RV := IntPower(10, -Digits);
8.   if AValue < 0 then
9.     Result := Int((AValue*RV) - 0.5)/RV
10.   else
11.     Result := Int((AValue*RV) + 0.5)/RV;
12. end;
13. {$endif}
14.  
15. {$ifdef FPC_HAS_TYPE_EXTENDED}
16. function SimpleRoundTo(const AValue: Extended; const Digits: TRoundToRange = -2): Extended;
17. var
18.   RV : Extended;
19. begin
20.   RV := IntPower(10, -Digits);
21.   if AValue < 0 then
22.     Result := Int((AValue*RV) - 0.5)/RV
23.   else
24.     Result := Int((AValue*RV) + 0.5)/RV;
25. end;
26. {$endif}

Why don't you simply check the sourcecode? The above IS the sourcecode.
I am very polite according to @Khrys, (that is a novelty), maybe I should change that.
You make not any effort yourself. You did not read the sourcecode in this case.
You make no meaningful effort. Simply stop programming, that is better for the rest of the world. You make yourself look silly.
(I  am a silly old man anyway, but I can actually program and I can also read and investigate: you do not have these additional properties. OTOH You can moan very good. I envy you )

[Hartmut]

Why don't you simply check the sourcecode?
...
You did not read the sourcecode in this case.

As so often you are wrong. Just yesterday I wrote:

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.

As so often you did not read, before you post something, complain and become rude.

You make not any effort yourself. You did not read the sourcecode in this case.
You make no meaningful effort. Simply stop programming, that is better for the rest of the world. You make yourself look silly.
(I  am a silly old man anyway, but I can actually program and I can also read and investigate: you do not have these additional properties.)

As so often you become impertinent.
The whole forum knows this.
Do a Forum Search for "Thaddy rude". You'll get 20 (twenty) matches!

[Thaddy]

I am not wrong and everybody else will agree.

[H₂SO₄]

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.

Yes:

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

1. lea  esp, [esp - 12]
2. fstp TBYTE PTR [esp]
3. call <fpc_int_real>

This is the disassembly of the 32-bit Windows code that passes the result of  x+k  to  fpc_int_real  with all 80 bits intact (as evidenced by the presence of  TBYTE PTR).

If ... you'd have understood by now that there's a global switch that influences every floating-point operation on a particular CPU core, and that on different platforms this switch is set to different positions at program startup.

I checked your Link, which is about a bug in Visual C++ 2008.
How should I know, whether this has any relevance to current Free Pascal on LINUX?
And even if I would have known: I found there nothing, what could answer any of the still open questions.

I linked to this specific answer in that thread (which I basically restated in the rest of my answer):

When the OS sets up a process (or thread) -- say, when your program launches, it is responsible for initializing the control registers to the default state (exactly what the default state is can vary somewhat from platform to platform). It is also responsible for saving and restoring state around context switches, so that your program doesn't bleed state into some other process that is also running on the system.

Beyond that, a compiler or a language runtime might modify the state before your code is executed.

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

1. procedure SysEntry(constref info: TEntryInformation);[public,alias:'FPC_SysEntry'];
2. begin
3.   SetupEntryInformation(info);
4. {$ifdef cpui386}
5.   Set8087CW(Default8087CW);
6. {$endif cpui386}
7.   info.PascalMain();
8. end;

I already knew, that there is a difference in float handling between 32bit and 64bit.
But this source does not answer any of the still open questions.

Wasn't your first question about wanting to understand why there are platform differences...?

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

Well? There you go.

You're stuck with that uncertainty unless you explicitly tell the CPU how to behave (via  SetRoundMode,  SetPrecisionMode  et al.).

I already examined SetRoundMode and SetPrecisionMode and reported the results above. Didn't you read that?

Yes, and you seem to have concluded that even though they solve the inconsistencies you still don't want to use them.

- 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"?

(Emphasis mine)

SetRoundMode  calls  Set8087CW,  so yes, it does have a side effect - across your entire program and every single floating-point operation it performs.
This function does not just influence  Round  or  SimpleRoundTo  in particular.

Maybe you are an expert for all that stuff - maybe you forgot, that there are people, who's knowledge about that is near zero.
With this kind of post you don't help me.

Your attitude irritates me. Seemingly eager to learn, you ask questions about a complicated topic only then to complain when - surprise - the answers are also quite complicated due to the nature of the topic. But instead of putting in even a modicum of effort to educate yourself, you expect to be spoonfed all the way.

Your knowledge is near zero? So was mine at some point, but at least I put in the effort to learn.

[Hartmut]

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.

Yes:

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

1. lea  esp, [esp - 12]
2. fstp TBYTE PTR [esp]
3. call <fpc_int_real>

This is the disassembly of the 32-bit Windows code that passes the result of  x+k  to  fpc_int_real  with all 80 bits intact (as evidenced by the presence of  TBYTE PTR).

Thanks a lot for this clear and helpful answer. As you see, I'm not only interested to find a solution (which I yet have 2 different since reply #4), I'm also interested to learn something and invest a lot of time for this, because this will become reply #23.

I linked to this specific answer in that thread (which I basically restated in the rest of my answer):

When the OS sets up a process (or thread) -- say, when your program launches, it is responsible for initializing the control registers to the default state (exactly what the default state is can vary somewhat from platform to platform). It is also responsible for saving and restoring state around context switches, so that your program doesn't bleed state into some other process that is also running on the system.
Beyond that, a compiler or a language runtime might modify the state before your code is executed.

As said, your Link was about Microsoft Visual C++ 2008 on Windows. I could not see, whether this has any relevance for my issue on Linux.

Wasn't your first question about wanting to understand why there are platform differences...?

You are right, my question was not precise enough. You showed me, *where* the switch is done, which causes the difference between 32bit and 64bit. My question was more intended to understand the differences in the *behaviour* between 32bit and 64 bit, when/how (internal) calculations are computed. That's on me asking not clear enough, sorry for confusion. But due to the answers of avk and you now I understand.

- 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"?

SetRoundMode  calls  Set8087CW,  so yes, it does have a side effect - across your entire program and every single floating-point operation it performs.
This function does not just influence  Round  or  SimpleRoundTo  in particular.

That sounds interesting.
I checked the official documentation for 'Set8087CW' in https://www.freepascal.org/docs-html/current/rtl/system/set8087cw.html which only says:

Set8087CW
Declaration
Source position: mathh.inc line 21
procedure Set8087CW(
  cw: Word);

And I checked the sources both in line 21 of <installdir>/fpcsrc/rtl/inc/mathh.inc and line 51 in <installdir>/fpcsrc/rtl/x86_64/math.inc but this is ASM-Code beyond my horizon.

BTW: I did this check for 'Set8087CW' also 2 days ago, before writing reply #6, before posting there the source of function SetRoundMode(). So your criticism

But instead of putting in even a modicum of effort to educate yourself, you expect to be spoonfed all the way.

is not applicable.

If you already know more about this "side affect" caused by Set8087CW() - what are those differences for "the entire program and every single floating-point operation it performs" - and if you want to write something about that, I would appreciate this.

But in this Forum nobody is forced to write something, if he does not want to. Even if he knows the answer, but only assumes, that for his criterion the questioner does not invest enough effort, to answer his questions by himself.
Since I found this rounding issue 6 days ago, I spent 100% of my freetime for it. First trying to understand and solve it by myself (in my 1st post you can see, how much I investigated and found out by myself), than starting and maintaining this Topic. Maybe you have unlimited freetime. Unfortunately I have not.

[Thaddy]

Well it is not the entire program if you save the round mode and restore it.
Clear demo:

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

1. {$mode objfpc}
2. uses math;
3. function RoundDown(const value:extended):integer;inline;
4. var
5.   OldMode:TFPURoundingMode;
6. begin
7.   OldMode := GetRoundMode;
8.   SetRoundMode(rmDown);
9.   Result := Round(value);
10.   SetRoundMode(OldMode);
11. end;
12.  
13. begin
14.   writeln(Round(5.55554):5);    // print 6
15.   writeln(RoundDown(5.5554):5); // prints 5
16.   writeln(Round(-5.55554):5);   // prints -6
17.   writeln(RoundDown(-5.5554):5);// prints -6 too, you may figure that one out...;)
18. end.

That does not take away from the fact that you will loose precision on 64 bit Windows, certainly with the +/- .5 trick. But that is explained.
Otherwise it is "Pearls before swine" all over again.
This shows you can manipulate round mode on a per routine basis.
If you can answer the question in the sourcecode, you have come a long way...

Btw: x86_64-win64 renders:

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

1. # [10] Result := Round(value);
2.         vmovsd  -8(%rbp),%xmm0  # a double!!, not extended.
3.         vcvtsd2si       %xmm0,%rax
4.         movl    %eax,-12(%rbp)

And these registers are not 80 bit.....Not what we already showed you for linux64 an 32 bit.
This is Windows64 ABI compliant and not the fpc compiler that dictates this.

[Hartmut]

@Hartmut, I use Delphi and the result is identical to @Thaddy test: works in WIN32 bit mode or WIN64 mode ad the same wtihout "BAD" signals.

Hello @LeP,
can you remember, whether you did your test on WIN64 (where the rounding issue did not occur) with my unchanged demo? Or did you call SetRoundMode(rmDown) before this test?

If you did your test on WIN64 with calling SetRoundMode(rmDown) before: could you please be so kind to repeat your test on WIN64 once without calling SetRoundMode(rmDown) before? Does the rounding issue then occur?
Thanks.

[LeP]

@Hartmut, I use Delphi and the result is identical to @Thaddy test: works in WIN32 bit mode or WIN64 mode ad the same wtihout "BAD" signals.

Hello @LeP,
can you remember, whether you did your test on WIN64 (where the rounding issue did not occur) with my unchanged demo? Or did you call SetRoundMode(rmDown) before this test?

If you did your test on WIN64 with calling SetRoundMode(rmDown) before: could you please be so kind to repeat your test on WIN64 once without calling SetRoundMode(rmDown) before? Does the rounding issue then occur?
Thanks.

I will test again with all vars ...

[Thaddy]

Yes, he did, as did I. Pearles before swine. I am out of here. You have the ffíng assembler output. Why do you question us?

[LeP]

Delphi x32 without touch your code (except for compatibility):

Delphi-Version 37.0, 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

Delphi x64 without touch your code (except for compatibility):

Delphi-Version 37.0, 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

With SetRoundMode(rmDown) it's all right without any bad signals.

[Hartmut]

Thank you very much LeP for those detailed results. It shows, that WIN and Linux behave completely identical!
 - without SetRoundMode(rmDown) 32bit is correct and 64bit is wrong (the latter was unknown before this test)
 - with SetRoundMode(rmDown) 64bit becomes correct too.