Jump to content

Sky Slate Blueberry Blackcurrant Watermelon Strawberry Orange Banana Apple Emerald Chocolate
Photo

Floating point comparison


  • Please log in to reply
No replies to this topic
Laszlo
  • Moderators
  • 4713 posts
  • Last active: Mar 31 2012 03:17 AM
  • Joined: 14 Feb 2005
There are several posts in the Forum concerning problems resulting from exact comparisons of floating point numbers. Here is a function fcmp(x,y, tolerance), which ignores small differences.

Floating point numbers are often computed with rounding errors, when equality is rarely exact. One could check if the numbers differ by a small error, which is cumbersome, because this tolerance should depend on the range of the numbers. Relative errors could also be used, but they get complicated around 0 or when x and y are far from each other.

A simpler approach is to look at the binary representation of floating point numbers. At proper rounding only the least significant few bits are inaccurate, so it is natural to give the tolerance for equality as the count in these bits. Interestingly, this works over different exponents, too. The binary representation of 1.99999999999999 and 2.00000000000001 are very close: Adding a small number to the first mantissa causes an overflow to the exponent, so we get the second number with its one larger exponent. It is described in more details by Bruce Dawson, for single precision floats.

For this floating point comparison we need 64 bit integer arithmetic, supported by AHK, except the little complication that the necessary constants overflow, but the resulting code will be slow. Therefore, the code below is written in C, and imported to AHK as a machine code function.
int fcmp(double A, double B, int maxUlps) { // 0 <= maxUlps < 2**53, tolerance in LS bits
    Int64 intDiff, aInt = *(Int64*)&A, bInt = *(Int64*)&B;
    if (aInt < 0) aInt = 0x8000000000000000 - aInt;
    if (bInt < 0) bInt = 0x8000000000000000 - bInt;
    intDiff = aInt - bInt;
    if (-maxUlps <= intDiff && intDiff <= maxUlps) return 0;
    return aInt < bInt ? -1 : 1;
}
The machine code is generated in the usual way, with the MCode function, and the fcmp function is called by a dllcall:
MCode(fcmp, "558bec83e4f883ec148b5510538b5d0c85db568b7508578b7d148974241889542410b9000000807f"
. "137c0485f6730d33c02bc68bd91b5d0c89442418837d14007f137c0485d2730d33c02bc28bf91b7d1489442410"
. "8b7424182b7424108b45188bcb1bcff7d8993bd17f187c043bc677128b4518993bca7f0a7c043bf0770433c0eb"
. "183bdf7f117c0a8b44241039442418730583c8ffeb0333c0405f5e5b8be55dc3")

Setformat Float, 0.18e
x := 1.0, y := 1.00000000000001 ; example
MsgBox % DllCall(&fcmp, "double",x, "double",y, "Int",50, "CDECL Int")

MCode(ByRef code, hex) { ; allocate memory and write Machine Code there
   VarSetCapacity(code,StrLen(hex)//2)
   Loop % StrLen(hex)//2
      NumPut("0x" . SubStr(hex,2*A_Index-1,2), code, A_Index-1, "Char")
}

If x = y within the tolerance tol in the least significant bits of the double precision floating point arguments, fcmp(x,y,tol) is 0. Otherwise, it is –1, or +1 if x < y or x > y, respectively. In the example above x and y are equal within a tolerance of 50, that is, only their least significant 6 bits are different of the 64. Don’t specify negative or too large tolerance, like 9007199254740992, because it causes internal overflow, and the result will not be accurate. A tolerance of 0 checks for exact equality.

Some examples:
x := 1/3, y := (1-x)/2 are equal, within a tolerance of 1 (not equal with 0 tolerance)
x := sqrt(2), y := 2/x are equal, within a tolerance of 1 (not equal with 0 tolerance)
x := tan(atan(99)), y := 99: x < y with a tolerance of 1..12, but x = y within a tolerance of 13.