> Everyone already knows floating-point operations are imprecise and it wouldn’t be fun to blog about.
Wellll... Yes and no. And I want to nitpick "FP ops are imprecise" because it's important sometimes.
It does come up in Lua - Lua uses 64-bit double floats for everything, _even array indexing_, because they have 53 bits of mantissa and they're guaranteed to represent all 32-bit integers with 100% precision.
I just opened Lua and got `2 ^ 32 == 4294967296.0` and `2 ^ 32 + 1 == 4294967297.0`. You can store 4 billion things in a Lua table and access them with double float indexes.
The usual "0.1 + 0.2 != 0.3" is _not_ imprecision. You could dedicate 1,000 bits to a float and still find some example of a number that's trivial to represent in decimal but repeats forever in binary.
While I'm on it, fixed-point and floating-point aren't magically different. Floats work better on very big values. Fixed-points are more predictable but run out of range easily when squaring numbers. This comes in 3D math when finding the length of vectors or normalizing vectors.
The PlayStation 1 didn't have jiggly vertices and warpy textures because of fixed-point. It had jiggly vertices because its GPU didn't have subpixel-precise rendering, and it only had 16 bits of precision. It had warpy textures because it had affine texture mapping. The N64's GPU had more bits of precision and it had perspective-correct texture mapping. There's no 16-bit floating-point format that would have saved the PS1 from wobbling.
Also floats aren't the cause of T-junctions or "sparklies" that are common in 3D game levels. That will happen in any system with finite precision, and if fixed-point would fix it, we'd use fixed-point.
On the other hand, decimal floating point would fail with 1/3 + 1/3 != 2/3.
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.
Thank you! I completely agree with this, floats are treated as a lot more magical than they actually are. I myself often rely on their exact guarantees for fun tricks (e.g. fast precise int<->float conversion: https://purplesyringa.moe/blog/fast-limited-range-conversion...). But while IEEE-754 is very precise, some subtleties arise when you add library functions to the mix -- many libm's and userland libraries don't guarantee 0.5 ulp precision for certain operations and don't document the guaranteed precision either, at which point you're left guessing and saying "well, I guess I should treat floats as magic in this case after all". I added "it's not imprecision" to step around this whole question.
> The usual "0.1 + 0.2 != 0.3" is _not_ imprecision.
It _is_ imprecision. Even if you cannot represent 0.3 exactly, you can consider the closest representable number, and that is different than the sum of the closest numbers to 0.1 and 0.2 because doing all of this introduced imprecision.
To expand on this, the imprecision arises the moment you write 0.1 (decimal) in source code and the compiler converts it to binary; the arithmetic itself is (mostly) exact. So as long as you use numbers that look "good enough" in base-2, floats behave very reasonably. The constantly made assumption that floats are imprecise is, in this sense, a user error -- the user shouldn't have used decimals.
> Lua uses 64-bit double floats for everything, _even array indexing_, because they have 53 bits of mantissa and they're guaranteed to represent all 32-bit integers with 100% precision.
Lua 5.3 introduced a proper integer type. However, JavaScript still represents all numbers as doubles!
There's a widespread misconception (not shared by the article author) that floating point arithmetic is "imprecise" in the sense that the result is off by some amount of random noise. On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.
For efficiency reasons, library functions (and sometimes, sadly, hardware implementations) don't always guarantee this (closest representable) exact result. That's fine if the function is then called "approximate_inverse_square_root" or something. Sometimes being off by some epsilon is a good tradeoff for 10x efficiency. I'm unhappy when such a tradeoff is smuggled into my math library without warning.
I checked Rust's source code to confirm the article's claim that it doesn't have a precise log function, and indeed, here's the implementation:
I'm sure other math libraries aren't better. But this is not a correct implementation of arbitrary-base logarithm. A function like that perhaps shouldn't even be offered in the standard library at all (since it's so trivial to begin with), or at the very least, not with that name. If a programmer wants to opt-in to a fast but slightly wrong value, they should do so explicitly, in my opinion.
> On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.
Not sure about Rust, but C++ only guarantees that for four arithmetic operations and `sqrt`. Not `pow`, not `log`. Not even when when all inputs/outputs are integers. Sometimes you can even get `(int)pow(10, 2) == 99`, which I believe is fully standard-compliant.
Yeah, those operations (+fma) are the only ones for which IEEE 754 mandates correct rounding.
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
It's still challenging to implement with hardware support but you don't notice because it challenges Intel instead of you. There is a reason that no instruction set since x87 implements transcendentals - and x87 did them in microcode.
Older numerical libraries such as cephes had to work with heterogeneous non-IEEE floating point implementations: you can still see remnants of VAX tests in the docs, for example. That, too, must've been quite a struggle
Thank you! This point has been driving me mad for the last couple of years.
While answering a sibling comment, I found a paper analyzing various libms for their precision: https://homepages.loria.fr/pzimmermann/papers/glibc238-20230... (2023). I only skimmed it, but it looks like only LLVM's libm guarantees 0.5 ulp for every supported single-precision operation, and everyone else is wildly off. That's a pretty good result, though -- it means there's at least one reasonably compliant libm :)
Another sibling comments says that guaranteeing 0.5 ulp for double-precision operations is nigh impossible (and then another says it is after all). I don't have the knowledge to confirm which is true, but it's possible that this is the best we can get.
Table--Maker's dilemma [0] says that we probably can't get correct rounding for some functions. That said, you can still calculate every function within 1 ULP by using extended precision (for example the double-double arithmetic, which uses two doubles to represent a single fp number) and rounding.
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.
I am curious what the correct algorithm would be, if for no other reason than to find out how on god's green earth it could be slower than using a division as an erstwhile shortcut. 8I
In this particular case, there never exists any need to compute a logarithm of an arbitrary base, so such a function must not exist in a library.
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.
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.
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.
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.
> CORE-MATH Mission: provide on-the-shelf high performance open-source mathematical functions with correct rounding that can be integrated into current mathematical libraries (GNU libc, Intel Math Library, AMD Libm, Newlib, OpenLibm, Musl, Apple Libm, llvm-libc, Microsoft libm, CUDA libm, ROCm)
That's beautiful. Here's the Git repo with implementation.
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).
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
Yeah, I simplified it a little. Though I must say I'm surprised pretty much every library I looked at uses the natural logarithm specifically, and not log2, which would seemingly be easier to compute with floats. Does anyone here know why, by any chance?
Not sure why they wouldn't use log2 internally, but the mathematical reason is that e^x is the function whose derivative is itself, so it and the natural logarithm are fundamental to lots of other theories/operations like Euler's formula.
The only advantage of the hyperbolic logarithm (a.k.a. natural logarithm) is that the multiplicative factor is unity when computing its derivative or primitive.
The binary logarithm can be computed faster and more accurately, except for arguments very close to 1.
It is a historical accident that most standard libraries have been defined in the past to include natural logarithms instead of binary logarithms. Moreover lazy implementers sometimes have written a binary logarithm function that uses a pre-existing hyperbolic logarithm function, instead of the direct implementation that would be much more efficient.
Derivatives and primitives are seldom computed, but the computational advantage of binary logarithms applies to each function evaluation.
Even when derivatives or primitives are computed, in almost all applications there already exists another multiplicative constant in the formula that must be computed, which can absorb in it the factor "ln 2", so there is no increased computational cost.
There exists absolutely no reason to ever use multiple bases with logarithms.
One can choose to always use only binary logarithms and there exists no problem that would need other kinds of logarithms.
Choosing a base for logarithms is equivalent with choosing a unit of measurement for them. Changing the base of logarithms is done by the same rules as changing the unit of measurement for any quantity. In the same way as there is no need to use multiple units of measurement for a quantity, there is no need to ever work with multiple bases for logarithms.
Thus there is no reason for a programming language to include a function with 2 arguments like "log(arg, base)" instead of having only a function with 1 argument, e.g. "log2" (and the 2 constants LN_2 and LOG2_E).
You can change base with a single multiplication by special constants if you accept 1ULP error, so working in one or the other doesn't really gain you much. But exp/logs are often implemented as a combination of a table lookup and polys, so you can just bake it into the constants even more easily.
I have the impression that creating monotonic floating point function is particularly hard. At least I know that std::lerp is designed to be monotonic in its t parameter (not even in a or b) and its implementation has a number of branches. Although it has other guarantees as well, so I don't know how much of the implementation complexity can be assigned to the monotonicity requirement.
I sort of needed to do that. I needed to compress data with an approximately geometric distribution, and as part of that process I needed to invert its CDF, which is `CDF = 1 - (1 - p)^x`. That translates to `x = log_(1 - p) (1 - CDF)`, which is variable over both the argument and the base. At that point I wondered how consistent the implementation of double-argument `log` is, since the encoder and the decoder need to agree about it, which led to this article.
The log base is a parameter, not a function, so that doesn't typecheck. Multivariate monotonicity isn't really a popular term, as far as I'm aware, so I'm not aware of good terminology for this. Maybe "log is non-monotonic with respect to the base"?
That's interesting! On my PC it reproduces under Lua 5.5 (`log_a x > log_b x`), but not under LuaJIT (`log_a x = log_b x`). I took a look at LuaJIT's implementation (https://github.com/LuaJIT/LuaJIT/blob/faaf663340347a78b22ed9...) and noticed that it always uses the `ln x / ln a` formula -- or, rather, `log_2 x / log_2 a`, which is just as correct I guess. Have you perhaps misinterpreted the results? (I do think it's valuable to add that this doesn't apply to LuaJIT though.)
Hm, perhaps, but I am confused now. In the article, you're comparing log_a x < log_b x, but now you're comparing log_a x > log_b x in your comment. To make it clear, I am running this code: https://bpa.st/ZCKA
`log_a x < log_b x` is the mathematically correct result. Due to inherent rounding from floats being a finite representation of real numbers, `log_a x <= log_b x` would be expected from any correct implementation using floats. So either `true true false` or `true false true` would be reasonable. (I made a mistake in the previous comment, LuaJIT returns `<` for me, just like in your comment, not `=`.)
The topic of the post is that in PHP and Lua (without LuaJIT), sometimes this inequality doesn't hold, and instead we get `log_a x > log_b x`, which is very incorrect and cannot be explained away by rounding.
> Everyone already knows floating-point operations are imprecise and it wouldn’t be fun to blog about.
Wellll... Yes and no. And I want to nitpick "FP ops are imprecise" because it's important sometimes.
It does come up in Lua - Lua uses 64-bit double floats for everything, _even array indexing_, because they have 53 bits of mantissa and they're guaranteed to represent all 32-bit integers with 100% precision.
I just opened Lua and got `2 ^ 32 == 4294967296.0` and `2 ^ 32 + 1 == 4294967297.0`. You can store 4 billion things in a Lua table and access them with double float indexes.
The usual "0.1 + 0.2 != 0.3" is _not_ imprecision. You could dedicate 1,000 bits to a float and still find some example of a number that's trivial to represent in decimal but repeats forever in binary.
While I'm on it, fixed-point and floating-point aren't magically different. Floats work better on very big values. Fixed-points are more predictable but run out of range easily when squaring numbers. This comes in 3D math when finding the length of vectors or normalizing vectors.
The PlayStation 1 didn't have jiggly vertices and warpy textures because of fixed-point. It had jiggly vertices because its GPU didn't have subpixel-precise rendering, and it only had 16 bits of precision. It had warpy textures because it had affine texture mapping. The N64's GPU had more bits of precision and it had perspective-correct texture mapping. There's no 16-bit floating-point format that would have saved the PS1 from wobbling.
Also floats aren't the cause of T-junctions or "sparklies" that are common in 3D game levels. That will happen in any system with finite precision, and if fixed-point would fix it, we'd use fixed-point.
> The PlayStation 1 didn't have jiggly vertices and warpy textures because of fixed-point
Indeed modern GPUs still use fixed-point rasterizers.
> Lua uses 64-bit double floats for everything
Not true in Lua >= 5.3. 64-bit signed integers are included in the number type now. And there is a math.ult function for treating them as unsigned.
And so 0x7fffffffffffffff is representable properly.
this is true but basically everyone uses LuaJIT these days
Pinned several language versions ago? No they don't.
Citation needed. You only need LuaJIT if performance is very important, which is not the case for many scripting purposes. I personally use Lua 5.3.
Well, I use PICO-8 LUA. Where are my double-precision floats and 64bit integers at? XD
> The usual "0.1 + 0.2 != 0.3" is _not_ imprecision.
Correct. That’s using the wrong base. Decimal floating point doesn’t have that problem.
On the other hand, decimal floating point would fail with 1/3 + 1/3 != 2/3.
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.
Thank you! I completely agree with this, floats are treated as a lot more magical than they actually are. I myself often rely on their exact guarantees for fun tricks (e.g. fast precise int<->float conversion: https://purplesyringa.moe/blog/fast-limited-range-conversion...). But while IEEE-754 is very precise, some subtleties arise when you add library functions to the mix -- many libm's and userland libraries don't guarantee 0.5 ulp precision for certain operations and don't document the guaranteed precision either, at which point you're left guessing and saying "well, I guess I should treat floats as magic in this case after all". I added "it's not imprecision" to step around this whole question.
> The usual "0.1 + 0.2 != 0.3" is _not_ imprecision.
It _is_ imprecision. Even if you cannot represent 0.3 exactly, you can consider the closest representable number, and that is different than the sum of the closest numbers to 0.1 and 0.2 because doing all of this introduced imprecision.
To expand on this, the imprecision arises the moment you write 0.1 (decimal) in source code and the compiler converts it to binary; the arithmetic itself is (mostly) exact. So as long as you use numbers that look "good enough" in base-2, floats behave very reasonably. The constantly made assumption that floats are imprecise is, in this sense, a user error -- the user shouldn't have used decimals.
> Lua uses 64-bit double floats for everything, _even array indexing_, because they have 53 bits of mantissa and they're guaranteed to represent all 32-bit integers with 100% precision.
Lua 5.3 introduced a proper integer type. However, JavaScript still represents all numbers as doubles!
There's a widespread misconception (not shared by the article author) that floating point arithmetic is "imprecise" in the sense that the result is off by some amount of random noise. On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.
For efficiency reasons, library functions (and sometimes, sadly, hardware implementations) don't always guarantee this (closest representable) exact result. That's fine if the function is then called "approximate_inverse_square_root" or something. Sometimes being off by some epsilon is a good tradeoff for 10x efficiency. I'm unhappy when such a tradeoff is smuggled into my math library without warning.
I checked Rust's source code to confirm the article's claim that it doesn't have a precise log function, and indeed, here's the implementation:
I'm sure other math libraries aren't better. But this is not a correct implementation of arbitrary-base logarithm. A function like that perhaps shouldn't even be offered in the standard library at all (since it's so trivial to begin with), or at the very least, not with that name. If a programmer wants to opt-in to a fast but slightly wrong value, they should do so explicitly, in my opinion.
> On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.
Not sure about Rust, but C++ only guarantees that for four arithmetic operations and `sqrt`. Not `pow`, not `log`. Not even when when all inputs/outputs are integers. Sometimes you can even get `(int)pow(10, 2) == 99`, which I believe is fully standard-compliant.
Yeah, those operations (+fma) are the only ones for which IEEE 754 mandates correct rounding.
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
It's still challenging to implement with hardware support but you don't notice because it challenges Intel instead of you. There is a reason that no instruction set since x87 implements transcendentals - and x87 did them in microcode.
True!
Older numerical libraries such as cephes had to work with heterogeneous non-IEEE floating point implementations: you can still see remnants of VAX tests in the docs, for example. That, too, must've been quite a struggle
Thank you! This point has been driving me mad for the last couple of years.
While answering a sibling comment, I found a paper analyzing various libms for their precision: https://homepages.loria.fr/pzimmermann/papers/glibc238-20230... (2023). I only skimmed it, but it looks like only LLVM's libm guarantees 0.5 ulp for every supported single-precision operation, and everyone else is wildly off. That's a pretty good result, though -- it means there's at least one reasonably compliant libm :)
Another sibling comments says that guaranteeing 0.5 ulp for double-precision operations is nigh impossible (and then another says it is after all). I don't have the knowledge to confirm which is true, but it's possible that this is the best we can get.
Table--Maker's dilemma [0] says that we probably can't get correct rounding for some functions. That said, you can still calculate every function within 1 ULP by using extended precision (for example the double-double arithmetic, which uses two doubles to represent a single fp number) and rounding.
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
I am curious what the correct algorithm would be, if for no other reason than to find out how on god's green earth it could be slower than using a division as an erstwhile shortcut. 8I
In this particular case, there never exists any need to compute a logarithm of an arbitrary base, so such a function must not exist in a library.
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.
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.
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.
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
My assumption is it's referring to mathematical tables and the historical makers of them. https://en.wikipedia.org/wiki/Mathematical_table
> 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/
> CORE-MATH Mission: provide on-the-shelf high performance open-source mathematical functions with correct rounding that can be integrated into current mathematical libraries (GNU libc, Intel Math Library, AMD Libm, Newlib, OpenLibm, Musl, Apple Libm, llvm-libc, Microsoft libm, CUDA libm, ROCm)
That's beautiful. Here's the Git repo with implementation.
https://gitlab.inria.fr/core-math/core-math/-/blob/master/RE...
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
Python
PHP
Go
Rust
Lua
---
[^1]: https://lilting.ch/en/articles/eml-single-operator-elementar...
I was prepared to read how its console logging function could reorder writes so that
could result in
or something, which would've been no less surprising.
I find it interesting the blog author tied PHP and Lua together. PHP uses a JIT from Lua. Is this related to the log issue?
> PHP uses a JIT from Lua.
Wow, TIL! I was certain this is a mistake, but apparently PHP uses DynAsm (https://wiki.php.net/rfc/jit), which was developed for LuaJIT (https://luajit.org/dynasm.html). Cool stuff!
To answer your question, probably not. I dated the PHP change that added this "optimization" back to 2014 (https://github.com/php/php-src/commit/b547e1358d3846fad4cd0c...), while DynAsm only started being used around 2019 (https://wiki.php.net/rfc/jit). I think this is just convergent evolution.
(I'm putting "optimization" in quotes because it changes semantics.)
They need not be base e. Any base will do as long as it is the same for both the numerator and the denominator
Yeah, I simplified it a little. Though I must say I'm surprised pretty much every library I looked at uses the natural logarithm specifically, and not log2, which would seemingly be easier to compute with floats. Does anyone here know why, by any chance?
Not sure why they wouldn't use log2 internally, but the mathematical reason is that e^x is the function whose derivative is itself, so it and the natural logarithm are fundamental to lots of other theories/operations like Euler's formula.
The only advantage of the hyperbolic logarithm (a.k.a. natural logarithm) is that the multiplicative factor is unity when computing its derivative or primitive.
The binary logarithm can be computed faster and more accurately, except for arguments very close to 1.
It is a historical accident that most standard libraries have been defined in the past to include natural logarithms instead of binary logarithms. Moreover lazy implementers sometimes have written a binary logarithm function that uses a pre-existing hyperbolic logarithm function, instead of the direct implementation that would be much more efficient.
Derivatives and primitives are seldom computed, but the computational advantage of binary logarithms applies to each function evaluation.
Even when derivatives or primitives are computed, in almost all applications there already exists another multiplicative constant in the formula that must be computed, which can absorb in it the factor "ln 2", so there is no increased computational cost.
There exists absolutely no reason to ever use multiple bases with logarithms.
One can choose to always use only binary logarithms and there exists no problem that would need other kinds of logarithms.
Choosing a base for logarithms is equivalent with choosing a unit of measurement for them. Changing the base of logarithms is done by the same rules as changing the unit of measurement for any quantity. In the same way as there is no need to use multiple units of measurement for a quantity, there is no need to ever work with multiple bases for logarithms.
Thus there is no reason for a programming language to include a function with 2 arguments like "log(arg, base)" instead of having only a function with 1 argument, e.g. "log2" (and the 2 constants LN_2 and LOG2_E).
You can change base with a single multiplication by special constants if you accept 1ULP error, so working in one or the other doesn't really gain you much. But exp/logs are often implemented as a combination of a table lookup and polys, so you can just bake it into the constants even more easily.
I have the impression that creating monotonic floating point function is particularly hard. At least I know that std::lerp is designed to be monotonic in its t parameter (not even in a or b) and its implementation has a number of branches. Although it has other guarantees as well, so I don't know how much of the implementation complexity can be assigned to the monotonicity requirement.
In the past there were a lot of standard libraries that computed the transcendental functions with big errors.
Now there are several libraries that are guaranteed to produce correctly rounded results for any FP64 arguments.
However, the correct libraries are somewhat slower than the libraries that provide correct results for most, but not for all input arguments.
Good find, but I'm curious if anyone is ever sweeping the base of the logarithm in real life.
I sort of needed to do that. I needed to compress data with an approximately geometric distribution, and as part of that process I needed to invert its CDF, which is `CDF = 1 - (1 - p)^x`. That translates to `x = log_(1 - p) (1 - CDF)`, which is variable over both the argument and the base. At that point I wondered how consistent the implementation of double-argument `log` is, since the encoder and the decoder need to agree about it, which led to this article.
Wouldn't a more applicable title be that the log base is non-monotonic?
The log base is a parameter, not a function, so that doesn't typecheck. Multivariate monotonicity isn't really a popular term, as far as I'm aware, so I'm not aware of good terminology for this. Maybe "log is non-monotonic with respect to the base"?
I've tested this with both LuaJIT and Lua 5.5. It appears to only affect LuaJIT. Worth mentioning in the article I'd say.
That's interesting! On my PC it reproduces under Lua 5.5 (`log_a x > log_b x`), but not under LuaJIT (`log_a x = log_b x`). I took a look at LuaJIT's implementation (https://github.com/LuaJIT/LuaJIT/blob/faaf663340347a78b22ed9...) and noticed that it always uses the `ln x / ln a` formula -- or, rather, `log_2 x / log_2 a`, which is just as correct I guess. Have you perhaps misinterpreted the results? (I do think it's valuable to add that this doesn't apply to LuaJIT though.)
Hm, perhaps, but I am confused now. In the article, you're comparing log_a x < log_b x, but now you're comparing log_a x > log_b x in your comment. To make it clear, I am running this code: https://bpa.st/ZCKA
Which prints, for LuaJIT: true true false
And for Lua 5.5: true false false
`log_a x < log_b x` is the mathematically correct result. Due to inherent rounding from floats being a finite representation of real numbers, `log_a x <= log_b x` would be expected from any correct implementation using floats. So either `true true false` or `true false true` would be reasonable. (I made a mistake in the previous comment, LuaJIT returns `<` for me, just like in your comment, not `=`.)
The topic of the post is that in PHP and Lua (without LuaJIT), sometimes this inequality doesn't hold, and instead we get `log_a x > log_b x`, which is very incorrect and cannot be explained away by rounding.
Does that make more sense?
Ahem. Clearly my math knowledge failing me here. Apologies, and thanks for clearing it up :)