# Dr. StrangeNumbers or: How I Learned to Stop Worrying and Love Floating Point Arithmetic

The following is the program.

test[t_, dt_] :=
Module[{}, For[ti = dt, ti <= t, ti = ti + dt, Print[ti];];
Print[MemberQ[{0.01, 0.02}, ti]]; Return[0];]
test[0.01, 0.001]

(* 0.001
...
False
0 *)


Obviously the result is wrong. Copy the above result, you will be surprised to see:

(* ...
0.009000000000000001
0.010000000000000002
... *)


Why is this so?

-
This is not a Mathematica question - it is a numerical analysis question: en.wikipedia.org/wiki/Round-off_error –  Mark McClure May 15 '12 at 12:41
It's not the For[] loop, really. Try 0.008 + 0.001 // InputForm. Welcome to the fun of floating-point arithmetic! –  Ｊ. Ｍ. May 15 '12 at 12:41
Good,Thank you! –  howard May 15 '12 at 13:16
What's the point of Module[{}, ...] with empty variable list? As far as I understand, it should not have any effect at all. –  celtschk May 16 '12 at 13:22
@celtschk I've used it a few times as a cool parentheses, reserved for future use... Then I stopped doing that –  Rojo May 16 '12 at 18:34

While it is a bad idea to compare floating point numbers, in this case I think something simpler is going on: you have an off-by-1 problem in your loop. See what your code does:

test[t_, dt_] := Module[{},
For[ti = dt, ti <= t, ti = ti + dt, Print[ti];];
Print[MemberQ[{0.01, 0.02}, ti]];
Return[0];
]


then, after

test[0.01, .001];


look at ti:

ti
(*0.011*)


If you use < rather than <= in the condition, the final value will indeed be what you expect it to be:

test2[t_, dt_] := Module[{},
For[ti = dt, ti < t, ti = ti + dt, Print[ti];];
Print[MemberQ[{0.01, 0.02}, ti]];
Return[0];
]

test2[0.01, .001]
ti


prints a list up to 0.009 and returns True.

In general, though, don't do things like that with floating-point numbers. Even if it works here, it will eventually fail. Observe:

(1 + $MaxMachineNumber) ==$MaxMachineNumber
(*True*)


or, as JM points out,

1 + $MachineEpsilon == 1 (*True*)  - The more traditional example is 1 +$MachineEpsilon == 1... ;) –  Ｊ. Ｍ. May 15 '12 at 12:59
@J.M. Tradition is for wimps :) –  acl May 15 '12 at 13:01
Bear in mind Equal applies an extra tolerance in Mathematica. The proper comparison is Block[{Internal$EqualTolerance = -Infinity}, 1 == 1 +$MachineEpsilon] (= False) vs. Block[{Internal$EqualTolerance = -Infinity}, 1 == 1 +$MachineEpsilon/2] (= True). –  Oleksandr R. May 15 '12 at 13:55
@OleksandrR. I didn't know that, thanks –  acl May 15 '12 at 14:05

This seems to work fine and is more in the spirit of the Mathematica way and perhaps illustrates why a Functional approach tends to lead to fewer bugs than the procedural approach:

test[t_, dt_] := {#, MemberQ[{0.01, 0.02}, #]} & /@ Range[dt, t, dt];


Which when run,

test[0.01, 0.001] // TableForm


gives:

-
Replacing TableForm with InputForm reveals that one does not totally avoid floating-point error... –  Ｊ. Ｍ. May 15 '12 at 13:55
@J.M. that's interesting, it makes me curious as to why that happens, the increment doesn't seem to be anywhere near the limit of precision of of a floating point number. Is this a general effect with FP numbers or more specific to MMA? –  image_doctor May 15 '12 at 15:26
it's general, see e.g. floating-point-gui.de –  acl May 15 '12 at 15:45
@acl, Thanks, surprising how in a great number of years that's never knowingly tripped me up, but here because MemberQ is being used as a criterion, rather than a < or > based test, it become obvious. –  image_doctor May 15 '12 at 16:00