12 ms·

Floats Are Weird

154 points by soopurman 3y ago

mjcohen 3y ago
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. The difference wasn't much, but it was there.
user3939382 3y ago
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.
- danhau 3y ago
  What libraries can you recommend?
  - user3939382 3y ago
    For JS for example https://github.com/MikeMcl/decimal.js https://github.com/MikeMcl/decimal.js The trick is 90% just realizing these types of libraries exist so if and when you run into these requirements you’re not reinventing the wheel on some of these issues.
098799 3y ago
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: 1.000000003456e0 - 1.000000002345e0 = 0.000000001111e0 = 1.111numericalnoise e-9 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.
- kolinko 3y ago
  That's not what the post says if I understand correctly - the post explains why in certain situations the "noise" disappears, and in other cases it doesn't. See comparison between f and g functions.
  - 098799 3y ago
    I see! yes, the magic is you can cancel the noise by repeating it twice: ``` In [1]: math.exp(1e-15)-1 Out[1]: 1.1102230246251565e-15 In [2]: math.log(math.exp(1e-15)) Out[2]: 1.110223024625156e-15 ``` risky business though, I imagine it's implementation dependent
    - Sharlin 3y ago
      It’s not (or shouldn’t be), it’s simply a result of math, as the article explains in length.
      - lifthrasiir 3y ago
        Many libm implementations don't have an accurate `log` or `exp` routine, so there does exist a risk. (Of course, it's also true that many of them also special-case `log(x) ~= x - 1` and `exp(x) ~= x + 1` for small enough `x`.)
        planede 3y ago
        The math hinges on that there is the same error for exp(x) at both places. So as long as exp(x) is deterministic then this should be alright.
        adgjlsfhk1 3y ago
        I don't know of any libm that have log or exp sufficiently inaccurate for this to break. Do you?
        3y ago
an1sotropy 3y ago
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".
- stabbles 3y ago
  I think the title is alright. The curious thing is that both `f` and `g` have catastrophic cancellation, so you expect both to be inaccurate, but magically `g` recovers from it. Alternative title: catastrophic cancellation cancels out
  - Someone 3y ago
    > I think the title is alright I disagree (for the title I see now, which is “Floats are weird”). The problem isn’t about floats, but about inexact computation and, as you say, catastrophic cancellation. If, for example, you do all your computations in 7.6 digit base 11 numbers (seven digits before the undecimal point, six after it), you get the same problem.
- topaz0 3y ago
  It's at least nice to see somebody reasoning through where error comes from instead of just throwing up their hands as if the error was inexplicable. I agree that the title is unfortunate though.
the-Black-Dread 3y ago
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.
- Dwedit 3y ago
  For a summation, add the smaller numbers first. Smaller as in 0.0000000000053, not like -5172365126.
  - Dban1 3y ago
    how about small as in ²³⁷
  - zokier 3y ago
    For summation just go with Kahan summation or some improved variant of it https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm
    - lifthrasiir 3y ago
      And if you can't remember that, you can use a pairwise summation like (((a+b)+(c+d))+((e+f)+(g+h)))+... which gives a decent error reduction in return.
- OscarCunningham 3y ago
  Floating point numbers have the highest precision when they're near 0. So if you know a number is going to be near 1, like exp(x) when x is small, then it's better to calculate the offset from 1 rather than the number itself. Programming languages make 'expm1' available for this purpose, which computes exp(x)-1.
- stabbles 3y ago
  Basically - don't subtract almost equal numbers - don't add numbers of vastly different magnitudes
