#### View Issue Details

ID Project Category View Status Date Submitted Last Update 0033167 FPC RTL public 2018-02-13 19:30 2019-12-09 13:09 wp normal minor have not tried new reopened 3.1.1 3.1.1 0033167: Calculation errors of floating point "mod" operation The unit math of fpc-trunk implements the "mod" operator for floating point numbers. This is coded as   operator mod(const a,b:float) c:float;inline;   begin     c:= a-b * Int(a/b);   end; However, since the integer part of ratio a/b may be off by 1 the result can be off by one divisor if the result is expected to be close to 0. Example:   a = 7.7   b = 1.1 Simple math tells that a mod b should be 0 here because 7.7/1.1 = 7.0. But due to the finite precision of floating point values the ratio actually turns out to be 6.9999999. Together with the Int() operation the result of the mod operation is off by 1 divisor. However, since the remainder of a division is always less than the divisor this error case can be detected by checking c and b for equality with the function SameValue which takes machine precision into account. Furthermore, the function abs should be used to cover the case of negative values.   operator mod(const a,b:float) c:float;inline;   begin     c:= a-b * Int(a/b);     if SameValue(abs(c), abs(b)) then c := 0.0;   end; Run this project: program Project1; uses   Math; operator mod(const a,b:Double) c:Double; inline; begin   c:= a-b * Int(a/b);   //if SameValue(abs(c), abs(b)) then c := 0; // Uncomment to see effect of fix end; procedure Test(Dividend, Divisor: double); var   res: Double; begin   WriteLn('Dividend: ', Dividend:10:5);   WriteLn('Divisor : ', Divisor:10:5);   WriteLn('Dividend/Divisor = ', Dividend/Divisor);   res := Dividend mod Divisor;   WriteLn('Result: ', res:10:5);   WriteLn; end; begin   Test(7.7, 1.1);   Test(-7.7, 1.1);   Test(7.7, -1.1);   Test(-7.7, -1.1);   {   Test(7.2, 1.1);   Test(7.2, -1.1);   Test(-7.2, 1.1);   Test(-7.2, -1.1);   }   ReadLn; end. I am not sure if negative values are handled by the fpc code correctly, anyway. I compared with Wolfram Alpha output:  7.2 mod 1.1 = 0.6 (Wolfram Alpha), same here  7.2 mod -1.1 = -0.5 (Wolfram Alpha), but 0.6 here  -7.2 mod 1.1 = 0.5 (Wolfram Alpha), but -0.6 here  -7.2 mod -1.1 = -0.6 (Wolfram Alpha), same here To me it looks as if there is no clear definition how to handle negative values: https://en.wikipedia.org/wiki/Modulo_operation No tags attached. 38332 Attached Files tw33167.pp (1,684 bytes) ```{\$ifndef FPUNONE} uses math; operator mod(const a,b:single) c:single;inline; begin c:= a-b * Int(a/b); if SameValue(abs(c),abs(b)) then c:=0.0; end; { procedure testextended; var a,b,c : extended; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; procedure testdouble; var a,b,c : double; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; } procedure testsingle; var a,b,c : single; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; begin testsingle; // testdouble; // testextended; {\$else FPUNONE} begin {\$endif FPUNONE} end. ``` tw33167.pp (1,684 bytes)

