
Math Keeps Changing - jashkenas
https://macwright.org/2020/02/14/math-keeps-changing.html
======
zozbot234
AIUI, the IEEE 754 standard does _not_ guarantee exactly-reproducible results
for transcendental operations, because of the "table-maker's dilemma" \- the
mathematical fact that given a result that's arbitrarily close to being
halfway in-between two floating-point values that differ by one ulp (unit in
the last place)[0], it may be computationally infeasible to determine the
"correctly" rounded result. So this outcome, while perhaps undesirable from an
intuitive POV, is not "wrong" in any real sense.

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

~~~
TAForObvReasons
There's a simpler explanation: IEEE 754 does not require the operations in
question to be defined. They are spelled out in Section 9 (Recommended
operations). There's some corner-case guidance (for example, divideByZero for
log(0)), but they are implementation-dependent and different JS engines
producing slightly different values are all in compliance with the actual
requirements.

~~~
titzer
It's worth pointing out that some programming languages do indeed dictate
precision for more advanced operations. E.g. Java requires Math.sin to be
precise to 1up:

[https://docs.oracle.com/javase/8/docs/api/java/lang/Math.htm...](https://docs.oracle.com/javase/8/docs/api/java/lang/Math.html#sin-
double-)

~~~
lifthrasiir
1 ulp precision means that there are at most two possible results from given
operations. Without an explicit guarantee (probably impossible, as this would
essentially pin the algorithm down) they are still not reproducible.

------
saagarjha
> JavaScript, on Twitter, gets a lot of heat for this behavior
    
    
      0.1 + 0.2 = 0.30000000000000004
    

I mean, there’s a lot of things JavaScript should get heat for, but this
really shouldn’t be one of them. This is inherently an issue with any language
that has floating point numbers.

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

~~~
skissane
> I mean, there’s a lot of things JavaScript should get heat for, but this
> really shouldn’t be one of them. This is inherently an issue with any
> language that has floating point numbers.

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

~~~
jcranmer
> 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.

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.

~~~
skissane
It doesn't completely fix the problem, but it makes the problem smaller.

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.

------
archgoon
JS Standards compliant implementations of some math functions:

[https://tc39.es/ecma262/#sec-math.sin](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);
      }

~~~
rcar
That's a pretty generous interpretation of

> implementation-dependent approximation to the sine of x

~~~
jcranmer
If you don't care what the error is, then 1 is a perfectly acceptable answer.
I have no clue what the error on that answer is, but you said you didn't care
about that, so...

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.

~~~
im3w1l
It sounds like you come from the C++ world, where the compiler is basically an
adversary playing gotcha games.

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.

~~~
jcranmer
Actually, my inspiration came from a story I heard (thirdhand) about a
numerical analysis professor who would give that response, set around 1970--
before C existed, let alone C++.

~~~
im3w1l
People obviously care about the error, _to some extent_ , in the javascript
world, and browser makers understand that and would never do something as
stupid as you propose even if it might technically be allowed.

What they might not care about is the provenance of your jokes.

~~~
jcranmer
Contrary to what you believe, no programming language community (not even
C/C++) would consider such a loose interpretation of the standard to be a
valid implementation.

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.

~~~
lmm
> Contrary to what you believe, no programming language community (not even
> C/C++) would consider such a loose interpretation of the standard to be a
> valid implementation.

Really? It doesn't seem obviously worse than deleting null pointer checks as
C/C++ compilers do, which is widely accepted.

~~~
jcranmer
What people forget is that, if the compiler didn't delete the null pointer
checks and executed the semantics exactly as the programmer wrote them... the
program is guaranteed to crash.

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.

~~~
lmm
Failure semantics are just as important as happy-path semantics. Good
languages are fail-stop.

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();
        }
    