eesmith 3y ago
For this specific case, use Python's expm1, see https://docs.python.org/3/library/math.html#math.expm1 https://docs.python.org/3/library/math.html#math.expm1 and history of expm1 at https://en.wikipedia.org/wiki/Exponential_function#expm1 https://en.wikipedia.org/wiki/Exponential_function#expm1 def f(x): return math.expm1(x)/x The expm1(x) means "exp(x) minus 1". >>> f(1e-15) 1.0000000000000007 The technique described gives 1.0000000000000004 which is 1 step smaller than the value computed via expm1(): >>> math.nextafter(f(1e-15), 0) 1.0000000000000004 Both are within 1 ulp of the more precise value from WolframAlpha: 1.0000000000000005000000000000001666666666666667083333333333333416...
- Sharlin 3y ago
  To be precise, it’s C’s `expm1`. Like other things in math.h, it’s just adopted as is to most popular languages, even including PHP and JS.
  - eesmith 3y ago
    To be precise, it's most likely Kahan's `expm1`, since Wikipedia attributes it to him. It adds that expm1 was first implemented on an HP calculator in the 1970s. It became part of BSD 4.3 due to Kahan and one of his students. The documentation at https://archive.org/details/unix-programmers-reference-manual-4.3-berkeley-software-distribution/page/n161/mode/2up?q=kahan+expm1 https://archive.org/details/unix-programmers-reference-manua... mentions "expm1" was available in HP BASIC, and mentions the Apple C compiler at the time used "exp1", which you can verify from the Apple Numerics Manual Second Edition 1988 at https://archive.org/details/apple-numerics-manual-second-edition-1988/page/n85/mode/2up?q=exp1+compiler https://archive.org/details/apple-numerics-manual-second-edi... That makes it BASIC's expm1 before it was C's. ;) While yes, Python's expm1 uses the platform libm's expm1 through C, Python's math module will avoid using the vendor libm for cases where the vendor implementation is known to have numerical issues. For one such example, on macOS the following: #include <math.h> #include <stdio.h> int main(void) { double x = -1.0; x = nextafter(x, 0.0); // -0.9999999999999999 printf("lgamma(%.19f) -> %.17f\n", x, lgamma(x)); return 0; } reports lgamma(-0.9999999999999998890) -> 36.25168935623486988 while Python, using its own lgamma() implementation: import math x = math.nextafter(-1.0, 0.0) print(f"lgamma({x}) -> {math.lgamma(x)}") reports: lgamma(-0.9999999999999999) -> 36.7368005696771 showing that Python's lgamma is not a pass-through to C's. FWIW, WolframAlpha says the value is 36.73680056... so Python's is more accurate. EDIT: Python's math.expm1 was added in version 3.2 back in February 20th, 2011, and did not depend on C's expm1. Instead, it used Kahan's method directly: /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed to avoid the significant loss of precision that arises from direct evaluation of the expression exp(x) - 1, for x near 0. */ double _Py_expm1(double x) { /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this also works fine for infinities and nans. For smaller x, we can use a method due to Kahan that achieves close to full accuracy. */ if (fabs(x) < 0.7) { double u; u = exp(x); if (u == 1.0) return x; else return (u - 1.0) * x / log(u); } else return exp(x) - 1.0; }
gonzus 3y ago
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.
- SideQuark 3y ago
  Please stop repeating this nonsense. It's purely based on ignorance and only pushes more people into it. This is universally repeated by people that have not written any modern financial software and who don't understand floats. If all you do is add US dollars and pennies, then maybe you can get by with integers. Once you do anything else in modern finance integers puke completely. There's a reason scientific computing doesn't use integers, but prefers floats/double - they are much easier to use to get the best answers per compute. For example, if you need to do anything with interest, which is fundamental, you'll soon find doubles are vastly better than bitsize equivalent integers, no matter what scaling/fixed-point/other tricks you employ. Add in currency differences (Yen to USD is a large multiplier, etc), any longer range calculations (50-100 year loans or flows), aggregation of varying items, derivatives, and on and on. Telling anyone in the field you're going to use integers will get you laughed at - the problem is not floats, the problem is the programmer hasn't taken a single class on numerical analysis and has not enough skill to do financial software. For example, one of the simplest things one needs to do is compute compound interest, say computing mortgage tables, with say principal P (left), annual (or weekly, or daily..) rate R, for N years, periodic payment m, and you want to compute payments and schedule. In reals, this is simple: each cycle you do something like r = (1+R/12) P -= m P = (1+r)*P Then you write out values you need. Trying to write this software with integers, no matter what scaling, fixed-point, shifting, and other tricks you employ, is going to be vastly more complex and error prone than simply using doubles. If you don't think so, pick a bit budget, say 32, or 64, for you base number type, and show me your code. Then I'll show you the naive one with doubles vastly outperforms your code. This continues through all of modern finance. So stop repeating this ignorance that one should use integers for money - that is purely a result of being ignorant about how to write robust numerical code, and only pushes people down a far worse rabbit hole when they hit issues with ints and try to patch them one at a time via ad-hoc hacks. There's a reason scientists use floating point, not ints, to do real numerical work.
