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)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 Multiprecision integers in hardware Reducing long division to medium division Reducing medium division to small division Normalisation The bug 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 Stronger bounds Doubling the quotient limbs Division by a constant How I got herePreparing 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 scratchMultiprecision integers in hardwareCryptographic 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: x=(xn−1,…,x0)b=∑i=0n−1xibi,0≤xi<b. 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 subc: x, y -> d, borrow # d = (x−y) mod b, borrow = 1 iff underflow mul: x, y -> (hi, lo) # x·y = hi·b + lo div: (hi, lo), y -> (q1, q0), r # hi·b + lo = q·y + r, 0 ≤ r < y, q = q1·b + q0 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, although some architectures treat this case differently as we will see later. Reducing long division to medium divisionOur task is to divide a multiprecision integer u by v, that is, find integers q and r such that u=q·v+r and 0≤r<v. Assume without loss of generality that the divisor v is an n-limb integer with non-zero top limb vn−1, and pad the dividend u with leading zeros until it has strictly more limbs than v; write n+m+1 for the number of limbs of u. We have another requirement which will prove natural in what follows: the highest n limbs of u, read as an n-limb integer, must be strictly smaller than v, (un+m,…,um+1)b<(vn−1,…,v0)b,equivalently⌊ubm+1⌋<v.If that is not already the case, appending a further zero to u will guarantee it. This assumption pins down the size of the quotient to m+1 limbs: ⌊ubm+1⌋<v⟺u<vbm+1⟺q=⌊uv⌋<bm+1.With the padding in place every operand has a fixed shape. Writing q and r for the quotient and remainder of u by v, we have: u = (u_{n+m}, u_{n+m-1}, ..., u_0)_b v = (v_{n-1}, v_{n-2}, ..., v_0)_b 0 < v_{n-1} q = (q_{m}, q_{m-1}, ..., q_0)_b u = q·v + r r = (r_{n-1}, r_{n-2}, ..., r_0)_b 0 ≤ r < v In this exposition we will use three different division algorithms, at three levels. The one we want is the n+m+1 by n limb, or the "long" n+m+1/n division. The one we have is the "short" 2/1 division instruction. To bridge them we use the "medium" n+1/n division. A natural way to compute the limbs of q is by going from the top down.3 The top limb is4 qm=⌊⌊u/v⌋bm⌋=⌊⌊u/bm⌋v⌋=⌊(un+m,…,um)b(vn−1,…,v0)b⌋.The numerator ⌊u/bm⌋ is simply the top n+1 limbs of u, so the top limb of the quotient qm is itself a quotient of an n+1/n division. And it fits in a single limb due to the hypothesis on the top n limbs of u being less than v (un+m,…,um+1)b<v⟺(un+m,…,um)b<bv,hence qm=⌊(un+m,…,um)bv⌋<b. Upon computing qm, we proceed by subtracting the multiple bmqmv of v from u, and continuing the algorithm with the updated u. This operation of updating u is equivalent to replacing the top n+1 limbs of u by Rm, the remainder corresponding to qm, which we know to fit in n limbs. Therefore the updated u will have one limb less, and we continue the algorithm to compute the remaining m limbs of q. While it might not be straightforward that this folding technique computes the correct answer, it can easily be deduced from observing the following algorithm and noting that the two invariants are satisfied at each entrance and exit of the loop: Algorithm 1: Long Division Input: u = (u_{n+m}, ..., u_0)_b, padded so top n limbs < v v = (v_{n-1}, ..., v_0)_b, v_{n-1} > 0 Start: r = u, q = 0 // r starts as u and shrinks to final remainder Loop: for k = m down to 0: u' = (r_{k+n}, ..., r_k)_b // top n+1 limbs of r at position k (q_k, R_k) = ⌊u'/v⌋, u' mod v // an n+1/n division q += q_k·b^k // k'th limb of the quotient r -= q_k·v·b^k // replaces (r_{k+n},...,r_k) by (0, R_k) Return: (q, r) Invariants: 1. u = q·v + r 2. top n limbs of r < v // from limb k+1 to n+k Running invariant 2 until k=0 leaves r<v, and so the pair (q,r) is exactly the quotient and the remainder. Long Division exampler is initialised to u, already shown. Each iteration: divide the highlighted top four words of r by v=314 (quotient digit to q); take them modulo 314 into r_k (r3,r2,r1,r0); clear those words of r; then move r_k up into r. After four steps q=5599 and r=207.Long Division exampleurvqrkr3r2r1r031455991758293/%/%/%/%183813202307b = 10 Step 1 / 15 Click Next / Prev or use ← →. Play runs through every step; Space toggles it. So an n+m+1/n long division costs m+1 medium n+1/n divisions, one per quotient limb. Reducing medium division to small divisionThe previous section left us with the medium division: u = (u_n, u_{n-1}, ..., u_0)_b, v = (v_{n-1}, ..., v_0)_b, (u_n, ..., u_1) < v, i.e., u/b < v, i.e., u < vb The division of u by v returns a quotient 0≤q=⌊uv⌋<b and remainder 0≤r<v. A natural step is to approximate the n+1/n division by means of a 2/1 division of the top limbs, which can be computed with the div instruction. Re-write the parameters as u = u''·b^{n-1} + u' u'' = (u_n, u_{n-1})_b 0 ≤ u'' < b², 0 ≤ u' < b^{n-1}, v = v''·b^{n-1} + v' v'' = v_{n-1} 0 < v'' < b, 0 ≤ v' < b^{n-1}, and set q^=⌊u″v″⌋ computed by a single div((u_{n},u_{n-1}), v_{n-1}). Note that q^ is the full two-limb quotient (q1,q0)b of that instruction. Our constraint bounds u against v·b, but not u″ against v″·b, so nothing stops q^ from exceeding a single limb. How good of a guess is q^? Can it overshoot badly? Can it undershoot? Both theorems below lean on the following simple fact about floors: for d=⌊m/n⌋ we have dn≤m≤(d+1)n−1. Theorem A (q≤q^): No undershoot q≤uv≤u″bn−1+u′v″bn−1≤(u″+1)bn−1−1v″bn−1<u″+1v″≤q^+1.The second step uses v≥v″bn−1, the third uses u′≤bn−1−1, and the last is the floor fact, in the form u″+1≤(q^+1)v″. Since q<q^+1 and both are integers, q≤q^. Theorem B' (q^<q+1+b/v″): Bounded overshoot q^≤u″v″=u″bn−1v″bn−1≤uv″bn−1≤(q+1)v−1v″bn−1<q+1+bv″.The third step is the floor fact for q=⌊u/v⌋. The last step expands v=v″bn−1+v′ and uses