Infinite loop in Burnikel-Ziegler division algorithm in _pylong.py

[Julien-Livet]

Bug report

Bug description:

I want to transpose the code in C++ but I discovered a problem in Python code.

Here is the sample test:

a = 561453*2**64 + 13205*2**32 + 1564
b = 1698*2**32 + 721
q = 330*2**32 + 2815252282
r = 414*2**32 + 1722637250
print("a", a)
print("b", b)
print("q", q)
print("r", r)
q_, r_ = int_divmod(a, b)
print("q_", q_)
print("r_", r_)

If I comment lines 435 and 436 for testing:

if a.bit_length() - n <= _DIV_LIMIT:
        return divmod(a, b)

I obtained an infinite loop and I don't understand why.
Is the algorithm work on an existing test?

CPython versions tested on:

3.12

Operating systems tested on:

Windows

[ZeroIntensity]

cc @skirpichev

[vstinner]

cc @sobolevn @tim-one @serhiy-storchaka

[eendebakpt]

As far as I understand it (but it would be nice to have this confirmed):

The lines you are commenting out are from _div2n1n which is defined in _pylong.py. That method corresponds to DivTwoDigitsByOne from the paper Fast Recursive Division, Burnikel-Ziegler, 1998. That paper states some conditions for the method and by commenting out the lines you are violating the "B has to be normalized". When I try the algorithm the _div2n1n recurses to a=1, b=1, n=1 where this no longer holds.

[skirpichev]

if a.bit_length() - n <= _DIV_LIMIT:
        return divmod(a, b)

This is actually just a base of recursion (see e.g. RecursiveDivision() function on page 7). Obviously, you should add something to terminate recursion. Either switch to school division when n is "small enough" (as in the current variant) or implement a special case when just two digits divided by one (I would guess this will be much more slow at outcome).

[tim-one]

Yes, when you comment out the base cases for a recursive algorithm, you're going to get an infinite loop. This was self-inflicted, and in no way a bug in Python, so I'm closing this.

What may be more surprising is why _DIV_LIMIT is so large. That's a matter of tuning. It was picked to maximize speed under CPython at the time, where the overheads of running the fancy recursive Python code are high enough that CPython's own "naïve division", coded in C, runs faster until the number of result bits needed is quite high. Under the PyPy implementation of Python, _DIV_LIMIT can be reduced substantially to advantage. Under C++, doubtless even more so.

But, regardless, you do need an entirely different algorithm to fall back on to handle base cases. Both for correctness and for speed.

To get the pragmatics across quickly, suppose you're dividing a million-bit integer into an int with a million and one bits. There are very few bits in the quotient, and essentially nothing could be faster than just trying plain subtraction a few times. The test in question is ensuring the quotient needs at least about 4 thousand bits before the overheads of the fancy method are overcome.