#### Activities

 2018-02-13 19:51 reporter   ~0106373 Last edited: 2018-02-13 20:00View 2 revisions Calling Floor() instead of Int() yields the same signs as Wolfram Alpha, and the abs() can be skipped as well, i.e:   operator mod(const a,b:float) c:float; inline;   begin     c:= a-b * Floor(a/b);     if SameValue(c, b) then c := 0.0;   end; 2018-02-13 22:00 administrator   ~0106374 Integer case: 72 mod -11=6 -72 mod 11=-6 So I think we shouldn't follow Wofram Alpha. 2018-02-13 22:38 reporter   ~0106376 I checked with Mathematica, Wolfram rounds towards negative infinity. That is why there is a difference when the result is negative: In[13]:= 72/-11 // N Out[13]= -6.54545 In[14]:= Quotient[72, -11] Out[14]= -7 IMHO, we should follow the currently active FPU rounding mode. 2018-02-13 23:19 administrator   ~0106377 How does the rounding mode influence the sign? 2018-02-14 08:32 reporter   ~0106383 It would, if the definition was   c:= a-b * Round(a/b); But I have since found out from 0030744 that the current Trunc/Round-to-zero is libc-compatible, so that's better. 2018-02-14 08:36 reporter   ~0106384 Last edited: 2018-02-14 08:38View 2 revisions I have written that code and tested against GNU C and Java indeed. 2018-02-14 09:07 reporter   ~0106385 Please, in this discussion about rounding, don't forget the main concern of this report: the error that "mod" may be off by one divisor if the expected result is 0. 2018-02-24 19:23 administrator   ~0106591 Fixed as proposed. 2018-02-26 01:14 reporter   ~0106635 Works only when the operands are "float" or "extended". "double" still yields the old error. -- Run the test project given under "Steps to reproduce" (remove the local operator mod). 2018-02-26 03:40 reporter   ~0106636 Last edited: 2018-02-26 03:43View 2 revisions Every calculator I have here tells me that Wolfram Alpha is wrong. I used the following and got all the same answers with my fx7400g, Windows calculor, and other ti's etc.  I suppose I could dig out my HP, Too. Function myMod(x,y:Currency):Currency; Begin  Result := X-Y*Int(X/Y); end;  I don't know what to say, are we so sure Wolfram Alpha is really that correct? 2018-02-26 07:46 reporter   ~0106637 Last edited: 2018-02-26 08:00View 4 revisions Note to wp: Float maps to the biggest available float size, i.e the highest available precision: - double on x86_64 - extended on i386-32 - double on ARM and AARCH64 This is the way the math unit handles float. I will run tests against GNU C and Java again to see if there is a incompatibility introduced... I am not sure about the fmod function overloads that were already in math before my patch. Compatibility with the behavior of GNU C(++) Delphi and Java is more important to me than Wolfram/Mathematica. On first sight the solution   if SameValue(abs(c), abs(b)) then c := 0.0; should work for all floats except when you pass mixed precision values. In the latter case there need to be overloaded operators for samevalue to work. That's because of the rounding of the least significant bits for different sizes. Java is always double, C++ either single or double. In the case of passing a single and a double the same value test will fail because the single is expanded to double with zero's, not e.g nine's 2018-02-26 08:06 reporter   ~0106638 Last edited: 2018-02-26 08:16View 2 revisions Note that in my original patch the results are correct for the *mathematically* correct fact that float mods are always positive, so I used *floor*. However that is not *computationally* correct since in most languages negative values for modulo are allowed. So maybe store sign and use floor() * sign(). I will spend some time on it today an put it on the mailing list to prevent bloat here. Since left-most determines sign for modulo (that is an existing rule) what I am referring to is:   operator mod(const a,b:float) c:float;inline;   begin     c:= a-b * floor(a/b) * Sign(a);   end; 2018-02-26 10:28 reporter   ~0106642 FMod-function has same bug. 2018-03-14 21:43 administrator   ~0107114 I tested the similar approach for double and single and at least for single it fails (see my attachement). So I think it is better just to limit the mod operator to float or even remove it again altogether, this seems to be a very slippery slope. @delfion: FMod shall be probably delphi compatible, so no need to change it. 2018-03-14 21:44 administrator tw33167.pp (1,684 bytes) ```{\$ifndef FPUNONE} uses math; operator mod(const a,b:single) c:single;inline; begin c:= a-b * Int(a/b); if SameValue(abs(c),abs(b)) then c:=0.0; end; { procedure testextended; var a,b,c : extended; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; procedure testdouble; var a,b,c : double; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; } procedure testsingle; var a,b,c : single; begin a:=7.7; b:=1.1; c:=a mod b; if not(SameValue(c,0.0)) then begin writeln(c); halt(1); end; a:=7.2; b:=-1.1; c:=a mod b; if not(SameValue(c,0.6)) then begin writeln(c); halt(1); end; a:=-7.2; b:=1.1; c:=a mod b; if not(SameValue(c,-0.6)) then begin writeln(c); halt(1); end; writeln('ok'); end; begin testsingle; // testdouble; // testextended; {\$else FPUNONE} begin {\$endif FPUNONE} end. ``` tw33167.pp (1,684 bytes) 2018-03-14 21:50 reporter   ~0107115 Last edited: 2018-03-14 21:53View 2 revisions Florian, If you think so, remove it. It is your child. I merely tried (and tested) a common feature for common languages. And it works as advertised: try the java or c++ code. It is compatible. I already explained that mod isn't mod... But then you must also remove the fmod's which behave exactly the same. 2018-03-14 21:52 reporter   ~0107116 @delphion: it is not a bug. It is a matter of interpretation consistent with other languages. 2018-03-14 21:56 reporter   ~0107117 Since it is compatible with major languages I suggest to leave it in. Both mod and fmod. 2018-03-15 09:07 reporter   ~0107123 > it is not a bug. It is a matter of interpretation consistent with other languages. What are you talking about, the rounding error, or the sign inconsistency? I don't care so much about the sign, but rounding must be correct for all data types, i.e. 7.7 mod 1.1 must be 0 independently of whether the input values are provided a single, double or extended. Otherwise this operator is useless and must be removed. 2019-12-06 10:46 reporter   ~0119656 Any news on this? I see that current trunk still exhibits the issue reported. 2019-12-07 03:04 reporter   ~0119671 Function MyFMod(A,B:Double):Double;  Var    cA,cB:Currency;  Begin    IF A < 0 Then cA := -A else cA := A;    If B < 0 Then cB := -B else cB := B;    Result := cA-Cb*Int(cA/cB);    If B < 0 Then Result := -Result;  End; This yields Zero for your 7.7, 1.1 test/ procedure TForm1.Button1Click(Sender: TObject); begin   Caption := MyFMod(7.7,-1.1).ToString; end; Testing with other values posted here, they are same except for the two 0.5, with my results its actually 0.5999999999 and its displayed as 0.6  The trick is to move this all into integers and perform the math that way.. I used Currency because of that, I suppose two side by side Int64 per digit can be processed to yield a long fraction it needed. 2019-12-09 13:05 reporter   ~0119707 Last edited: 2019-12-09 13:09View 3 revisions To me the compatibility with the GCC compiler suite and Java is essential. I wrote that - original - code based on math.fmod(), nothing more (and suggested sign() as a possible solution). Any other solution would be to adapt fmod too and file reports against the GCC compiler suite. I would keep mod and fmod as-is. Modulo is already slow enough and does not need more comparisons.

