Knuth's Algorithm D

also: long division, multiprecision division

Long division on multi-limb integers. The loop is straightforward; the quotient digit is not. Normalizing so the divisor's top limb is at least half the base bounds the estimate's error at two, which a short correction loop removes, with a rare add-back when the multiply-and-subtract goes negative.

Knuth gives it as Algorithm D in <i>Seminumerical Algorithms</i> 4.3.1, and the reason it takes a page rather than a line is the quotient digit. Guessing it from the leading limbs is cheap; the guess can be wrong, and the algorithm is mostly about bounding and repairing that error.

Normalize first, by shifting both operands left until the divisor's top limb is at least half the base. That bounds the error in the estimate

q^= min (⌊ \frac{u _{j + n} β + u _{j + n - 1}}{v _{n - 1}} ⌋, β - 1)

to at most two too large, so a short correction loop fixes it. After the multiply-and-subtract step goes negative, a single add-back restores it, and that case is rare enough that Knuth notes it happens with probability about $2/ β$ .

The correction, which is the part that is easy to get wrong and hard to notice.

(* qhat is at most 2 too large after normalization *)
while !qhat >= base
      || !qhat * v.(n - 2) > (!rhat * base) + u.(j + n - 2) do
  decr qhat;
  rhat := !rhat + v.(n - 1);
  if !rhat >= base then qhat := 0 (* force loop exit *)
done

Two divisions from the Part 1 bignum, cross-checked against Python. The second is the add-back case.

10^40 / 123456789012345678901
  q = 81000000729000006634
  r = 6661773269766170766
2^127 / (2^64 + 1)
  q = 9223372036854775807
  r = 9223372036854775809
both satisfy q*b + r = a

Everything above it in a CAS depends on this being exactly right. A rational is only in lowest terms because a GCD said so, and a GCD is a chain of divisions.