Log is non-monotonic in PHP and Lua

The general rule is that, in base b, you have finite positional ("decimal") representations of numbers p/q where q is a divisor of bˆn for some n. So e.g. 1/10 doesn't have a finite binary representation because 10 = 2 * 5, and that factor 5 ensures 10 is never a divisor of a power of 2. Likewise, 3 is coprime with 10, so you can't represent 1/3 in decimal.

Calculating correct rounding for other functions without hardware support is both challenging to implement and is often really slow. And that's before one gets into special functions. Or, even worse, inverses of special functions. Sometimes you might be lucky when those don't over/underflow around the edges of the domain

But I don't think anyone does this in practice. Simulations don't need perfect precision, anything involving real data has to deal with measurement errors which are usually much larger than f64 ULP scales, and then there's configurable arbitrary-precision for when it's really needed.

[0]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT, which coincidentally is also about the log function

If someone had the weird idea to write some formula with logarithms in a random base they should convert that formula to use binary logarithms, because there is no advantage in using another base.

Nonetheless, if one would want to include such a function, or any other function, one could reduce the argument to some small range using mathematical properties of the function and then some polynomial approximation should be determined, in order to use it to compute the function values.

If you determine a polynomial approximation for a function, whichever function it is, you can compute the entire function faster and more accurately than when you compute the function as a composition of other functions, one or more of which are computed by their own polynomial approximations.

In order to ensure that the function is rounded correctly, one must determine how many polynomial terms are needed and how big must be the numbers used for storing the coefficients.

For a function with a single argument, like "log2", it is much easier to do an exhaustive search to find all the argument values that would need more polynomial terms or more bits in the coefficients, to guarantee correct rounding.

For functions with 2 arguments, it can be difficult or impossible to find all problematic argument pairs, unless the function has some special mathematical properties that would allow a formal proof about the accuracy needed in polynomial evaluation to ensure correct rounding.

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

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

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

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

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

Which prints, for LuaJIT: true true false

And for Lua 5.5: true false false