udiv128by64to64default not adding back

[haberg-1]

The function “udiv128by64to64default” in file udivmodti4.c on the link below does not add back after the reductions of q1 and q0, which causes an error of 1 with probability 2/𝛽, where 𝛽 is the base.

https://raw.githubusercontent.com/llvm/llvm-project/refs/heads/main/compiler-rt/lib/builtins/udivmodti4.c

[efriedma-quic]

Can you give an example of an input that doesn't produce the correct result?

[haberg-1]

See the 64-32 bit implementation below, which has the add-back step, but also examples where it is required. In the reference you give in your implementation file, it is the steps D5–D6 that are missing.

https://raw.githubusercontent.com/hcs0/Hackers-Delight/master/divmnu64.c.txt

[efriedma-quic]

Please give an example of an input to __udivmodti4, for the sake of anyone looking at this who isn't deeply familiar with the underlying math.

CC @legrosbuffle

[haberg-1]

I can't run your function. It certainly occurs when the divisor is 3 words, as in the given example. I have found no reference saying that it cannot occur when the divisor is 2 words, which is what you are implicitly assuming. Please let me know if you find such a reference.

[efriedma-quic]

I'm not saying you're wrong, I just don't have a day to invest into proving the existence of a theoretical issue.

https://godbolt.org/z/7cMMnoYsG if you want to try.

[haberg-1]

It is a must-have theoretical issue only if you do not implement the add-back. You might check with GCC; they may use the same method. When I test my own code, which covers the general case and implementation differs from yours, I use randomly generated numbers in the millions at a time, among other methods. I could not see how to do that on your link.

[haberg-1]

You might check with the author of the code:
https://danlark.org/2020/06/14/128-bit-division/

[efriedma-quic]

Couldn't figure out the right GitHub handle at first, but dug a bit more... CC @danlark1 . (Author of 68c011a .)

[danlark1]

I'll take a look but if you can provide an example, it would be great. If there is a mistake, I'd be curious to know how it evaded so much testing and such a wide usage for AArch64, for example

[haberg-1]

All examples seem to use at least a 3-word divisor. You might check your code yourself against binary division, which is easy to implement though slower, and which may have been what you had before. Random numbers can be generated by the example on this page:
https://en.cppreference.com/w/cpp/numeric/random/uniform_int_distribution.html

It is in fact good if you can avoid this check as it is on such a low programming level, but it is currently unclear if it is theoretically correct.

[haberg-1]

The theory is easy:

If the high bit of the divisor is set, the quotient with its top word can overflow with at most 2 over the full division that is sought for. Then there are various methods to capture and adjust this overflow. The one in your code captures all 2 overflows but misses the 1 overflows with low probability.

Assuring that the high bit of the divisor is set can be achieved by shifting left both dividend and divisor (if not larger than the divisor) with the same amount, and as it is a multiplication with the same factor, it does not change the quotient, but the remainder, which also gets this factor, must be shifted back.

[danlark1]

All examples seem to use at least a 3-word divisor. You might check your code yourself against binary division, which is easy to implement though slower, and which may have been what you had before. Random numbers can be generated by the example on this page: https://en.cppreference.com/w/cpp/numeric/random/uniform_int_distribution.html

It is in fact good if you can avoid this check as it is on such a low programming level, but it is currently unclear if it is theoretically correct.

I'd politely ask to help finding a counterexample together. So far I was not able to reproduce the issue. One difference can be is that the function we have has a precondition to fit the quotient in 64 bits, it can be the case we don't need to check this overflow because of that