Calculation errors of floating point "mod" operation
Original Reporter info from Mantis: wp @wpam
-
Reporter name:
Original Reporter info from Mantis: wp @wpam
- Reporter name:
Description:
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;
Steps to reproduce:
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.
Additional information:
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
Mantis conversion info:
- Mantis ID: 33167
- Version: 3.1.1
- Fixed in version: 3.1.1
- Fixed in revision: 38332 (#8a2cf56d)