Floating Number Woes

The following snippet evaluates to 0 in Mathematica:

fp = 1/2 (fc + Sqrt[fc^2 - 2*fz^2]);
fm = 1/2*(fc - Sqrt[fc^2 - 2*fz^2]);
fp^2 + fm^2 + fz^2 - fc^2 // Simplify

But when you run

fc = 30000000.;
fz = 4000000.;

and evaluate the first bit of code again - the result will be non-zero (on my machine: 0.125).

It took me embarrassingly long to figure out that this was a floating point problem. I have computed this sort of code for the last 5 years, and somehow I must have gotten lucky, because I always computet the last line slightly different:

Sqrt[fp^2 + fm^2 + fz^2] - fc

The way that floating point numbers are stored makes the second version much, much more benign - but for me, such floating points woes are not intuitive yet.

I am somewhat sad that it's 2015 and these problems still sneak up on us. Mathematica should have warned be about this. We've had "modern" interval arithmetic since the 1960s - machines have the power to warn when they might be outputting garbage! (Edit: See bottom of page for why and when Mathematica warns me)

Indeed, you can easily use Intervals in Mathematica (and in Matlab, it's a popular package) to show that the result "0.125" is such garbage:

In Mathematica, you can define arbitrary intervals using the Interval function. When this function is given a single argument, it gives you the smallest Interval that can be represented with standard machine precision. Therefore, define

fc = Interval[30000000.];
fz = Interval[4000000.];

and rerun the first snippet - the result is Interval[{-1.75, 1.875}]. This interval is bigger than the actual error, but this is a deliberate trade-off to guarantee that the true result is within the interval, without having to use brutally complicated math for the proof. You can get a smaller interval by using arbitrary precision numbers, for instance with 200 digits:

fc = Interval[30000000.`200];
fz = Interval[4000000.`200];

Oh well, now I know to be more careful. And instead of complaining I should probably start setting a good example by using Interval arithmetic from the start. It's actually quite simple!

Edit: I just learned why Mathematica did not warn about the garbage that it returned. When using machine-precision numbers (ordinary floating numbers like 1.2), then it does not do any fancy interval arithmetic or precision control with them. This makes the calculations reasonably fast. When using arbitrary precision numbers (like 1.2`20), then Mathematica does reason about the precision and warns about garbage. Normal machine precision is a precision of about 16 decimal digits. When specifying all numbers as 4000000`16, and so on, then Mathematica warns me. This is a speed/convenience trade-off - still, I wish Mathematica was smart enough to see that my calculation costs no time whatsoever and simply rerun it with arbitrary precision numbers to give me a warning...

Bottom line: Instead of using Interval, one can also use arbitrary precision numbers and trust that (hopefully) Mathematica gives more sensible output.