Asymptotically faster divmod and str(long)

[fredrikj]

BPO	3451
Nosy	@tim-one, @mdickinson, @vstinner, @cfbolz, @ericvsmith, @IIIIllllIIIIllllIIIIllllIIIIllllIIIIll, @serhiy-storchaka
Files	div.py: python implementation of Newton division div_timings.txt: Comparison of performance between builtin and Newton-based divmod and str fast_div.py: Fast recursive division diff: division with Burnikel , Ziegler algorithm longobject2.diff longformat.diff: str(n) recursive algorithm longformat_BZ.diff: str(n) and BZ division algorithm

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = None
closed_at = <Date 2021-04-13.15:47:42.644>
created_at = <Date 2008-07-27.01:11:04.477>
labels = ['interpreter-core', 'performance']
title = 'Asymptotically faster divmod and str(long)'
updated_at = <Date 2022-01-30.01:41:10.989>
user = 'https://bugs.python.org/fredrikj'

bugs.python.org fields:

activity = <Date 2022-01-30.01:41:10.989>
actor = 'tim.peters'
assignee = 'none'
closed = True
closed_date = <Date 2021-04-13.15:47:42.644>
closer = 'vstinner'
components = ['Interpreter Core']
creation = <Date 2008-07-27.01:11:04.477>
creator = 'fredrikj'
dependencies = []
files = ['10990', '10991', '11060', '11423', '11773', '13496', '13497']
hgrepos = []
issue_num = 3451
keywords = ['patch']
message_count = 26.0
messages = ['70309', '70313', '70316', '70325', '70326', '70702', '70704', '70719', '70720', '70776', '71010', '72772', '73698', '74679', '84105', '84702', '84704', '103028', '190407', '221270', '390940', '390966', '390968', '390971', '390972', '412121']
nosy_count = 10.0
nosy_names = ['tim.peters', 'mark.dickinson', 'pernici', 'vstinner', 'Carl.Friedrich.Bolz', 'eric.smith', 'fredrikj', 'steve21', 'xuanji', 'serhiy.storchaka']
pr_nums = []
priority = 'normal'
resolution = 'wont fix'
stage = 'resolved'
status = 'closed'
superseder = None
type = 'performance'
url = 'https://bugs.python.org/issue3451'
versions = ['Python 3.2']

[fredrikj]

A few weeks ago, I blogged about taking advantage of Karatsuba
multiplication and Newton's method to divide big integers quickly (some
of you may have read it, as it was posted to Daily Python-URL among
other places):

http://fredrik-j.blogspot.com/2008/07/making-division-in-python-faster.html

To summarize, this method blows the builtin division out of the water
already at ~(2000 digits)/(1000 digits).

The primary catch is that the result of the Newton division may be
slightly wrong (typically 1 ulp). However, a single extra multiplication
and a subtraction at the end allows one to compute a remainder, and
since the remainder must satisfy 0 <= r < q, the error is easily
corrected. From a quick test, the cost of the extra multiplication seems
to move the break-even point with the builtin division up to around
5000/2500 digits.

A pure Python implementation of divmod, with error correction based on
the remainder, can be found in this file:

http://www.dd.chalmers.se/~frejohl/code/libarith/libarith.py

(See the function idivmod)

Of particular note is that fast divmod gives a fast way to do radix
conversion, by recursively splitting the number in half. The function
numeral (see same .py file) does this, e.g:

>>> from time import clock
>>> a = 2**1257787-1
>>> t1=clock(); s1=str(a); t2=clock(); t2-t1
109.08065159119386
>>> t1=clock(); s2=numeral(a); t2=clock(); t2-t1
7.1394886780303182
>>> s1 == s2
True

(This recursive algorithm, by the way, is actually faster than str()
even with the slow builtin divmod.)

Would there be any interest in porting these algorithms to C and using
them in the standard Python long implementation?

There are likely some problems that I have overlooked. A critical review
will be most welcome.

[mdickinson]

Would there be any interest in porting these algorithms to C and using
them in the standard Python long implementation?

Yes, definitely! Though it depends a little bit how much complication
is involved.

A subquadratic algorithm for converting strings to longs might also fit
well here.

Now I'm just waiting for you to propose an implementation of integer
square root :-).

[fredrikj]

Yes, definitely! Though it depends a little bit how much complication
is involved.

