That's a neat trick, "cancelling" out catastrophic cancellation. Of course it doesn't work anymore once x is smaller than machine epsilon, whereas expm1(x) / x will continue to work.<p>In general herbie is pretty good at suggesting the right functions / rearrangements. For this example, it finds expm1: <a href="https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e00f46723916.f54ecc7459e2826747d647a0d98ec9e41c611ae9/graph.html" rel="nofollow">https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e...</a>
I guess "floats are weird" is a catchier title than "numerical computing is an acquired skill based in part on understanding the various consequences of value representation density being inversely proportional to the absolute value".
Ages ago, when I was in college, a "numerical methods" course was more or less obligatory... has that changed?<p>Anyway, toying with a slide rule is a fun way of making these things very obvious, because you are essentially working with something like 3-digit floating point. When you find yourself subtracting 1.234 from 1.235, but you know that both the "4" and the "5" are really fuzzy around the edges, it becomes obvious that the resulting 0.001 does not really mean much. Similarly, when you take 1.234e+6 (again with a really fuzzy "4") and add 2, you realize that 1.234002e+6 is not really a sane answer.
For this specific case, use Python's expm1, see <a href="https://docs.python.org/3/library/math.html#math.expm1" rel="nofollow">https://docs.python.org/3/library/math.html#math.expm1</a> and history of expm1 at <a href="https://en.wikipedia.org/wiki/Exponential_function#expm1" rel="nofollow">https://en.wikipedia.org/wiki/Exponential_function#expm1</a><p><pre><code> def f(x):
return math.expm1(x)/x
</code></pre>
The expm1(x) means "exp(x) minus 1".<p><pre><code> >>> f(1e-15)
1.0000000000000007
</code></pre>
The technique described gives<p><pre><code> 1.0000000000000004
</code></pre>
which is 1 step smaller than the value computed via expm1():<p><pre><code> >>> math.nextafter(f(1e-15), 0)
1.0000000000000004
</code></pre>
Both are within 1 ulp of the more precise value from WolframAlpha:<p><pre><code> 1.0000000000000005000000000000001666666666666667083333333333333416...</code></pre>
Re money systems: I have found that problems arise no matter whether floats are used or something else: The system can be too exact.<p>My current approach is that I try to guess (or ask) how a customer understands or checks a statement/invoice. The customer usually takes a calculator and enters the rounded numbers they see. I then try to do all programmatic calculations exactly the same way: Operation - round - operation - round etc, with rounding-steps at exactly the same places where these values will be visible to the customer, and with the same number of decimal places (can vary). Sometimes that is not possible of course, which is usually solved with a final row called 'rounding errors in the customer's favor').<p>For my (smaller) systems I find floats sufficient and easier to work with compared to the alternatives, which tend to make code very verbose and annoying to read.
That's called the stability of the algorithm. The first algorithm is not stable whereas the second with the logarithm is.<p>Another example:
compare the two algorithms, that are actually pure-mathematically equal:<p><pre><code> def fc(x):
return sqrt(x+1) - sqrt(x)
def fd(x):
return 1/(sqrt(x+1) + sqrt(x))
</code></pre>
For large values of x (1e16 for example) plot the results and see the difference.<p><pre><code> xs = np.linspace(10\*14,10\*16,10000)
plt.figure(figsize=(8, 6), dpi=120)
plt.plot(xs,[fc(x) for x in xs])
plt.plot(xs,[fd(x) for x in xs])</code></pre>
i think it is our educational system that is weird.<p>we teach students from arithmetic to differential equations assuming they are working with Real Numbers and never really talking about what that means or what a Field is. Then we present computation as something that can easily be done with Real numbers. Which is a lie because we don't even know if Pi plus E is rational or not.<p>Then we present computers as tools of computation, handwaving away the fact that Computer floating point math is not using real numbers, and is not a Field nor anything even remotely like a Field. Floats are not even closed under basic arithmetic operations. And god help someone who learned Geometry on graph paper then tries to apply it to a computer with Floats because Floats cannot even represent a uniform grid in space, they vary in resolution by definition depending on distance from origin.<p>So in the end by trying to gloss over the details of how Real Numbers work, and how Computers work, we have actually misrepresented both Real Numbers and Computers to generations of people, for going on 80 years now. And we keep having to have these articles, or "your calculator is wrong" videos go viral every few months/years because people were not taught basic facts about Reals and Computers.
This is a whole scientific field, called numerical analysis: <a href="https://en.wikipedia.org/wiki/Numerical_analysis" rel="nofollow">https://en.wikipedia.org/wiki/Numerical_analysis</a>
It’s not so much that floats are weird, but that they are not reals. Most of the problems happen due to confusion about this fact.<p>Reals are much weirder/more interesting than floats, imo.
Then what's the best way to handle these cases? Are there any set of rules we should use while implementing the mathematical equations dealing with limits in floating points.
A number of years back I was programming in ADA, a language I generally liked, and had a case where the compiled version of a floating point constant differed from the value generated when the same string was entered at run time. Evidently, the compiler and the run time used different routines to convert strings to floats.<p>The difference wasn't much, but it was there.
Nothing weird about it. It should be obvious that subtracting two floats that are very close to each other results in a loss of numerical precision:<p>1.000000003456e0 - 1.000000002345e0 = 0.000000001111e0 = 1.111numericalnoise e-9<p>It's exactly the same issue here. `math.exp(1e-15)` is `1.000000000000001`. If you subtract 1, you get 1 significant digit and numerical noise.
I decided years ago that the next time I hear someone suggesting we use floats / doubles to represent money amounts, I am going to punch them in the face.
I’ve had weird float bugs that I didn’t have time for and fixed by operating on them as strings. Nowadays I’ve found pretty good libraries in common languages designed to handle the weird edges.
Hotter take: Floats are bad and shouldn't be used.<p>It's clear that even experienced programmers do NOT understand how they work, how they don't work, and the pile of edge cases.<p>Use something saner, like binary coded decimal, larger types, and use floats as a last and very approximate resort.