How to divide without dividing or division optimization for compilers(s)
If you've never tried to see how C++ code is turned into Assembly code by the compiler, you'll be in for a lot of surprises, and you don't need to look at some complicated source code full of templates or other complex constructs: consider the following snippet:
uint8_t div10(uint8_t x)
{
return x/10;
}Of course, I have already done this, and I will present the results right here, although I advise you to visit this wonderful resource yourself https://godbolt.org/ – set there, for example, x86-64 gcc 14.1 and make sure that the results are extremely interesting:
div10(unsigned char):
push rbp
mov rbp, rsp
mov eax, edi
mov BYTE PTR [rbp-4], al
movzx eax, BYTE PTR [rbp-4]
mov edx, -51
mul dl
shr ax, 8
shr al, 3
pop rbp
retIndeed, what is not visible at all in this piece of code is the instructions divwhich undoubtedly exists and actually performs division on the x86 architecture. But where did that magic constant come from, and it's negative! Well, okay, in fact it's positive and equal to 205 (since 205+51=256 = mod 256), so that question is closed, but how does it all work?
Basic principle
The way it works is that after multiplying by 205 we do a total right shift of 11 places, which is equivalent to dividing by 2^11 = 2048 and truncating the remainder. The last part is very important. Note that 205/2048 = 0.10009765625, in other words, slightly more than 0.1, if you multiply the last number on a calculator (or in Python) by 255 you get 25.524, in other words, after truncating the fractional part, this is the correct answer.
But could we take a number slightly less than 0.1? No – trivially multiplying by 10 we would get “slightly less than 1”, and simply discarding the fractional part we would get 0 – not quite the result you expect when dividing 10 by 10. And how much slightly larger should the number be for the trick to work? For uint8_t the maximum number is 255, so the number should differ from 1/10 by no more than 1/256. And with any number (at least from the base type uint8_t this is possible) – yes, with any.
Here, I will remind you about such mathematical functions as floor, ceil and trunc: we work only with positive numbers, therefore, without unnecessary complications, floor=trunc and simply discards the fractional part, and ceil always rounds up. Let d be our divisor, N is our digit capacity (we have 8), we want to get an expression of the form: m / 2 ^ (N + k) such that it is slightly greater than 1/d (here we think about everything in real numbers), or to be more precise:
m / (2 ^ (N + k)) – (1 / d) >= 1 / 2^N (just a generalization of the previous paragraph).
How to find the right numbers
Statement: m = ceil(2^(ceil(log(d)) + N) / d), k = ceil(log(d)) – how to understand all this, what is written here? … It can be understood, for example, in this way: the number at the bottom, according to the condition, is a power of 2, and this power is definitely not less than 2^N, then let's assume that I somehow found k – how can I now, given the given, select m? I already know the denominator – 2 ^ (N + k), I want to select an integer so that it is slightly greater than 1/d, what if I simply take trunc(2 ^ (N + k) / d)?
Let's use numbers from the example: I want to divide uint8_t by 10: i.e. N=8, d = 10, and for now let's take k equal to 2, for example, then trunc(2 ^ (N + k) / d) = trunk(2 ^ 10 / 10) = trunk (1024 / 10) = 102, and the whole expression m / (2 ^ (N + k)) = 102/ 1024 = 0.099609375, in general, close, but a little less than we need. But ceil will always be greater because: ceil(x) >= x for positive numbers. I haven't made this caveat yet, but I will, that d > 1 and d is not an exact power of two, that is, we won't have exact divisions here, the second caveat trunc (x/y) in C++ is just a regular integer division.
So, I hope by now it has become clear that m in the form I am looking for really approximates 1/d from above. Now let's see why I chose such k: ceil(2^(ceil(log(d)) + N) / d, here when doing the external ceil I add to the numerator a number no more than d – because the remainder can be neither greater nor equal to d, and it is clear that 2^(ceil(log(d)) > d.
I think it's time to show the code:
template<typename InputInteger, typename OutputInteger>
std::pair<OutputInteger, uint8_t> getDisionMultiplier(InputInteger divisor)
{
if (!divisor)
{
throw std::invalid_argument("Division by zero is impossible");
}
if (divisor == 1)
{
return {1,0};
}
constexpr uint8_t n = sizeof(InputInteger) * CHAR_BIT;
const double log_d_temp = std::log2(static_cast<double>(divisor));
const uint8_t log_d = std::ceil(log_d_temp);
if (log_d == std::floor(log_d_temp))
{
return {1, log_d};
}
OutputInteger res = std::ceil(static_cast<double>(static_cast<OutputInteger>(1) << (log_d + n)) / double(divisor));
return {res, n + log_d};
}
// somewhere in the main function
for(uint8_t divisor = 1; divisor > 0; divisor++)
{
auto [multiplier, shift] = getDisionMultiplier<uint8_t, uint16_t>(divisor);
for(uint8_t numenator = 1; numenator > 0; numenator++)
{
uint32_t res = static_cast<uint32_t>(numenator * multiplier) >> shift;
if (res != numenator / divisor)
{
std::cout << "panic: did something went wrong?" << std::endl;
}
}
}It is probably obvious that we will never see the phrase about panic – does that mean that’s it? It works and it’s time to finish the article? – No.
What does a compiler do?
The cospiler definitely does it differently, indeed, we deduce that the code given above gives for d=10:
auto p = getDisionMultiplier<uint8_t, uint16_t>(static_cast<uint8_t>(10));
std::cout << p.first << " " << (uint16_t)(p.second) << std::endl;410 12
And from the Assemly piece from the beginning of the article it is clear that it should have been 205 and 11… The thing is that gcc wants to get a constant m of the same size as d and uses a more cunning algorithm for this (if you look closely at my code above – I prudently used the uint16_t type).
The gcc algorithm is based on calculating a pair of constants m_low = trunc(2^(ceil(log(d)) + N) / d) and m_high = trunc(2^(ceil(log(d)) + N) + 2^(ceil(log(d)) / d), it is important to understand that m_low cannot be a real m (I showed this above), and m_high can (sorry, I will not try to describe this here, like all the remaining calculations), but m_high, generally speaking, is even greater than m from my naive algorithm. Yes, but both m_low and m_high are definitely less than 2^(N + 1), that is, for our case it would be 9 bits, and also exactly in integers m_high > m_low (without equality). What does this algorithm do next to make an 8-bit number from a 9-bit number? That's right, it will simply shift m_high one place to the right: it is important here understand that in this way it does trunc(m_high/2), and of course in parallel, it will decrease k (the number of shifts to the right in the final code), but … if the number is odd, it is not the same, is it? Yes, in this case there is a risk that we will start approximating from below… therefore, the compiler does this and trunc(m_low / 2) and compares these two numbers, because m_low is initially too small – if m_high has degraded to it – then it is impossible to decrease m_high and the corresponding shift further. And here is the code that does this:
template<typename InputInteger>
std::tuple<InputInteger, uint8_t, bool> getDisionMultiplier(InputInteger divisor)
{
if (!divisor)
{
throw std::invalid_argument("Division by zero is impossible");
}
if (divisor == 1)
{
return {1, 0, false};
}
constexpr uint8_t n = sizeof(InputInteger) * CHAR_BIT;
const double log_d_temp = std::log2(static_cast<double>(divisor));
const uint8_t log_d = std::ceil(log_d_temp);
if (log_d == std::floor(log_d_temp))
{
return {1, log_d, false};
}
uint64_t temp_low = (1UL << (log_d + n));
uint64_t temp_hight = (1UL << log_d) | (1UL << (log_d + n));
temp_hight /= divisor;
temp_low /= divisor;
uint8_t additionla_shift = log_d;
while (additionla_shift)
{
if (temp_low /2 >= temp_hight/2)
{
break;
}
temp_low /= 2;
temp_hight /= 2;
--additionla_shift;
}
return {temp_hight, n + additionla_shift, temp_hight > std::numeric_limits<uint8_t>::max()};
}
// somewhere in the main function
auto [coeff, shift, _] = getDisionMultiplier(static_cast<uint8_t>(10));
std::cout << (uint16_t)coeff << " " << (uint16_t)(shift) << std::endl;205 11
Success! The coefficients matched! Is that all? Alas: there is no guarantee that the loop that reduces temp_hight will not exit immediately after the first iteration. That is, there is no guarantee of getting an 8-bit number, but then we can have a modulo slice in return. Gcc can successfully use this sliced value – but this is clearly a separate topic.
And if you are really interested or just can't wait, then I will leave the links that I relied on
