View Issue Details

IDProjectCategoryView StatusLast Update
0033167FPCRTLpublic2019-12-09 13:09
ReporterwpAssigned To 
PrioritynormalSeverityminorReproducibilityhave not tried
Status newResolutionreopened 
Product Version3.1.1Product Build 
Target VersionFixed in Version3.1.1 
Summary0033167: Calculation errors of floating point "mod" operation
DescriptionThe 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;

Steps To ReproduceRun 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.
Additional InformationI 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
TagsNo tags attached.
Fixed in Revision38332
FPCOldBugId
FPCTarget
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

wp

2018-02-13 19:51

reporter   ~0106373

Last edited: 2018-02-13 20:00

View 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;

Florian

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.

Martok

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.

Florian

2018-02-13 23:19

administrator   ~0106377

How does the rounding mode influence the sign?

Martok

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.

Thaddy de Koning

2018-02-14 08:36

reporter   ~0106384

Last edited: 2018-02-14 08:38

View 2 revisions

I have written that code and tested against GNU C and Java indeed.

wp

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.

Florian

2018-02-24 19:23

administrator   ~0106591

Fixed as proposed.

wp

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).

jamie philbrook

2018-02-26 03:40

reporter   ~0106636

Last edited: 2018-02-26 03:43

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

Thaddy de Koning

2018-02-26 07:46

reporter   ~0106637

Last edited: 2018-02-26 08:00

View 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

Thaddy de Koning

2018-02-26 08:06

reporter   ~0106638

Last edited: 2018-02-26 08:16

View 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;

delfion

2018-02-26 10:28

reporter   ~0106642

FMod-function has same bug.

Florian

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.

Florian

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)

Thaddy de Koning

2018-03-14 21:50

reporter   ~0107115

Last edited: 2018-03-14 21:53

View 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.

Thaddy de Koning

2018-03-14 21:52

reporter   ~0107116

@delphion: it is not a bug. It is a matter of interpretation consistent with other languages.

Thaddy de Koning

2018-03-14 21:56

reporter   ~0107117

Since it is compatible with major languages I suggest to leave it in. Both mod and fmod.

wp

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.

wp

2019-12-06 10:46

reporter   ~0119656

Any news on this? I see that current trunk still exhibits the issue reported.

jamie philbrook

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.

Thaddy de Koning

2019-12-09 13:05

reporter   ~0119707

Last edited: 2019-12-09 13:09

View 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

Date Modified Username Field Change
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:36 Thaddy de Koning Note Added: 0106384
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:46 Thaddy de Koning Note Added: 0106637
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:06 Thaddy de Koning Note Added: 0106638
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 =>
2018-03-14 21:50 Thaddy de Koning Note Added: 0107115
2018-03-14 21:52 Thaddy de Koning Note Added: 0107116
2018-03-14 21:53 Thaddy de Koning Note Edited: 0107115 View Revisions
2018-03-14 21:56 Thaddy de Koning Note Added: 0107117
2018-03-15 09:07 wp Note Added: 0107123
2018-03-15 09:07 wp Status feedback => new
2019-12-06 10:46 wp Note Added: 0119656
2019-12-07 03:04 jamie philbrook Note Added: 0119671
2019-12-09 13:05 Thaddy de Koning Note Added: 0119707
2019-12-09 13:08 Thaddy de Koning Note Edited: 0119707 View Revisions
2019-12-09 13:09 Thaddy de Koning Note Edited: 0119707 View Revisions