- 3y ago
cgadski 3y ago
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. In general herbie is pretty good at suggesting the right functions / rearrangements. For this example, it finds expm1: https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e00f46723916.f54ecc7459e2826747d647a0d98ec9e41c611ae9/graph.html https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e...
- zokier 3y ago
  Gotta love that Herbie also provides just `1` as alternative. "66.8% accurate, 305.0× speedup" But the expression you get if you disable "numerics" ruleset (which includes expm1) is pretty wild also: (if (<= (exp x) 4504410275303423/4503599627370496) (+ 1 (+ (* 1/24 (pow x 3)) (* x (+ 1/2 (* x 1/6))))) (/ (+ (exp x) -1) (log (exp x)))))
- aragilar 3y ago
  That doesn't surprise me, as before floats are involved you need to transform the series into something you can compute (doing this by hand is a good way to understand what's going on underneath). https://developers.redhat.com/blog/2015/01/02/improving-math-performance-in-glibc https://developers.redhat.com/blog/2015/01/02/improving-math... gives a basic intro to how these functions are actually implemented, and we can if you compute exp(x) - 1 near 0 that 1 + (small stuff) - 1 is always going to be a problem (it'd be nice if we had a sufficiently smart compiler that could understand how to inline these things, but then it'd probably do surprising things as well).
Craig_Rodriguez 3y ago
[dead]
pierat 3y ago
Hotter take: Floats are bad and shouldn't be used. It's clear that even experienced programmers do NOT understand how they work, how they don't work, and the pile of edge cases. Use something saner, like binary coded decimal, larger types, and use floats as a last and very approximate resort.
- zokier 3y ago
  As if even experienced programmers understand anything about numerics, regardless of number format used. And all numbers have pitfalls, reals are not representable on computers, not all reals are even computable. Floats at least are well-documented and predictable, and there is 30+ years of literature on how to do error analysis for them.
- aragilar 3y ago
  No, it's clear people don't know how to compute (because they're not taught how to). The vast majority of issues that people encounter with floating point aren't floating point issues, it's "how do I perform this calculation without infinite space and time", and these issues occur whether or not you do it by hand, or use a machine to do it for you. This issue becomes really obvious when you teach early-undergrad science labs, because people conflate the number of decimal points used with the number of significant figures used. This can be seen between the definition of the speed of light (which has in effect infinite precision because we defined what it is) as 299792458 m/s, vs the gravitational constant with 0.000 000 000 066 74 m^3 /kg /s^2 which is generally regarded as the most inaccurate and hardest to measure physical constant (usually G * M can be measured more accurately, so you'd rather use that).
rho4 3y ago
Re money systems: I have found that problems arise no matter whether floats are used or something else: The system can be too exact. 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'). 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.
- fellerts 3y ago
  I thought "money systems" always worked in whole units of some fraction of your currency, e.g. "millicents" or similar, and shunned floating point arithmetic like the plague. > The system can be too exact Indeed, high precision, low accuracy! Fixed-point arithmetic can be both precise and accurate if everyone agrees on the same system. Maybe we don't, and that's the problem?
  - gpderetta 3y ago
    In Europe, after MIFID II, equities are priced with tick sizes (i.e. the minium whole unit) that depend on the magnitude of the price, i.e are de-jure (decimal) floating points that behave weirdly.
  - doubloon 3y ago
    money systems work with whatever tool is convenient to the operations staff which these days is typically Microsoft Excel and Visual Basic, neither of which come with a built in Decimal type. this will typically get transmuted through multiple layers of javascript, xml, back end java, cobol, SQL, python, some domain specific language, and back again.
  - nyrikki 3y ago
    IBM partially has a hold on banking because they have supported decimal floating point. By switching binary floating-point from radix 2 to radix 10 (decimal exponents) they get rid of some of the problems. IEEE has had these Decimal types for a while, but IBM is the only company shipping CPUs with them. Fixed decimal is still used where appropriate, but Decimal floating point is a critical feature for some processes. While personal taxes typically drop the cents, tariffs and other taxes are often small percentages as an example. There are software implementations if you don't need performance. But Decimal64 has a lot more represented values and it is easier to trap on some edge cases. As Decimal floating point types were added into C23 we will see if hardware support follows. Note this is also why HP calculators were more accurate for a given bit size, because the engineering ones used radix 10 floats.
