Figuring out round, floor and ceil with integer division(blog.pkh.me)
blog.pkh.me
Figuring out round, floor and ceil with integer division
http://blog.pkh.me/p/36-figuring-out-round%2C-floor-and-ceil-with-integer-division.html
27 comments
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 in the past because of roundings not following the usual convention ("banker's rounding").
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.
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.
Ok sorry, I think that's me misinterpreting what banker's rounding mean. But I know the regulator will expect round(0.5) = 1, not 0.
Banker's rounding is "round to even", because it limits the amount of bias introduced by the rounding.
> But I know the regulator will expect round(0.5) = 1, not 0.
That's not necessarily the most helpful, as there's at least two tie-breakings which provide that, but differ in their treatment of negative numbers (round up or round to +∞, and round away from zero or round to ∞).
> But I know the regulator will expect round(0.5) = 1, not 0.
That's not necessarily the most helpful, as there's at least two tie-breakings which provide that, but differ in their treatment of negative numbers (round up or round to +∞, and round away from zero or round to ∞).
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.
So, consider the source of the complaint.
So, consider the source of the complaint.
The biggest annoyance with P3’s round is that you can’t specify a rounding mode.
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”.
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”.
You can specify a rounding mode if you use decimals instead of floats
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!
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.
Edit: Actually further tests with INT_MIN / - 1 did lead to a crash. Maybe the first result was due to constant folding or something.
Like all signed overflow, it's undefined behavior in C and C++, so what exactly happens is up to the gods of optimizing compilers.
In case divisor is always positive, you can use these variants:
int divF(int a, int b) { return a / b - (a % b < 0); }
int divC(int a, int b) { return a / b + (a % b > 0); }
> 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.
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.
I thought the same might be true, and a theorem prover agrees with us: https://gist.github.com/jepler/eda8401d1b613c44376dc8148fcf9...
Thank you for pointing this out. I added a note about this.
> 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.
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.
Floating point operations are not deterministic in other situations: https://randomascii.wordpress.com/2013/07/16/floating-point-...
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.
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.
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.
Thanks to the author for describing everything so clearly!
I wonder if there is a use case for this instead of fixed point math which is more precise than floating point math
As always, Hacker's Delight is a wonderful book for all integer arithmetic tricks and operations.
til: "ceil" is a verb
I'm intensely interested in this type of stuff, but the article throws me off right away. If I were do undertake this project, I wouldn't "define round the way POSIX does" and "integer division the way C11 does". Instead I'd reinspect the choices POSIX made and the choices C11 made, and then create routines that work however you need it to, i.e. that are configurable and let you control how it works.
Rounding--which btw I prefer to think of as rounding to a decimal place, not rounding to an integer--for example, can introduce biases because there are 10 digits, 0 through 9, but "11 spots on the dial", from round-down-to-0 to round-up-to-0 and therefore 5 is "the exact center" and always rounding 5 in one direction is going to bias your results in one direction, depending on your dataset. The article is breezing past this stuff.
I'm going to keep reading because I'm sure I'll learn some interesting things, but at the end of the first page my skin is crawling because I can see ambiguities in the examples I'm being shown and I keep thinking "wut?".
Rounding--which btw I prefer to think of as rounding to a decimal place, not rounding to an integer--for example, can introduce biases because there are 10 digits, 0 through 9, but "11 spots on the dial", from round-down-to-0 to round-up-to-0 and therefore 5 is "the exact center" and always rounding 5 in one direction is going to bias your results in one direction, depending on your dataset. The article is breezing past this stuff.
I'm going to keep reading because I'm sure I'll learn some interesting things, but at the end of the first page my skin is crawling because I can see ambiguities in the examples I'm being shown and I keep thinking "wut?".
> The article is breezing past this stuff.
That is because the author has a very specific goal which they explain upfront, and the exploration is done in furtherance of that goal. It is not a general treaty on rounding.
That is because the author has a very specific goal which they explain upfront, and the exploration is done in furtherance of that goal. It is not a general treaty on rounding.
This paper gives a very useful treatment of integer division and clarifies the bits you're finding confusing. https://www.microsoft.com/en-us/research/publication/divisio...
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