Math Keeps Changing
macwright.org
macwright.org
[0] Assuming round-to-nearest, of course. For directed rounding, the problematic case is a result that's arbitrarily close to some exactly-representable value, but not yet determined to be "apart" from it or not in the rounding sense.
https://docs.oracle.com/javase/8/docs/api/java/lang/Math.htm...
0.1 + 0.2 = 0.30000000000000004
> And this is the way that compiled languages have always worked: when you compile a C program, the methods you import from math.h are included in the compiled binary.
Depends on your inlining level and choice of linking.
It is inherently an issue with binary floating point[1]. It is not an issue with decimal floating point. In decimal floating point, 0.1 + 0.2 = 0.3 exactly.
JavaScript could have used decimal floating point instead – but that would be slower, since it usually has to be emulated in software – mainstream platforms like x86 and ARM lack decimal floating point hardware, it is only found on rare platforms such as POWER, System z, and Fujitsu SPARC – and even with hardware acceleration it is still going to be slower. JavaScript, like most languages, prioritises speed over obviousness in this area.
Alternatively, JavaScript could have supported both decimal and binary floating point, and beginners could stick to the more obvious former, and those needing high performance could switch to the later – but its designers decided (probably correctly) that the added complexity wasn't worth it. (Also, having two types of floating point numbers in a language would likely confuse beginners just as much as binary floating point does.)
[1] And hexadecimal floating point too, but that is rarely encountered outside of IBM mainframes, and also various other obscure obsolete systems such as Interdata 8/32 minicomputers
Decimal floating point doesn't fix the problem, it just makes the simplest example work. If you did ⅓ + ⅓ + ⅓ in decimal floating point, you'd still get 0.99999 instead of 1.
The gap between naive expectations of how it should work, and how it actually works, exists for both binary and decimal floating point, but the gap is bigger in the case of binary than it is in the case of decimal.
Also, the problem that ⅓ + ⅓ + ⅓ < 1 is not unique to floating point – the same issue occurs with fixed point arithmetic (unless your scale factor is a multiple of 3, which is not a common choice). In the general case, if one wants to avoid this problem, one needs to use rationals.
I think the issue with ⅓ + ⅓ + ⅓ is easier to explain to beginners – it requires an understanding of rounding and that ⅓ cannot be represented exactly as a decimal, which in many countries the school curriculum pays more attention to than it does to other number bases, so many people have a better intuitive understanding of the former than the later.
This is mentioned in the paragraph that is immediately after the text you excerpted.
It's not an issue for languages that have floating point numbers, it's an issue for languages that make them the default (or worse, the only option). There are plenty of other languages that let you use arbitrary-precision integers or decimals as well as the usual fixed-size floats and integers.
When you're writing a program intended to provide calculations to the end user (such as a calculator) you really ought to use arbitrary-precision decimals so that you never run into this problem.
Almost anything you'd want to use a non-integer number for is going to be difficult to nail down to 5 or 6 decimal digits of precision on (honestly, even 3 is pushing it). Double-precision floats gives you 11 decimal digits of precision--ample room for accumulated round-off error to not affect the final result, unless you've got numerical instabilities.
I'd be interested in a calculator that used arbitrary-precision decimals to compute arcsines, logs, etc. But I wouldn't recommend anyone create or use one.
The purity of the approach is tainted (in a helpful way, not a bad way) by a note at the end:
> The most recent versions of our calculator explicitly track rational multiples of π and some other common irrational constants. This allows us to compute a rational result for sin(π/6) in radian mode, as we already did for sin(30°).
When I mentioned arcsines above, that wasn't a coincidence -- radians are convenient, but they're also impossible to measure exactly, and when your computation goes through a trig function, you may unavoidably lose the ability to state your results exactly (except by invoking a function name like "arcsine").
That isn't as much of a theoretical obstacle as it might sound like -- you can't state square roots exactly either, except by referring to the concept "square root". But they can be a little more conceptually direct -- trig functions are inextricably bound up with pi, and it's hard -- impossible -- to recognize pi when it appears in your input, unless it appears as the symbol π.
They address that in the very next sentence...
https://tc39.es/ecma262/#sec-math.sin
function sin(x) {
if (x == -Infinity || x == Infinity) return NaN;
return x;
}
function cos(x) {
if (x == 0) return 1;
return sin(x);
}
function sqrt(x) {
if (x < 0) return NaN;
return x;
}
function exp(x) {
if (x == -Infinity) return 0;
if (x == Infinity) return x;
return cos(x);
}Firstly, you are missing ToNumber(): “Each of the following Math object functions applies the ToNumber abstract operation to each of its arguments (in left-to-right order if there is more than one). If ToNumber returns an abrupt completion, that Completion Record is immediately returned. Otherwise, the function performs a computation on the resulting Number value(s). The value returned by each function is a Number.”.
Secondly, no one would call those reasonable approximations, let alone approximations of standard C library functions: “these functions are intended to compute approximations to the results of familiar mathematical functions, but some latitude is allowed in the choice of approximation algorithms. The general intent is that an implementer should be able to use the same mathematical library for ECMAScript on a given hardware platform that is available to C programmers on that platform.”
> implementation-dependent approximation to the sine of x
That said, sin(x) being about x for "sufficiently small" x (about 0.1 or less) is a very well approximation. The error bound is given by the Maclaurin expansion of sin(x), which is x - x³/3! + x⁵/5! - …, so sin(x) = x has an error bound of x³/6.
I can help! The error is at most 2. ;D
> The error bound is given by the Maclaurin expansion of sin(x), which is x - x³/3! + x⁵/5! - …, so sin(x) = x has an error bound of x³/6.
I want to note that this is only true for small x. In general, the Maclaurin sequence is not absolutely decreasing (that is, the sequence formed by the absolute value of each term is not strictly decreasing -- x⁵/5! > x³/3! when x > 5), so the first unused term is not a strict bound on the total error.
But the web has a different culture where webpage and user agent collaborate. Scripts try to account for browser bugs, and browsers try to render even syntactically invalid pages.
What they might not care about is the provenance of your jokes.
The point of the parable is that, since the standard would seem to admit such a ludicrously incorrect implementation, the standard itself must be wrong.
I will point out, since you appear to believe contrarily, that my experience is that the JS and browser engine community have historically been far more accepting and blasé about loose specification than the C/C++ community.
Really? It doesn't seem obviously worse than deleting null pointer checks as C/C++ compilers do, which is widely accepted.
Most of the complaints of undefined behavior are actually people who write code with nonsense semantics, and are then surprised to find that the compiler comes out with a different nonsense semantics than what they intended.
The case I'm thinking of was some linux kernel code that looked roughly like:
void checkBeforeDangerousThing(SOMETHING* detailsPtr) {
SOMETHING details = *detailsPtr;
if(detailsPtr == null) return;
if(!isValid(details)) return;
doDangerousThing();
}
- For many transcendental functions, 1 ulp error bound is reasonable but not necessarily reproducible, and 0.5 ulp (optimal) error bound algorithm exists but is significantly slower, to the point that Kahan [1] has once pointed out that some versions of MATLAB hard-coded log10(10^n) cases to give a false illusion of accuracy.
- The optimal algorithm is, fortunately, computationally tractable. Many functions have been thoroughly analyzed for the closest outputs to IEEE 754 binary64 values and therefore the number of bits required to accurately round them, and that number of bits is small enough (~150 bits) that all trig/log/exp/pow functions can be computed without bignums.
- For a long time CRlibm from INRIA was the only library implementing optimal algorithms, and it was licensed in GPL. A successor to CRlibm, Metalibm, was in development with no licensing statements available, but nowadays a version of Metalibm [2] is licensed in MIT so you can probably make use of them.
- If I recall correctly, CRlibm had way too big code footprint (several hundreds of KBs when compiled) that I couldn't really use that anyway. The main culprit seems to be five sets of functions for each possible rounding mode. I haven't checked yet if a new Metalibm is still big.
- There are many approximation algorithms and libraries if you don't want the full accuracy. One notable example is Sleef [3] that has multiple functions for desired accuracy (from 0.5001 ulp to 3.5 ulp). Assuming that we will eventually have a set of small and portable libms necessary, I believe a variation of Sleef's approach is desirable for future programming languages: a (compile-time) argument for accuracy, all the way from the optimal (0.5 ulp) to the crude (2^20 ulp?).
[1] https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
The author's preferred solution is to use a cross-platform library that always gives the same answers. To which I ask the HN audience: what are actual valid use cases for that?
I'm trying to imagine where it is required that re-running the same computations on different platforms yields not just the same results in digits 1-15, but the same in the final digits 16-20. I mean, that is just an insane level of accuracy. And really... instead of relying on any result in those last 5 digits (accuracy which is obviously not there), would it be better practice just to truncate anyways?
The second is just a single bit (...2F68 vs ...2F67) so may be simply an accuracy increase in the underlying algorithm.
> instead of relying on any result in those last 5 digits (accuracy which is obviously not there), would it be better practice just to truncate anyways?
No. You have absolutely no idea what your truncation just did to your calculation.
This comes up where you have to be careful inverting a matrix because you can get massive cancellation that makes the accuracy of the lower bits quite significant.
In addition, if you are running a Monte Carlo simulation, you may have introduced a bias into your simulation.
IEEE 754 was written by numerical analysts--quite often in obtuse ways that complicate the hardware implementation dramatically. You have to understand the full mathematical rationale behind it before you tamper with it. And, even then, you do so at your own peril.
Calculating Math.tanh(0.1)
// Node 4
0.099667994624955*90234*
// Node 6
0.099667994624955*81907*
So I guess I don't understand the point you're making?
Sure I don't understand what truncation is doing... but all the digits after truncation appear to be arbitrarily different between platforms anyways, so I don't see what there is to rely on there? If it's a Monte Carlo simulation as you suggest, there's already potentially arbitrary bias in the last 5 digits, no?
I'm just suggesting that if you want repeatable results cross-platform with different implementations, doesn't truncation seem reasonable? Or perhaps rounding, to avoid bias?
0x3FB983D7795F_4140 == 0.09966799462495590234 0x3FB983D7795F_413A == 0.09966799462495581907
0x3FA2F684BDA1_2F68 == 0.03703703703703703498 0x3FA2F684BDA1_2F67 == 0.03703703703703702804
> but all the digits after truncation appear to be arbitrarily different between platforms anyways
Not random at all. If you follow a properly characterized algorithm and execute it properly on IEEE-754 hardware, the result will be the same as if you carried it out on infinite-precision hardware unless you get unlucky and hit a Table Maker's Dilemma instance. And, in that instance, the only difference between your algorithm and a slightly different algorithm should be 1ULP (unit in the last place).
5 ULPs is a bug or a garbage library. Both are possible when Javascript is involved.
> I'm just suggesting that if you want repeatable results cross-platform with different implementations, doesn't truncation seem reasonable? Or perhaps rounding, to avoid bias?
Yes. We call that IEEE-754--it has round-to-nearest-even and a guard bit.
Note by William Kahan from Berkeley about exponential functions and keeping the accuracy to about 1ULP.
"A Logarithm Too Clever by Half" https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
I wouldn't say that. There is definitely a trade-off between accuracy and speed, and that trade-off varies for different mathematical functions. And it's not unreasonable for some people to say "I can tolerate losing 10 ULPs if it means the code will run twice as fast." None of the current math library implementations do a good job of letting users explore those trade-offs, and many of them don't even bother to give a claimed accuracy result.
(Note that Sun's [1] java.lang.Math does give ULP requirements on math functions, and Math.tanh has a 2.5 ULP accuracy requirement--atan2, sinh, cost, and tanh are the only functions with a >1 ULP accuracy bound).
[1] And yes, I'm calling out Sun here and not Oracle, because it was Sun who designed the math library and reported the error bounds, and the standard math library that people use when they're not shelling out to libm is Sun's fdlibm implementation anyways.
The arbitrary bias you'd find is probably less than you'd find by rounding (truncating or otherwise).
The problem with truncation to make things repeatable cross-platform is it assumes "cross-platform" to be "any platform where the Nth to last digit is correct" when the standard explicitly makes no such guarantee about the Nth to last digit.
So, are we bemoaning that we don't have an standard for the language? Akin to common lisp? For my part, that would be awesome. Though, you are likely to just relearn everything they did.
There is a standard. [0]
[0] https://www.ecma-international.org/publications/standards/Ec...
I would not be shocked to find I am off on that claim.
And IEEE 754 has quite a number of edge cases that people find surprising, including things like operations only have to be reproducible up to a certain epsilon - an epsilon that is only defined to a minimum range in JavaScript.
I think this mindset is good for getting things done efficiently though. Unfortunately, getting things done efficiently has been going out of fashion lately.
I’m working my way through this book:
https://www.amazon.com/Data-Science-Scratch-Principles-Pytho...
by rewriting the examples in Swift:
https://github.com/melling/data-science-from-scratch-swift
I do like the author’s idea of creating his Simple Statistics library:
https://macwright.org/2012/06/26/simple-statistics.html
You definitely learn much more by building something, even if it’s less than perfect.