Item — HN Reader

Recently I was reading an article on the EML operator (exp-min-log) that uses exponentiation and logarithm to build elementary math functions including arithmetic operations. There was a table of results from testing across languages.

  Language         Result of 2 x 3    Error
  ---
  Node.js v25.3.0  6.000000000000000  0
  Python 3.9.6     6.000000000000001  8.88e-16
  PHP 8.5.1        6.000000000000001  8.88e-16
  Go 1.26.2        6.000000000000000  0
  Rust             6.000000000000001  8.88e-16

It was speculated that this miniscule margin of error, 1 ULP (unit in the last place), is likely due to the difference in how log() is implemented by the language. Supposedly Python, PHP, and Rust use LLVM's libm (C math library) but maybe Go and Node.js internally compile to CPU instructions directly? There was no evidence presented, so I was skeptical of this explanation.

2 comments

Python, PHP, and Rust do indeed use libm. Not LLVM's libm specifically, though, but rather just about anything they can find in runtime, and those libraries can have suboptimal accuracy.

Node uses V8, which notably implements a ton of math manually, because it needs the results to be deterministic across devices. It seemingly uses LLVM's libm, which IIRC promises 0.5 ulp for most operations. (https://github.com/v8/v8/blob/f24c62fbc342d616032734b714116e...)

Go avoids dynamic linking, so they also have their own implementation. (https://github.com/golang/go/blob/543ead71a8e7acc2bd6f326327...) They only promise 1 ulp, but I guess in this specific case it works out better than approximations used by the default libm on their system by pure chance?

Nice, thanks for tracking down those links, fascinating. Apparently being off by 1 ULP is within IEEE 754 floating-point spec, and "effectively negligible".

Correctly rounding transcendental functions is very difficult because of something called the table maker's dilemma, so the standard didn't want to impose a potentially extreme performance cost if they didn't have to.

Correctly rounded single precision functions at a reasonable cost are much easier these days. LLVM 18+ implements them, but the gnu libraries have been slower to improve.

"The Table-Maker's Dilemma"

https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT

This is the best source I could find, funny enough the author doesn't explain the name "table maker".

I suppose it's from the parable of a table maker who finds that one leg is too long, so they sand down that leg, only to find that another leg is too long, so they sand that down... slowly sanding away all the legs.

As a simple programmer I don't understand the dilemma. I do understand that sin, tan, sqrt, etc., take one argument, so you can guarantee a certain precision for them, at least in 32-bit floats, whereas a log in arbitrary base has two inputs, so its input space is huge and it's hard to guarantee anything, especially for 64-bit doubles.

If someone could write up a good explanation that could be its own HN post

> the author doesn't explain the name "table maker".

I suppose that parable fits, but also there is the dilema of someone making a table of (results) of transcendental functions. In the before times, most people got their sines and cosines and logs and what not by looking them up in tables. And someone has to make those tables. The example given where rounding to some specific number of digits depends on how many digits you calculate to is deeply unsatisfying, especially if you wanted to make a table for other people to use to lookup the results of functions.

The dilemma is hypothetical (because we now know the exact precision requirement for those functions over all 64-bit doubles) but let's assume that we don't know that and also we happen to use decimals instead. What if, some f(x) is known to be, say, 123.4999...9997 with 1,000 fractional digits? You have to calculate at least 1,000 decimal digits of f(x) in order to correctly round, before that we only know it's between 123.4999...999 and 123.5000...001. There seems no mathematical theory that specially prevents this annoying possibility, hence the dilemma.

> because we now know the exact precision requirement for those functions over all 64-bit doubles

Do you have a link or something on this?

I'm on a mobile right now so I can't give the exact paper link but you can learn about the CORE-MATH project [1].

[1] https://core-math.gitlabpages.inria.fr/

My assumption is it's referring to mathematical tables and the historical makers of them. https://en.wikipedia.org/wiki/Mathematical_table

I simply assumed the tables in question were log tables. (No, not that kind of log.)

I assume these results are all amd64. Is this the same on arm64?

The original article [^1] didn't specify the CPU, but I gathered the code examples (and added one for Lua, untested) if anyone is curious to try with arm64. The test case to compare the margin of error was eml_mul(2, 3).

JavaScript

    const eml = (x, y) => Math.exp(x) - Math.log(y);
    const emlLn = (x) => eml(1, eml(eml(1, x), 1));
    const emlMul = (x, y) => eml(emlLn(x) + emlLn(y), 1);
    const emlAdd = (x, y) => emlLn(eml(x, 1) * eml(y, 1));

Python

    import math

    def eml(x, y):
        return math.exp(x) - math.log(y)
    def eml_ln(x):
        return eml(1, eml(eml(1, x), 1))
    def eml_mul(x, y):
        return eml(eml_ln(x) + eml_ln(y), 1)
    def eml_add(x, y):
        return eml_ln(eml(x, 1) * eml(y, 1))

PHP

    function eml(float $x, float $y): float {
        return exp($x) - log($y);
    }
    function eml_ln(float $x): float {
        return eml(1, eml(eml(1, $x), 1));
    }
    function eml_mul(float $x, float $y): float {
        return eml(eml_ln($x) + eml_ln($y), 1);
    }
    function eml_add(float $x, float $y): float {
        return eml_ln(eml($x, 1) * eml($y, 1));
    }

    func eml(x, y float64) float64 {
        return math.Exp(x) - math.Log(y)
    }
    func emlLn(x float64) float64 {
        return eml(1, eml(eml(1, x), 1))
    }
    func emlMul(x, y float64) float64 {
        return eml(emlLn(x)+emlLn(y), 1)
    }
    func emlAdd(x, y float64) float64 {
        return emlLn(eml(x, 1) * eml(y, 1))
    }

Rust

    fn eml(x: f64, y: f64) -> f64 {
        x.exp() - y.ln()
    }   
    fn eml_ln(x: f64) -> f64 {
        eml(1.0, eml(eml(1.0, x), 1.0))
    }
    fn eml_mul(x: f64, y: f64) -> f64 {
        eml(eml_ln(x) + eml_ln(y), 1.0)
    }
    fn eml_add(x: f64, y: f64) -> f64 {
        eml_ln(eml(x, 1.0) * eml(y, 1.0))
    }

Lua

    local math = require("math")

    function eml(x, y)
        return math.exp(x) - math.log(y)
    end
    function eml_ln(x)
        return eml(1, eml(eml(1, x), 1))
    end
    function eml_mul(x, y)
        return eml(eml_ln(x) + eml_ln(y), 1)
    end
    function eml_add(x, y)
        return eml_ln(eml(x, 1) * eml(y, 1))
    end

---

[^1]: https://lilting.ch/en/articles/eml-single-operator-elementar...