layer8 3y ago
This is a whole scientific field, called numerical analysis: https://en.wikipedia.org/wiki/Numerical_analysis https://en.wikipedia.org/wiki/Numerical_analysis
ska 3y ago
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. Reals are much weirder/more interesting than floats, imo.
- mathgradthrow 3y ago
  the axiomatization of floating point arithmetic is much larger than the axiomatization of real arithmetic, so it is reasonable to describe the former as more complicated.
  - ska 3y ago
    Complexity of the definition isn’t particularly interesting, the properties of the object being defined can be. The reals are stranger than most people realize.
    - JohnKemeny 3y ago
      What's a strange thing about reals?
      - floxy 3y ago
        "Why should I believe in a real number if I can’t calculate it, if I can’t prove what its bits are, and if I can’t even refer to it? And each of these things happens with probability one! The real line from 0 to 1 looks more and more like a Swiss cheese, more and more like a stunningly black high-mountain sky studded with pin-pricks of light." ...see chapter 5 of Meta Math! by Gregory Chaitin: https://arxiv.org/pdf/math/0404335.pdf https://arxiv.org/pdf/math/0404335.pdf
      - adgjlsfhk1 3y ago
        to start, there are an uncountable infinity of them so 100% are not representable. Also when notated typically, there are 2 representations of lots of numbers.
      - ska 3y ago
cratermoon 3y ago
Here's a REALLY weird one: https://latkin.org/blog/2014/11/22/mullers-recurrence-roundoff-gone-wrong/ https://latkin.org/blog/2014/11/22/mullers-recurrence-roundo...
Pinus 3y ago
Ages ago, when I was in college, a "numerical methods" course was more or less obligatory... has that changed? 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.
- bee_rider 3y ago
  I assume it depends on the department. CS students are afraid of non-integer values, engineering students are afraid of codes longer than 200 lines. Just kidding, mostly. I do somewhat wonder if numerical methods have suffered from their success; like everybody uses LAPACK via Numpy but understanding it under the hood, that job always belongs to the next department over wherever you are, haha.
- jmrm 3y ago
  AFAIK in Europe numerical methods is mandatory in most Telco/EE degrees, but not in CS/SWE. The opposite happens with discrete maths btw
  - anthk 3y ago
    Here in Spain if you can't grasp numerical methods in CS you can't even think on earning a degree.
    - rerdavies 3y ago
      Here in Canada, CS is part of the Math Department at the University of Waterloo, but part of the Philosophy Department at University of Toronto. Waterloo graduates take a Numerical Methods course in second year. I fear for the philosophers.
- ArcticLandfall 3y ago
doubloon 3y ago
i think it is our educational system that is weird. 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. 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. 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.
- galdosdi 3y ago
  Great points, let me riff on that though -- what's really "weird" (but totally understandable in the historical context of underpowered computers until recently) is that we think of floats as the default representation of a real at all, when more reasonable modern default choices for all purposes outside of high performance applications would either be a rational composed of bigints, or a bigdecimal, or maybe a generating function for the coefficients of a Taylor series expansion or something like that if you really want to get nuts. In any case you can have as much precision as you can afford, although for anything but rationals it'll be an approximation always, and the question is just how many digits you care for. A lot of choices but the only one that's really absurd is the actual default of floats which are terrible for every purpose except their original raison detre, which was performance in a time when computers were so slow that buying a floating point unit to put in your computer was a normal desirable upgrade for a PC owner (and still crucial for high thoroughput computing, video games, graphics rendering, ML, etc) I do remember being taught the difference between rationals and irrationals in high school. Another aside... The really crazy difference is between the definable and non definable numbers. The first which includes everything from pi to every busy beaver number is countably finite just like the integers, only the nondefinables (numbers whose shortest mathematical definition would not be finite or not exist) give the reals their uncountably infinite nature
architectonic 3y ago
That's called the stability of the algorithm. The first algorithm is not stable whereas the second with the logarithm is. Another example: compare the two algorithms, that are actually pure-mathematically equal: def fc(x): return sqrt(x+1) - sqrt(x) def fd(x): return 1/(sqrt(x+1) + sqrt(x)) For large values of x (1e16 for example) plot the results and see the difference. 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])
echoangle 3y ago
Kind of OT: Why do so many people say “phenomena” when they want to use the singular case?