The intended semantics of that are pretty clear. There are two defensible
implementations for this code if detailsPtr is null (hence why it was left
undefined in the C standard): either it should return cleanly without calling
doDangerousThing, or it should trap on the invalid null pointer dereference
(possibly asynchronously). What no reasonable language community would
consider to be a reasonable implementation - what only an adversarial reading
of the standard reveals is even an option - is to execute doDangerousThing()
when passed a null pointer.

------
lifthrasiir
I had once surveyed algorithms and libraries that implement accurate and
reasonably efficient mathematical operations. The status quo is a bit of mess.

\- 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](https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT)

[2] [https://github.com/kalray/metalibm](https://github.com/kalray/metalibm)

[3] [https://sleef.org/purec.xhtml](https://sleef.org/purec.xhtml)

------
crazygringo
Seems like the correct answer is exactly what the author admits at the end:
use an "epsilon" so that tests only compare against the expected result to a
certain number of decimal places.

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?

~~~
bsder
The first change is quite a bit more than a single bit (...413A vs ...4140),
that seems like a bug.

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.

~~~
crazygringo
I'm very confused by your reply... none of those example numbers you're using
(e.g. 2F68) are in the article. The article gives examples that are far more
different:

    
    
      Calculating Math.tanh(0.1)
      // Node 4
      0.099667994624955*90234*
      // Node 6
      0.099667994624955*81907*
    

The article is mostly about trigonometric functions etc., which the article
itself states IEEE 754 does _not_ define expected behavior for.

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?

~~~
bsder
Sorry, those are the hexadecimal representation of the underlying IEEE-754
binary floating point number:

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](https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT)

~~~
jcranmer
> 5 ULPs is a bug or a garbage library.

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.

------
graycat
As I recall, way back in the years of IBM's language PL/I, there was some high
quality documentation of the numerical methods of the code for the usual
special functions. Since the functions haven't changed and numerical methods
have been studied for literally centuries, this old IBM documentation might
still be relevant. When I saw that documentation, I guessed that it would be
good to have way over the horizon. Might be able to find that documentation
somewhere and see if it would still be helpful.

------
taeric
This feels more like complaining that JavaScript keeps changing. Because it
does. You can't expect what you wrote five years ago to run the same today.

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.

~~~
shakna
> So, are we bemoaning that we don't have an standard for the language?

There is a standard. [0]

[0] [https://www.ecma-
international.org/publications/standards/Ec...](https://www.ecma-
international.org/publications/standards/Ecma-262.htm)

~~~
taeric
Apologies, I meant a standard that was closer to common lisp. Since I believe
it is one of the more comprehensive standards.

I would not be shocked to find I am off on that claim.

~~~
shakna
This is more about the numeric tower, which common lisp does quite well,
whilst JavaScript leans pretty heavily onto IEEE 754. (rev 2008 in the above
linked spec).

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.

~~~
taeric
Right, and my assertion was supposed to be on just how extensive the common
lisp spec is.

------
cryptica
I had a similar problem. I have a very efficient garbage collector; if my mind
can't link the data to some concrete real-world purpose, then it will be
archived and eventually overwritten.

I think this mindset is good for getting things done efficiently though.
Unfortunately, getting things done efficiently has been going out of fashion
lately.

------
32gbsd
I never understood why PHP gets so much hate but Javascript gets a pass. maybe
times are a changing. But it might even be too late for PHP.

------
m4r35n357
Expecting a cobbled together mess to do accurate arithmetic. There's your
problem.

------
melling
With Data Science and Machine Learning, my interest in math has greatly
increased too. There’s definitely more of a need.

I’m working my way through this book:

[https://www.amazon.com/Data-Science-Scratch-Principles-
Pytho...](https://www.amazon.com/Data-Science-Scratch-Principles-
Python/dp/1492041130/ref=nodl_)

by rewriting the examples in Swift:

[https://github.com/melling/data-science-from-scratch-
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](https://macwright.org/2012/06/26/simple-statistics.html)

You definitely learn much more by building something, even if it’s less than
perfect.

