4

I'm running a program that evaluates equalities after I substitute all symbols for numbers. The equalities on floats are giving me trouble. My program is reaching a point where something like the following is being evaluated: 2.304793200735844`*^-8 == 2.30479302310016`*^-8 . This gives False, despite the difference being 1.77636*10^-15, which is below the Chop threshold of 10^-10. Running Chop on both numbers obviously doesn't fix the problem, since neither number is individually below the threshold.

Chop[2.304793200735844`*^-8 - 2.30479302310016`*^-8 ]==0 gives the correct result, but I can't easily make this change as not all comparisons are between floats or even lists of floats (Chop[a-b] == 0 or {0,0,0,...} is not always a valid way to test equality between a and b).

I want to avoid overcomplicating my program: I don't want to do a ReplaceAll on a_==b_ /; (a and b are numbers) -> Chop[a-b]==0 as this is very complicated, requires something like Hold, and doesn't resolve the same exact issues which arise in the case of float inequality comparisons. There's gotta be a better way! If there's not, what's the best solution?

EDIT: My solution was to modify Equal (==) for float comparison, and to update List equality to use this new float equality. It's extremely simple, but I think it might slightly slow things down:

Unprotect[Equal]; (* this lets us edit Equal (==) *)
Equal[a_Real,b_Real]:=Chop[a-b]==0  (* change comparison for reals *)
Equal[a_List?RealQ,b_List?RealQ]:=
  Dimensions@a==Dimensions@b&&AllTrue[MapThread[#1==#2&,{a,b}],TrueQ] (* update comparison for real lists to use the new real comparison *)
Protect[Equal];

(* True iff r is a real or any nested List containing only reals *)
RealQ[r_]:=DeveloperMachineRealQ[r]||DeveloperRealQ[r]||(MatchQ[r,_List]&&AllTrue[r,RealQ]);
Lucas Mumbo
  • 578
  • 2
  • 11
  • 2
    From the documentation for Equal, "Approximate numbers with machine precision or higher are considered equal if they differ in at most their last seven binary digits (roughly their last two decimal digits)." If you want a different tolerance, use Abs[x1 - x2] < tolerance – Bob Hanlon Mar 08 '22 at 04:05
  • 3
    @Bob's rec is a standard floating-point programming approach for absolute error. For relative error, use something like Abs[x1 - x2]/Max@Abs[{x1,x2}] < relTolerance. This last may be programmed in Mathematica by setting Block[{Internal`$EqualTolerance = Log10[2*relTolerance/$MachineEpsilon]}, x1 == x2]. Or Block[{Internal`$EqualTolerance = Log10[2*relTolerance/$MachineEpsilon]}, code] in which each == between real numbers will be compared with relTolerance relative tolerance. For your example numbers, relTolerance should be about 7.70723 * 10^-8 or greater. – Michael E2 Mar 08 '22 at 05:03
  • 3
    See https://mathematica.stackexchange.com/a/35094/4999 and related Q&A – Michael E2 Mar 08 '22 at 05:19

1 Answers1

5

I put this in a comment, but then noticed the OP had answered in the question. The relative error of the example numbers is 7.707229332944118`*^-8, which is pretty high, imo. Nonetheless, there is a way to handle relative error: Internal`$EqualTolerance. To set it to how many bits to ignore -- it's really Log10 of 2^bits -- compare it the relative tolerance to one ulp = $MachineEpsilon/2 (roughly).

Block[{Internal`$EqualTolerance = 
   Log10[2*7.707229332944118`*^-8/$MachineEpsilon]},
 2.304793200735844`*^-8 == 2.30479302310016`*^-8]

(*  True  *)


Block[{Internal$EqualTolerance = Log10[1.999999*7.707229332944118^-8/$MachineEpsilon]},
 2.304793200735844`^-8 == 2.30479302310016`*^-8]


(*  False  *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Great! So the answer is Block[{Internal`$EqualTolerance=Log10[2^(number of ignored bits)/$MachineEpsilon]}, ...]? – Lucas Mumbo Mar 08 '22 at 06:07
  • 2
    @AndrewP. Rather, Block[{Internal`$EqualTolerance=Log10[2^(number of ignored bits)]}, ...]. I coded Block[{Internal`$EqualTolerance=Log10[2*relTolerance/$MachineEpsilon]}, ...]. The $EqualTolerance (= number of ignored digits) is set equal to the difference between the logs of the relative tolerance and one ulp. If you know how many bits or digits to ignore, you don't need to bother with $MachineEpsilon. (And if you do happen to use arbitrary-precision numbers then the epsilon depends on the Precision[] of the numbers.) – Michael E2 Mar 08 '22 at 06:45
  • Thanks so much - this is really helpful. Is there a place to find additional definitions of this stuff? I can't find any on the Mathematica website. – Lucas Mumbo Mar 08 '22 at 06:59