154 points by soopurman 3y ago

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. The difference wasn't much, but it was there.

user3939382 3y ago

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.

danhau 3y ago
What libraries can you recommend?
- user3939382 3y ago
  For JS for example https://github.com/MikeMcl/decimal.js https://github.com/MikeMcl/decimal.js The trick is 90% just realizing these types of libraries exist so if and when you run into these requirements you’re not reinventing the wheel on some of these issues.

098799 3y ago

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: 1.000000003456e0 - 1.000000002345e0 = 0.000000001111e0 = 1.111numericalnoise e-9 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.

kolinko 3y ago
That's not what the post says if I understand correctly - the post explains why in certain situations the "noise" disappears, and in other cases it doesn't. See comparison between f and g functions.
- 098799 3y ago
  I see! yes, the magic is you can cancel the noise by repeating it twice: ``` In [1]: math.exp(1e-15)-1 Out[1]: 1.1102230246251565e-15 In [2]: math.log(math.exp(1e-15)) Out[2]: 1.110223024625156e-15 ``` risky business though, I imagine it's implementation dependent
  - Sharlin 3y ago
    It’s not (or shouldn’t be), it’s simply a result of math, as the article explains in length.
    - lifthrasiir 3y ago
      Many libm implementations don't have an accurate `log` or `exp` routine, so there does exist a risk. (Of course, it's also true that many of them also special-case `log(x) ~= x - 1` and `exp(x) ~= x + 1` for small enough `x`.)
      - planede 3y ago
        The math hinges on that there is the same error for exp(x) at both places. So as long as exp(x) is deterministic then this should be alright.
      - adgjlsfhk1 3y ago
        I don't know of any libm that have log or exp sufficiently inaccurate for this to break. Do you?
        3y ago

an1sotropy 3y ago

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".

stabbles 3y ago
I think the title is alright. The curious thing is that both `f` and `g` have catastrophic cancellation, so you expect both to be inaccurate, but magically `g` recovers from it. Alternative title: catastrophic cancellation cancels out
- Someone 3y ago
  > I think the title is alright I disagree (for the title I see now, which is “Floats are weird”). The problem isn’t about floats, but about inexact computation and, as you say, catastrophic cancellation. If, for example, you do all your computations in 7.6 digit base 11 numbers (seven digits before the undecimal point, six after it), you get the same problem.
topaz0 3y ago
It's at least nice to see somebody reasoning through where error comes from instead of just throwing up their hands as if the error was inexplicable. I agree that the title is unfortunate though.

the-Black-Dread 3y ago

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.

Dwedit 3y ago
For a summation, add the smaller numbers first. Smaller as in 0.0000000000053, not like -5172365126.
- Dban1 3y ago
  how about small as in ²³⁷
- zokier 3y ago
  For summation just go with Kahan summation or some improved variant of it https://en.wikipedia.org/wiki/Kahan_summation_algorithm https://en.wikipedia.org/wiki/Kahan_summation_algorithm
  - lifthrasiir 3y ago
    And if you can't remember that, you can use a pairwise summation like (((a+b)+(c+d))+((e+f)+(g+h)))+... which gives a decent error reduction in return.
OscarCunningham 3y ago
Floating point numbers have the highest precision when they're near 0. So if you know a number is going to be near 1, like exp(x) when x is small, then it's better to calculate the offset from 1 rather than the number itself. Programming languages make 'expm1' available for this purpose, which computes exp(x)-1.
stabbles 3y ago
Basically - don't subtract almost equal numbers - don't add numbers of vastly different magnitudes

eesmith 3y ago

For this specific case, use Python's expm1, see https://docs.python.org/3/library/math.html#math.expm1 https://docs.python.org/3/library/math.html#math.expm1 and history of expm1 at https://en.wikipedia.org/wiki/Exponential_function#expm1 https://en.wikipedia.org/wiki/Exponential_function#expm1 def f(x): return math.expm1(x)/x The expm1(x) means "exp(x) minus 1". >>> f(1e-15) 1.0000000000000007 The technique described gives 1.0000000000000004 which is 1 step smaller than the value computed via expm1(): >>> math.nextafter(f(1e-15), 0) 1.0000000000000004 Both are within 1 ulp of the more precise value from WolframAlpha: 1.0000000000000005000000000000001666666666666667083333333333333416...