#### Issue History

2018-02-13 19:30 wp New Issue
2018-02-13 19:51 wp Note Added: 0106373
2018-02-13 20:00 wp Note Edited: 0106373 View Revisions
2018-02-13 22:00 Florian Note Added: 0106374
2018-02-13 22:38 Martok Note Added: 0106376
2018-02-13 23:19 Florian Note Added: 0106377
2018-02-14 08:32 Martok Note Added: 0106383
2018-02-14 08:38 Thaddy de Koning Note Edited: 0106384 View Revisions
2018-02-14 09:07 wp Note Added: 0106385
2018-02-24 19:23 Florian Fixed in Revision => 38332
2018-02-24 19:23 Florian Note Added: 0106591
2018-02-24 19:23 Florian Status new => resolved
2018-02-24 19:23 Florian Fixed in Version => 3.1.1
2018-02-24 19:23 Florian Resolution open => fixed
2018-02-24 19:23 Florian Assigned To => Florian
2018-02-26 01:14 wp Note Added: 0106635
2018-02-26 01:14 wp Status resolved => feedback
2018-02-26 01:14 wp Resolution fixed => reopened
2018-02-26 03:40 jamie philbrook Note Added: 0106636
2018-02-26 03:43 jamie philbrook Note Edited: 0106636 View Revisions
2018-02-26 07:50 Thaddy de Koning Note Edited: 0106637 View Revisions
2018-02-26 07:55 Thaddy de Koning Note Edited: 0106637 View Revisions
2018-02-26 08:00 Thaddy de Koning Note Edited: 0106637 View Revisions
2018-02-26 08:16 Thaddy de Koning Note Edited: 0106638 View Revisions
2018-02-26 10:28 delfion Note Added: 0106642
2018-03-14 21:43 Florian Note Added: 0107114
2018-03-14 21:44 Florian File Added: tw33167.pp
2018-03-14 21:44 Florian Assigned To Florian =>