I was surprised by the observation that Python rounds towards even numbers, so round(0.5) = 0, but round(1.5) = 2.
This is indeed clearly documented [1], I guess I never looked closely enough. I found some discussion on the dev-python list [2] which shows at least I'm not the only one surprised by this!
That's good to know and kind of dangerous. I have seen banking regulators getting upset in the past because of roundings not following the usual convention ("banker's rounding"). You will have a hard time explaining it's an unexpected behaviour of some programming language rather than the management of the bank trying to be cute with the numbers.
I have seen banking regulators getting upset because they (effectively) assumed infinite precision until the "final rounding" (their verbiage). This was at a major financial services provider with exactly zero numerical stability analysts.
A fun fact about integer division is that on x86 dividing by zero is CPU exception 0. Except, it's not just for dividing by zero, it's whenever a division is impossible to perform.
So, for example what happens if you divide -1 by -INT_MIN? As you probably know, Abs(INT_MIN) is larger than INT_MAX, and so it is not possible to perform -1 / INT_MIN, as well as -1LL / INT64_MIN. It will crash your program, your emulator/sandbox and your server. So, be careful!
This had me scratching my head for a while. Shouldn't -1 / INT_MIN just be 0? Indeed it should. INT_MIN / -1 is what causes the issue. I also tested doing that the latter just to be sure, and on my system, performing that division led to a result of INT_MIN rather than an exception.
Edit: Actually further tests with INT_MIN / - 1 did lead to a crash. Maybe the first result was due to constant folding or something.
> We could stop right here but this suffers from overflow limitations if translated into C.
FWIW, the final version also suffers from integer overflow limitations. If the difference between an and INT_MIN/INT_MAX (depending on whether you floor or ceil) is <= b/2, you will have integer overflow.
> transforming a float based algorithm to integers in order to make it bit-exact
Modern CPUs, and most programming languages, guarantee bit-exactness for floats most of the times, because they implement IEEE 754.
Relatively easy to break (approximate instructions like rsqrtps, FMA instructions, compiler switches like -ffast-math or /fp:fast, MXCSR rounding modes) but I think fixing bugs is probably easier than refactoring everything from FP into fixed-point.
BTW, multiplication/division/rounding instructions are much faster for floats. I think a fixed-point version will only be a performance win when the integers are 8 or 16 bits wide, and the algorithm is vectorizable.
This can be terribly hard to manage in tests if you don't have a threshold infrastructure in place (which is hard and sometimes impossible to setup if you don't want to preserve huge references).
Also, and this is a bit off topic, even if I remained in the float domain for my use case, I would have to re-implement the lib math routines anyway (in my case cbrt and pow) because they are not implemented the same in all the libc, making them again victim to non determinism.
So far I observed better performance with integer arithmetic anyway, but that's not the reason that motivated me.
Also note that there are other situations were you may want to keep integers; typically in the kernel, or on arch were the floating point is not optimal or even available.
> Floating point operations are not deterministic in other situations
That article was written in 2013. A decade has passed, and in modern world these situations became rather rare.
x87 FPU is deprecated and not used by 64-bit programs. Returns values are passed in XMM0 register. Temporary values are either FP32 or FP64. Modern C++ compilers evaluate FP expression from left to right, and they don't use FMA unless explicitly allowed with -ffast-math or /fp:fast compiler switch. Programmers rarely using estimate instructions because modern CPUs complete precise division and square root in 10-15 cycles, little profit from the faster approximations.
The only variable thing which remains is MXCSR register in the thread state.
> I would have to re-implement the lib math routines anyway (in my case cbrt and pow)
Or you can copy-paste them from OpenBSD. The license is permissive, they are implemented in C without assembly, and the code quality is pretty good.
Readit News
This is indeed clearly documented [1], I guess I never looked closely enough. I found some discussion on the dev-python list [2] which shows at least I'm not the only one surprised by this!
[1] https://docs.python.org/3/library/functions.html#round
[2] https://groups.google.com/g/dev-python/c/VNf8TABiB9k
That… is exactly the behaviour python implements, and what GP was surprised by…
It is also, incidentally, the IEEE 754 recommendation.
Python 2 used to round away from 0.
So, consider the source of the complaint.
Well that and it returns an integer when rounding to 0 digits, but that’s less of an issue since the langage also switched to “true division”.
So, for example what happens if you divide -1 by -INT_MIN? As you probably know, Abs(INT_MIN) is larger than INT_MAX, and so it is not possible to perform -1 / INT_MIN, as well as -1LL / INT64_MIN. It will crash your program, your emulator/sandbox and your server. So, be careful!
Edit: Actually further tests with INT_MIN / - 1 did lead to a crash. Maybe the first result was due to constant folding or something.
FWIW, the final version also suffers from integer overflow limitations. If the difference between an and INT_MIN/INT_MAX (depending on whether you floor or ceil) is <= b/2, you will have integer overflow.
Modern CPUs, and most programming languages, guarantee bit-exactness for floats most of the times, because they implement IEEE 754.
Relatively easy to break (approximate instructions like rsqrtps, FMA instructions, compiler switches like -ffast-math or /fp:fast, MXCSR rounding modes) but I think fixing bugs is probably easier than refactoring everything from FP into fixed-point.
BTW, multiplication/division/rounding instructions are much faster for floats. I think a fixed-point version will only be a performance win when the integers are 8 or 16 bits wide, and the algorithm is vectorizable.
This can be terribly hard to manage in tests if you don't have a threshold infrastructure in place (which is hard and sometimes impossible to setup if you don't want to preserve huge references).
Also, and this is a bit off topic, even if I remained in the float domain for my use case, I would have to re-implement the lib math routines anyway (in my case cbrt and pow) because they are not implemented the same in all the libc, making them again victim to non determinism.
So far I observed better performance with integer arithmetic anyway, but that's not the reason that motivated me.
Also note that there are other situations were you may want to keep integers; typically in the kernel, or on arch were the floating point is not optimal or even available.
That article was written in 2013. A decade has passed, and in modern world these situations became rather rare.
x87 FPU is deprecated and not used by 64-bit programs. Returns values are passed in XMM0 register. Temporary values are either FP32 or FP64. Modern C++ compilers evaluate FP expression from left to right, and they don't use FMA unless explicitly allowed with -ffast-math or /fp:fast compiler switch. Programmers rarely using estimate instructions because modern CPUs complete precise division and square root in 10-15 cycles, little profit from the faster approximations.
The only variable thing which remains is MXCSR register in the thread state.
> I would have to re-implement the lib math routines anyway (in my case cbrt and pow)
Or you can copy-paste them from OpenBSD. The license is permissive, they are implemented in C without assembly, and the code quality is pretty good.