Sharlin 3y ago
To be precise, it’s C’s `expm1`. Like other things in math.h, it’s just adopted as is to most popular languages, even including PHP and JS.
- eesmith 3y ago
  To be precise, it's most likely Kahan's `expm1`, since Wikipedia attributes it to him. It adds that expm1 was first implemented on an HP calculator in the 1970s. It became part of BSD 4.3 due to Kahan and one of his students. The documentation at https://archive.org/details/unix-programmers-reference-manual-4.3-berkeley-software-distribution/page/n161/mode/2up?q=kahan+expm1 https://archive.org/details/unix-programmers-reference-manua... mentions "expm1" was available in HP BASIC, and mentions the Apple C compiler at the time used "exp1", which you can verify from the Apple Numerics Manual Second Edition 1988 at https://archive.org/details/apple-numerics-manual-second-edition-1988/page/n85/mode/2up?q=exp1+compiler https://archive.org/details/apple-numerics-manual-second-edi... That makes it BASIC's expm1 before it was C's. ;) While yes, Python's expm1 uses the platform libm's expm1 through C, Python's math module will avoid using the vendor libm for cases where the vendor implementation is known to have numerical issues. For one such example, on macOS the following: #include <math.h> #include <stdio.h> int main(void) { double x = -1.0; x = nextafter(x, 0.0); // -0.9999999999999999 printf("lgamma(%.19f) -> %.17f\n", x, lgamma(x)); return 0; } reports lgamma(-0.9999999999999998890) -> 36.25168935623486988 while Python, using its own lgamma() implementation: import math x = math.nextafter(-1.0, 0.0) print(f"lgamma({x}) -> {math.lgamma(x)}") reports: lgamma(-0.9999999999999999) -> 36.7368005696771 showing that Python's lgamma is not a pass-through to C's. FWIW, WolframAlpha says the value is 36.73680056... so Python's is more accurate. EDIT: Python's math.expm1 was added in version 3.2 back in February 20th, 2011, and did not depend on C's expm1. Instead, it used Kahan's method directly: /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed to avoid the significant loss of precision that arises from direct evaluation of the expression exp(x) - 1, for x near 0. */ double _Py_expm1(double x) { /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this also works fine for infinities and nans. For smaller x, we can use a method due to Kahan that achieves close to full accuracy. */ if (fabs(x) < 0.7) { double u; u = exp(x); if (u == 1.0) return x; else return (u - 1.0) * x / log(u); } else return exp(x) - 1.0; }

gonzus 3y ago

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.

SideQuark 3y ago
Please stop repeating this nonsense. It's purely based on ignorance and only pushes more people into it. This is universally repeated by people that have not written any modern financial software and who don't understand floats. If all you do is add US dollars and pennies, then maybe you can get by with integers. Once you do anything else in modern finance integers puke completely. There's a reason scientific computing doesn't use integers, but prefers floats/double - they are much easier to use to get the best answers per compute. For example, if you need to do anything with interest, which is fundamental, you'll soon find doubles are vastly better than bitsize equivalent integers, no matter what scaling/fixed-point/other tricks you employ. Add in currency differences (Yen to USD is a large multiplier, etc), any longer range calculations (50-100 year loans or flows), aggregation of varying items, derivatives, and on and on. Telling anyone in the field you're going to use integers will get you laughed at - the problem is not floats, the problem is the programmer hasn't taken a single class on numerical analysis and has not enough skill to do financial software. For example, one of the simplest things one needs to do is compute compound interest, say computing mortgage tables, with say principal P (left), annual (or weekly, or daily..) rate R, for N years, periodic payment m, and you want to compute payments and schedule. In reals, this is simple: each cycle you do something like r = (1+R/12) P -= m P = (1+r)*P Then you write out values you need. Trying to write this software with integers, no matter what scaling, fixed-point, shifting, and other tricks you employ, is going to be vastly more complex and error prone than simply using doubles. If you don't think so, pick a bit budget, say 32, or 64, for you base number type, and show me your code. Then I'll show you the naive one with doubles vastly outperforms your code. This continues through all of modern finance. So stop repeating this ignorance that one should use integers for money - that is purely a result of being ignorant about how to write robust numerical code, and only pushes people down a far worse rabbit hole when they hit issues with ints and try to patch them one at a time via ad-hoc hacks. There's a reason scientists use floating point, not ints, to do real numerical work.
3y ago

