A long division story – Novak Kaluđerović - Kolja

A long division story

13 Aug, 2026

How I received a theorem in Knuth's "The Art of Computer Programming" after finding a decades-old bug in Algorithm D (+ a "bug" in llvm)<br>I was implementing Algorithm D, the well-known long division algorithm from Knuth's "The Art of Computer Programming", and I stumbled upon an issue that I couldn't let go. The correctness of the algorithm relied on Theorem B, and its proof bugged me. It felt unnatural, it took a very convoluted path to proving a simple statement, and it isolated a special case which was not a corner case and that seemed unrelated to the problem at hand. There was something odd about it, so I tried to prove the theorem myself, and I failed. However, the failure handed me a counterexample to Algorithm D that had passed as correct for decades, and with it a theorem on the correctness of the algorithm carrying my name.

In this post I'll give some background, then cover long division from scratch for those who want it, share my thoughts on how the bug came to be and how it stayed hidden so long, and finish with a preview of more modern ways to implement long division. In the process of writing this blog I also found a "bug" in this algorithm's implementation in llvm, and I will expand on that too. If you're only interested in the bug you can jump directly to The bug.

Contents

How I got here

Long division from scratch<br>Multiprecision integers in hardware

Reducing long division to medium division

Reducing medium division to small division<br>Normalisation

The bug<br>How did it stay hidden for decades

Can it be exploited

The llvm "bug"

AI didn't find it

The check

A little trit more<br>Stronger bounds

Doubling the quotient limbs

Division by a constant

How I got here<br>Preparing for an interview, I decided to do a small project: build a little library for arithmetic over prime fields. This meant fixed-size multiprecision integers, arithmetic operations, some field operations, constant-time, constant memory access, all in all a starting point for modern cryptographic protocols.

As I worked on implementing it, the process turned into a game with one rule: avoid division at all costs. You can get almost all the way there, and to the best of my knowledge cryptographic libraries never execute the division instruction at runtime. Whenever a divide is needed, it is generally substituted by a multiplication followed by a bit of shuffling.1

Why do we go so far out of our way to avoid division? Well, multiplication is very simple, in fact it can be thought of as an axiom of the natural numbers. Division is far more complicated. Firstly, it's not everywhere defined: we can't divide by zero. But we can't divide 5 by 2 either! What we actually have is "division with remainder", a more complicated operation which returns two answers: the quotient and the remainder, the smallest non-negative difference between the dividend and a multiple of the divisor. The issue hides in what "smallest" exactly means, why we choose this particular definition, and why the notion of size enters the picture at all. One may choose differently, say zero-centred remainders. But we could go further and choose a different size function, which gives rise to a different division algorithm altogether.2

With all that in mind it's no wonder that the theoretical complication transfers into practice. A multiply instruction costs a cycle or two on modern machines and fully pipelines, while a divide can cost up to twenty cycles and usually doesn't pipeline.

Eventually the only gap left in the multiprecision implementation was the multiprecision division algorithm. So I did the obvious thing and sat down to implement long division, and for reference I used Donald Knuth's "The Art of Computer Programming" Vol. II, Third Edition, Algorithm 4.3.1D.

Let's look at how long division actually works.

Long division from scratch<br>Multiprecision integers in hardware<br>Cryptographic integers run to hundreds or thousands of bits, well past a single register, so we store them in base b, one limb per machine word:<br>x=(xn−1,…,x0)b=∑i=0n−1xibi,0≤xib.

The main building blocks of multiprecision arithmetic algorithms are the four primitive instructions that operate over single/double limbs:

addc: x, y -> s, carry # s = (x+y) mod b, carry = 1 iff overflow<br>subc: x, y -> d, borrow # d = (x−y) mod b, borrow = 1 iff underflow<br>mul: x, y -> (hi, lo) # x·y = hi·b + lo<br>div: (hi, lo), y -> (q1, q0), r # hi·b + lo = q·y + r, 0 ≤ r

The first three are unremarkable, but the division stands as the odd one out. While multiplying two single-limb multiplicands always returns a two-limb product, the quotient of a two-limb dividend over a one-limb divisor does not always fit in a single limb, so we use a two-limb quotient. In addition to that, the remainder shows up as a necessary byproduct. And in division by zero we assume undefined behaviour, i.e. that q and r may take any value,...