Not that much, I think, the pure Python version being just a few lines
of code (plus a few more lines of boilerplate). But I'm not too familiar
with converting Python to C myself, so I may underestimate the
difficulties. It might get a little more complicated if you want to
minimize temporary allocations, for example.

Now I'm just waiting for you to propose an implementation of integer
square root :-).

Yes, there is also code for fast (O(M(n)) integer square roots in
libarith.py. This function would be much less useful overall as a
builtin, though.

[fredrikj]

For your convenience, I have split the division and numeral code off to
a standalone .py file which I'm attaching here.

I also updated the remainder logic to use the default divmod instead of
repeated subtraction, which ensures worst-case performance similar to
the builtin divmod in the very unlikely case that div_newton would
return junk.

[fredrikj]

And here some timings on my laptop.

Both int->str and balanced division become faster somewhere between 1000
and 2000 digits. This is after editing the division benchmark to use
random numbers instead of exact powers of ten (interestingly causing a
bit of slowdown). String conversion might be a little slower for lengths
that aren't close to powers of two.

[mdickinson]

There's also the recursive division algorithm due
to Burnikel and Ziegler; this might be worth a look.
I think it's the same asymptotic complexity (constant
times karatsuba multiplication complexity), but may
turn out to be faster for one reason or another.

I had a Python implementation of this somewhere;
I'll see if I can dig it out.

Assigning this to me so that it doesn't get lost or
forgotten; but note that I don't intend to do
anything about it before 2.6/3.0 final. If anyone else
wants to take it off my hands before then, feel free.

[mdickinson]

See also bpo-1746088.

[mdickinson]

Here's a pure Python implementation of the Burnikel and Ziegler recursive
division algorithm. I've no idea whether it's faster or slower than
Newton, but it might be worth a look. It depends heavily on bit
operations, which ought to be much faster when coded in C. (Some of the
shifts would be completely unnecessary---replaced by changes in indexing
instead.)

The original paper describing the algorithm is available here:

http://cr.yp.to/bib/1998/burnikel.ps

[mdickinson]

Oops. Wrong file. Here's the right one.

[fredrikj]

Indeed, that seems to be much faster. In the mean time, do you mind if I
steal the code? :-)

[mdickinson]

Indeed, that seems to be much faster. In the mean time, do you mind if I
steal the code? :-)

Not at all!

I guess constant factors could well appear and/or disappear when
recoding in C; it might well be worth trying to implement both
methods to see which is faster.

There are some semi-obvious tricks available that might speed up the
recursive version. For one thing, all shifts should be in multiples
of 15 bits (because longs are stored in base 2**15). For another, it
ought to be possible to avoid the single-bit normalization shifts
every time the number of bits ('n') is odd---one can do a single
shift at the beginning of the calculation instead.

[pernici]

I have translated in C the algorithm for long division by
Burnikel and Ziegler (BZ), using the Python version fast_div.py
and the suggestions by Mark.

Here is a benchmark for divmod(p. q), p = 7**np, q = 5**nq
digits = q_digits = p_digits/2; the time taken by Python division
is normalized to 1
tests on Debian, BZ with Python-2.6b3
Pentium 1.60GHz Athlon XP 2600+ Athlon 64 3800+
digits BZ fast_div BZ fast_div BZ fast_div
500 1.01 1.27 1.00 1.18 1.00 1.21
700 0.88 1.22 0.76 1.08 0.81 1.14
1000 0.82 1.17 0.72 1.04 0.76 1.10
2000 0.66 0.85 0.55 0.73 0.59 0.78
4000 0.51 0.62 0.43 0.52 0.45 0.56
10000 0.32 0.38 0.31 0.37 0.33 0.39
20000 0.24 0.25 0.23 0.25 0.25 0.27
100000 0.14 0.14 0.11 0.11 0.12 0.12

BZ starts being faster than Python division around 2000 bits, apart
from two cases:
for q with less than 4000 bits and p much smaller than q**2
x_divrem is faster than divmod_pos;
for p very much larger than q**2 x_divrem becomes
again faster.
I put a bound in long_divrem to fall back to x_divrem in those cases,
based on tests on my laptop Pentium 1.60GHz.

The treatment of exceptions is very rudimentary;
I use a simple macro STUB_EXC to return NULL,
but without releasing memory.
Help for doing this properly is welcome.

Please find attached the diff to the longobject.c file in Python-2.6b3 .

Mario