cgadski 3y ago

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. In general herbie is pretty good at suggesting the right functions / rearrangements. For this example, it finds expm1: https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e00f46723916.f54ecc7459e2826747d647a0d98ec9e41c611ae9/graph.html https://herbie.uwplse.org/demo/378a0682ec4d3e735c790a58fa97e...

zokier 3y ago
Gotta love that Herbie also provides just `1` as alternative. "66.8% accurate, 305.0× speedup" But the expression you get if you disable "numerics" ruleset (which includes expm1) is pretty wild also: (if (<= (exp x) 4504410275303423/4503599627370496) (+ 1 (+ (* 1/24 (pow x 3)) (* x (+ 1/2 (* x 1/6))))) (/ (+ (exp x) -1) (log (exp x)))))
aragilar 3y ago
That doesn't surprise me, as before floats are involved you need to transform the series into something you can compute (doing this by hand is a good way to understand what's going on underneath). https://developers.redhat.com/blog/2015/01/02/improving-math-performance-in-glibc https://developers.redhat.com/blog/2015/01/02/improving-math... gives a basic intro to how these functions are actually implemented, and we can if you compute exp(x) - 1 near 0 that 1 + (small stuff) - 1 is always going to be a problem (it'd be nice if we had a sufficiently smart compiler that could understand how to inline these things, but then it'd probably do surprising things as well).

Craig_Rodriguez 3y ago

[dead]

pierat 3y ago

Hotter take: Floats are bad and shouldn't be used. It's clear that even experienced programmers do NOT understand how they work, how they don't work, and the pile of edge cases. Use something saner, like binary coded decimal, larger types, and use floats as a last and very approximate resort.

zokier 3y ago
As if even experienced programmers understand anything about numerics, regardless of number format used. And all numbers have pitfalls, reals are not representable on computers, not all reals are even computable. Floats at least are well-documented and predictable, and there is 30+ years of literature on how to do error analysis for them.
aragilar 3y ago
No, it's clear people don't know how to compute (because they're not taught how to). The vast majority of issues that people encounter with floating point aren't floating point issues, it's "how do I perform this calculation without infinite space and time", and these issues occur whether or not you do it by hand, or use a machine to do it for you. This issue becomes really obvious when you teach early-undergrad science labs, because people conflate the number of decimal points used with the number of significant figures used. This can be seen between the definition of the speed of light (which has in effect infinite precision because we defined what it is) as 299792458 m/s, vs the gravitational constant with 0.000 000 000 066 74 m^3 /kg /s^2 which is generally regarded as the most inaccurate and hardest to measure physical constant (usually G * M can be measured more accurately, so you'd rather use that).

rho4 3y ago

Re money systems: I have found that problems arise no matter whether floats are used or something else: The system can be too exact. 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'). 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.

fellerts 3y ago
I thought "money systems" always worked in whole units of some fraction of your currency, e.g. "millicents" or similar, and shunned floating point arithmetic like the plague. > The system can be too exact Indeed, high precision, low accuracy! Fixed-point arithmetic can be both precise and accurate if everyone agrees on the same system. Maybe we don't, and that's the problem?
- gpderetta 3y ago
  In Europe, after MIFID II, equities are priced with tick sizes (i.e. the minium whole unit) that depend on the magnitude of the price, i.e are de-jure (decimal) floating points that behave weirdly.
- doubloon 3y ago
  money systems work with whatever tool is convenient to the operations staff which these days is typically Microsoft Excel and Visual Basic, neither of which come with a built in Decimal type. this will typically get transmuted through multiple layers of javascript, xml, back end java, cobol, SQL, python, some domain specific language, and back again.
- nyrikki 3y ago
  IBM partially has a hold on banking because they have supported decimal floating point. By switching binary floating-point from radix 2 to radix 10 (decimal exponents) they get rid of some of the problems. IEEE has had these Decimal types for a while, but IBM is the only company shipping CPUs with them. Fixed decimal is still used where appropriate, but Decimal floating point is a critical feature for some processes. While personal taxes typically drop the cents, tariffs and other taxes are often small percentages as an example. There are software implementations if you don't need performance. But Decimal64 has a lot more represented values and it is easier to trap on some edge cases. As Decimal floating point types were added into C23 we will see if hardware support follows. Note this is also why HP calculators were more accurate for a given bit size, because the engineering ones used radix 10 floats.

layer8 3y ago

This is a whole scientific field, called numerical analysis: https://en.wikipedia.org/wiki/Numerical_analysis https://en.wikipedia.org/wiki/Numerical_analysis

ska 3y ago

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. Reals are much weirder/more interesting than floats, imo.

mathgradthrow 3y ago
the axiomatization of floating point arithmetic is much larger than the axiomatization of real arithmetic, so it is reasonable to describe the former as more complicated.
- ska 3y ago
  Complexity of the definition isn’t particularly interesting, the properties of the object being defined can be. The reals are stranger than most people realize.
  - JohnKemeny 3y ago
    What's a strange thing about reals?
    - floxy 3y ago
      "Why should I believe in a real number if I can’t calculate it, if I can’t prove what its bits are, and if I can’t even refer to it? And each of these things happens with probability one! The real line from 0 to 1 looks more and more like a Swiss cheese, more and more like a stunningly black high-mountain sky studded with pin-pricks of light." ...see chapter 5 of Meta Math! by Gregory Chaitin: https://arxiv.org/pdf/math/0404335.pdf https://arxiv.org/pdf/math/0404335.pdf
    - adgjlsfhk1 3y ago
      to start, there are an uncountable infinity of them so 100% are not representable. Also when notated typically, there are 2 representations of lots of numbers.
    - ska 3y ago

cratermoon 3y ago

Here's a REALLY weird one: https://latkin.org/blog/2014/11/22/mullers-recurrence-roundoff-gone-wrong/ https://latkin.org/blog/2014/11/22/mullers-recurrence-roundo...

Pinus 3y ago

Ages ago, when I was in college, a "numerical methods" course was more or less obligatory... has that changed? 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.

bee_rider 3y ago
I assume it depends on the department. CS students are afraid of non-integer values, engineering students are afraid of codes longer than 200 lines. Just kidding, mostly. I do somewhat wonder if numerical methods have suffered from their success; like everybody uses LAPACK via Numpy but understanding it under the hood, that job always belongs to the next department over wherever you are, haha.
jmrm 3y ago
AFAIK in Europe numerical methods is mandatory in most Telco/EE degrees, but not in CS/SWE. The opposite happens with discrete maths btw
- anthk 3y ago
  Here in Spain if you can't grasp numerical methods in CS you can't even think on earning a degree.
  - rerdavies 3y ago
    Here in Canada, CS is part of the Math Department at the University of Waterloo, but part of the Philosophy Department at University of Toronto. Waterloo graduates take a Numerical Methods course in second year. I fear for the philosophers.
ArcticLandfall 3y ago

doubloon 3y ago

i think it is our educational system that is weird. 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. 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. 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.

galdosdi 3y ago
Great points, let me riff on that though -- what's really "weird" (but totally understandable in the historical context of underpowered computers until recently) is that we think of floats as the default representation of a real at all, when more reasonable modern default choices for all purposes outside of high performance applications would either be a rational composed of bigints, or a bigdecimal, or maybe a generating function for the coefficients of a Taylor series expansion or something like that if you really want to get nuts. In any case you can have as much precision as you can afford, although for anything but rationals it'll be an approximation always, and the question is just how many digits you care for. A lot of choices but the only one that's really absurd is the actual default of floats which are terrible for every purpose except their original raison detre, which was performance in a time when computers were so slow that buying a floating point unit to put in your computer was a normal desirable upgrade for a PC owner (and still crucial for high thoroughput computing, video games, graphics rendering, ML, etc) I do remember being taught the difference between rationals and irrationals in high school. Another aside... The really crazy difference is between the definable and non definable numbers. The first which includes everything from pi to every busy beaver number is countably finite just like the integers, only the nondefinables (numbers whose shortest mathematical definition would not be finite or not exist) give the reals their uncountably infinite nature

architectonic 3y ago

That's called the stability of the algorithm. The first algorithm is not stable whereas the second with the logarithm is. Another example: compare the two algorithms, that are actually pure-mathematically equal: def fc(x): return sqrt(x+1) - sqrt(x) def fd(x): return 1/(sqrt(x+1) + sqrt(x)) For large values of x (1e16 for example) plot the results and see the difference. 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])

echoangle 3y ago

Kind of OT: Why do so many people say “phenomena” when they want to use the singular case?