Een verhaal over staartdeling

Hoe ik een theorema kreeg in Knuth's "The Art of Computer Programming" nadat ik een decennia oude bug in Algoritme D vond (plus een "bug" in llvm).

Ik was Algoritme D aan het implementeren, het bekende algoritme voor staartdeling uit Knuth's "The Art of Computer Programming", toen ik stuitte op een probleem dat ik niet kon loslaten. De correctheid van het algoritme vertrouwde op Theorema B, en het bewijs daarvan stoorde me. Het voelde onnatuurlijk; het nam een zeer omslachtige weg om een eenvoudige bewering te bewijzen en isoleerde een speciaal geval dat geen randgeval was en dat niet gerelateerd leek aan het probleem in kwestie. Er was iets vreemds aan, dus ik probeerde het theorema zelf te bewijzen, en dat mislukte. Echter, die mislukking leverde me een tegenvoorbeeld op voor Algoritme D dat decennia lang als correct was beschouwd, en daarmee een theorema over de correctheid van het algoritme dat mijn naam draagt.

In dit artikel geef ik wat achtergrondinformatie, behandel ik staartdeling vanaf nul voor wie dat wil, deel ik mijn gedachten over hoe de bug is ontstaan en waarom deze zo lang verborgen is gebleven, en sluit ik af met een vooruitblik op modernere manieren om staartdeling te implementeren. Tijdens het schrijven van deze blog vond ik ook een "bug" in de implementatie van dit algoritme in llvm, waar ik verder op in zal gaan.

Hoe ik hier kwam

Ter voorbereiding op een sollicitatiegesprek besloot ik een klein project te doen: een bibliotheek bouwen voor rekenkunde over priemvelden. Dit betekende multiprecisie-geheelgetallen met een vaste grootte, rekenkundige bewerkingen, enkele veldoperaties, constante tijd en constante geheugenbenadering — kortom, een startpunt voor moderne cryptografische protocollen.

Terwijl ik aan de implementatie werkte, veranderde het proces in een spel met één regel: vermijd deling te allen tijde. Je kunt bijna het hele proces zo inrichten, en voor zover ik weet, voeren cryptografische bibliotheken tijdens runtime nooit de delingsinstructie uit. Wanneer er een deling nodig is, wordt deze over het algemeen vervangen door een vermenigvuldiging gevolgd door wat bit-verschuivingen.[^1]

Waarom doen we zoveel moeite om deling te vermijden? Vermenigvuldiging is immers zeer eenvoudig; het kan zelfs worden beschouwd als een axioma van de natuurlijke getallen. Deling is veel ingewikkelder. Ten eerste is het niet overal gedefinieerd: we kunnen niet delen door nul. Maar we kunnen 5 ook niet delen door 2! Wat we eigenlijk hebben is "deling met rest", een complexere operatie die twee antwoorden teruggeeft: het quotiënt en de rest (het kleinste niet-negatieve verschil tussen het dividend en een veelvoud van de divisor). Het probleem schuilt in wat "kleinste" precies betekent, waarom we deze specifieke definitie kiezen, en waarom het begrip 'grootte' überhaupt een rol speelt. Men zou anders kunnen kiezen, bijvoorbeeld nul-gecentreerde resten. Maar we zouden nog verder kunnen gaan en een andere groottefunctie kiezen, wat leidt tot een geheel ander delingsalgoritme.[^2]

Met dit alles in het achterhoofd is het niet vreemd dat de theoretische complicatie zich vertaalt naar de praktijk. Een vermenigvuldigingsinstructie kost op moderne machines een cyclus of twee en is volledig pipelinebaar, terwijl een deling tot twintig cycli kan kosten en meestal niet pipelinebaar is.

Uiteindelijk was de enige openstaande post in de multiprecisie-implementatie het algoritme voor multiprecisie-deling. Dus deed ik het voor de hand liggende en ging ik aan de slag met het implementeren van staartdeling. Als referentie gebruikte ik Donald Knuth's "The Art of Computer Programming" Vol. II, Derde Editie, Algoritme 4.3.1D.

Staartdeling vanaf nul

Multiprecisie-geheelgetallen in hardware

Cryptografische geheltallen lopen op tot honderden of duizenden bits, wat ver voorbij een enkel register gaat. Daarom slaan we ze op in basis $b$, met één limb per machinewoord:

$x = (x{n-1}, \dots, x0)b = \sum{i=0}^{n-1} xi b^i, \quad 0 \le xi < b$

De belangrijkste bouwstenen van multiprecisie-rekenalgoritmen zijn de vier primitieve instructies die werken op enkele of dubbele limbs:

• addc: $x, y \rightarrow s, \text{carry}$ (waarbij $s = (x+y) \pmod b$, $\text{carry} = 1$ bij overflow)
• subc: $x, y \rightarrow d, \text{borrow}$ (waarbij $d = (x-y) \pmod b$, $\text{borrow} = 1$ bij underflow)
• mul: $x, y \rightarrow (hi, lo)$ (waarbij $x \cdot y = hi \cdot b + lo$)
• div: $(hi, lo), y \rightarrow (q1, q0), r$ (waarbij $hi \cdot b + lo = q \cdot y + r, \quad 0 \le r < y, \quad q = q1 \cdot b + q0$)

De eerste drie zijn onopvallend, maar de deling is de vreemde eend in de bijt. Terwijl het vermenigvuldigen van twee single-limb multiplicanden altijd een two-limb product oplevert, past het quotiënt van een two-limb dividend gedeeld door een one-limb divisor niet altijd in een enkele limb; daarom gebruiken we een two-limb quotiënt. Daarnaast verschijnt de rest als een noodzakelijk bijproduct. Bij deling door nul gaan we uit van ongedefinieerd gedrag, hoewel sommige architecturen dit anders behandelen.

Lange deling reduceren tot middelgrote deling

Onze taak is om een multiprecisie-geheelgetal $u$ te delen door $v$, oftewel geheltallen $q$ en $r$ vinden zodanig dat $u = q \cdot v + r$ en $0 \le r < v$. Stel zonder verlies van algemeenheid dat de divisor $v$ een $n$-limb geheelgetal is met een niet-nul top-limb $v_{n-1}$, en vul het dividend $u$ aan met leidende nullen totdat het strikt meer limbs heeft dan $v$; we noteren $n+m+1$ voor het aantal limbs van $u$.

We hebben een aanvullende vereiste: de hoogste $n$ limbs van $u$, gelezen als een $n$-limb geheelgetal, moeten strikt kleiner zijn dan $v$:

$(u{n+m}, \dots, u{m+1})b < (v{n-1}, \dots, v0)b, \quad \text{equivalent aan } \lfloor u/b^{m+1} \rfloor < v$

Indien dit niet het geval is, garandeert het toevoegen van een extra nul aan $u$ dit. Deze aanname bepaalt de grootte van het quotiënt op $m+1$ limbs:

$\lfloor u/b^{m+1} \rfloor < v \iff u < v b^{m+1} \iff q = \lfloor u/v \rfloor < b^{m+1}$

Met de padding op zijn plek heeft elke operand een vaste vorm:

• $u = (u{n+m}, u{n+m-1}, \dots, u0)b$
• $v = (v{n-1}, v{n-2}, \dots, v0)b \quad (0 < v_{n-1})$
• $q = (qm, q{m-1}, \dots, q0)b \quad (u = q \cdot v + r)$
• $r = (r{n-1}, r{n-2}, \dots, r0)b \quad (0 \le r < v)$

We gebruiken drie verschillende delingsalgoritmen op drie niveaus. Het doel is de "lange" $(n+m+1)/n$ deling. Wat we hardwarematig hebben is de "korte" $2/1$ delingsinstructie. Om deze te overbruggen, gebruiken we de "middelgrote" $(n+1)/n$ deling.

Een natuurlijke manier om de limbs van $q$ te berekenen is van boven naar beneden.[^3] De top-limb is:

$qm = \lfloor \lfloor u/v \rfloor / b^m \rfloor = \lfloor \lfloor u/b^m \rfloor / v \rfloor = \lfloor (u{n+m}, \dots, um)b / (v{n-1}, \dots, v0)_b \rfloor$

De teller $\lfloor u/b^m \rfloor$ is simpelweg de top $n+1$ limbs van $u$. De top-limb van het quotiënt $qm$ is dus zelf een quotiënt van een $(n+1)/n$ deling. Dankzij de hypothese dat de top $n$ limbs van $u$ kleiner zijn dan $v$, past $qm$ in een enkele limb:

$(u{n+m}, \dots, u{m+1})b < v \iff (u{n+m}, \dots, um)b < bv$, dus $qm = \lfloor (u{n+m}, \dots, um)b / v \rfloor < b$.

Na het berekenen van $qm$ trekken we het veelvoud $b^m qm v$ af van $u$ en gaan we verder met het bijgewerkte $u$. Deze operatie is equivalent aan het vervangen van de top $n+1$ limbs van $u$ door $Rm$, de rest die hoort bij $qm$, waarvan we weten dat deze in $n$ limbs past.

Algoritme 1: Lange Deling

• Input: $u$ (gepadt zodat top $n$ limbs $< v$), $v$ (met $v_{n-1} > 0$)
• Start: $r = u, q = 0$
• Loop: Voor $k$ van $m$ naar 0:

1. $u' = (r{k+n}, \dots, rk)_b$ (top $n+1$ limbs van $r$ op positie $k$)
2. $(qk, Rk) = \lfloor u'/v \rfloor, u' \pmod v$ (een $(n+1)/n$ deling)
3. $q += q_k \cdot b^k$
4. $r -= qk \cdot v \cdot b^k$ (vervangt $(r{k+n}, \dots, rk)$ door $(0, Rk)$)

• Return: $(q, r)$

Een lange deling kost dus $m+1$ middelgrote $(n+1)/n$ delingen.

Middelgrote deling reduceren tot kleine deling

De middelgrote deling is als volgt: $u = (un, u{n-1}, \dots, u0)b$, met $\lfloor u/b \rfloor < v$ (dus $u < vb$). De deling van $u$ door $v$ geeft een quotiënt $0 \le q = \lfloor u/v \rfloor < b$ en rest $0 \le r < v$.

Een natuurlijke stap is om de $(n+1)/n$ deling te benaderen via een $2/1$ deling van de top-limbs. We herschrijven de parameters als: $u = u'' \cdot b^{n-1} + u', \quad u'' = (un, u{n-1})b \quad (0 \le u'' < b^2, 0 \le u' < b^{n-1})$ $v = v'' \cdot b^{n-1} + v', \quad v'' = v{n-1} \quad (0 < v'' < b, 0 \le v' < b^{n-1})$

We stellen $\hat{q} = \lfloor u''/v'' \rfloor$, berekend met een enkele div((un, u{n-1}), v_{n-1}) instructie.

Theorema A (geen onderschatting): $q \le \hat{q}$. (Bewijs: $q \le u/v \le u'' b^{n-1}/v'' b^{n-1} < (u''+1)/v'' \le \hat{q}+1$. Omdat $q < \hat{q}+1$, geldt $q \le \hat{q}$).

Theorema B' (begrensde overschatting): $\hat{q} < q + 1 + b/v''$.

Gecombineerd hebben we de grens: $q \le \hat{q} < q + 1 + b/v''$. Voor kleine $v''$ is deze grens niet erg nauwkeurig. Dit lossen we op via normalisatie.

Normalisatie

Het quotiënt blijft ongewijzigd wanneer beide operanden met dezelfde factor worden geschaald: $q = \lfloor u/v \rfloor = \lfloor uf/vf \rfloor$. We zoeken een factor $f$ die $v''$ groot maakt. De keuze $f = \lfloor b/(v''+1) \rfloor$ werkt. De nieuwe top-limb voldoet dan aan $v'' \ge \lfloor b/2 \rfloor$. Een dergelijke $v$ noemen we genormaliseerd.

Met $v'' \ge \lfloor b/2 \rfloor$ wordt Theorema B': Theorema B: $\hat{q} \le q + 3$.

De procedure voor $(n+1)/n$ deling is nu: bereken $\hat{q}$, test de gok door $u - \hat{q}v$ te vormen. Als de aftrekking underflowt, voegen we $v$ opnieuw toe en verlagen we $\hat{q}$ met één. We herhalen deze correctie maximaal drie keer.

Algoritme 2: Middelgrote Deling

• Input: $u$ (met $\lfloor u/b \rfloor < v$), $v$ (genormaliseerd: $v_{n-1} \ge \lfloor b/2 \rfloor$)
• Trial: $(\hat{q}, \hat{r}) = \text{div}((un, u{n-1}), v_{n-1})$
• Fix: $r = u - \hat{q} \cdot v$
• Correctie:
• Als $r < 0$: $\hat{q} -= 1, r += v$
• Als $r < 0$: $\hat{q} -= 1, r += v$
• Als $r < 0$: $\hat{q} -= 1, r += v$
• Return: $(\hat{q}, r)$

De bug

In Algoritme D (stap D3) stelt Knuth een methode voor om de trial-quotiënt te verbeteren voordat de fix-fase begint. Dit heeft twee voordelen:

1. $\hat{q}$ wordt aangescherpt naar de grens $\hat{q} \le q+1$. Dit betekent dat er maximaal één correctiestap nodig is in plaats van drie.
2. Knuth stelt dat $\hat{q}$ aan het einde van stap D3 in een enkele limb past, waardoor de daaropvolgende vermenigvuldiging $\hat{q} \cdot v$ een $1 \times n$ vermenigvuldiging is in plaats van $2 \times n$, wat de prestaties verbetert.

Algoritme 3: Middelgrote Deling (Knuth)

• D3:

1. $(\hat{q}, \hat{r}) = \text{div}((un, u{n-1}), v_{n-1})$
2. Als $\hat{q} \ge b$ of $\hat{q} \cdot v{n-2} > b \cdot \hat{r} + u{n-2}$:

$\hat{q} -= 1; \hat{r} += v_{n-1}$;

1. Als $(\hat{r} < b)$ en ($\hat{q} \ge b$ of $\hat{q} \cdot v{n-2} > b \cdot \hat{r} + u{n-2}$):

$\hat{q} -= 1; \hat{r} += v_{n-1}$;

• D4: $r = u - \hat{q} \cdot v$ (nu een $1 \times n$ vermenigvuldiging)
• D5: Als $r < 0$: $\hat{q} -= 1, r += v$

Een oplettende lezer zal een probleem zien. In Theorema B hadden we de grens $\hat{q} \le q+3$. Stap D3 voert echter maximaal twee correcties uit, waarna wordt verwacht dat $\hat{q}$ in één limb past. In het ergste geval ($q = b-1$ en $\hat{q} = q+3 = b+2$) zullen twee correcties niet genoeg zijn om $\hat{q}$ in een enkele limb te dwingen.

Dit scenario gebeurt daadwerkelijk. Het kleinste voorbeeld waarin Algoritme D een fout antwoord geeft is:

• $b = 3$
• $u = (1,2,0,0)_3 = 45$
• $v = (1,2,1)_3 = 16$
• $q = \lfloor 45/16 \rfloor = 2$
• $\hat{q} = \lfloor (1,2)3 / (1)3 \rfloor = \lfloor 5/1 \rfloor = 5$

Na twee correcties in stap D3 is $\hat{q}$ nog steeds een two-limb waarde: $\hat{q} = (1,0)_3$. De $1 \times n$ vermenigvuldiging ziet vanaf dat punt alleen de lage limb van $\hat{q}$ (die nul is), trekt dus niets af, en de fouten planten zich voort.

Hoe is dit decennia lang verborgen gebleven?

Opvallend genoeg vertoont Program D, Knuth's implementatie in MIX6 assembly, deze fout niet. Het programma is correct, maar het algoritme is fout. De reden ligt in de implementatie van de delingsinstructie.

Er zijn minstens twee manieren om een delingsinstructie te implementeren:

1. divarm: $(hi, lo), (y) \rightarrow (q1, q_0), r$ (geeft het volledige two-limb quotiënt terug).
2. divx86: $(hi, lo), (y) \rightarrow (q0), r, F$ (geeft een single-limb quotiënt terug en zet een vlag bij overflow).

Knuth's MIX-machine gebruikte de verzadigde quotiënt $\hat{q}_{x86}$. Voor deze waarde waren de theorema's in de eerste en tweede editie van het boek correct.

In de jaren '90 introduceerde Knuth de MMIX-machine, wat leidde tot een revisie van de boeken. Op 28 september 1995 werd Algoritme D aangepast. De wijziging verving $\hat{q}{x86}$ door $\hat{q}{arm}$, maar de Theorema's A en B werden niet aangepast; zij bewezen nog steeds alleen de grens voor $\hat{q}{x86}$. Voor de meeste inputs was dit geen probleem, maar in het degenerate geval $\hat{q}{arm} \ge b$ was het gat $\hat{q}{arm} - q$ technisch gezien onbegrensd. De bug treedt alleen op wanneer $b$ oneven is, $v{n-1} = (b-1)/2$, $un = (b-1)/2$ en $u{n-1} = b-1$.

Kan dit worden misbruikt?

Niet echt, tenzij je een Setun (een Sovjet-ternaire computer) gebruikt. Omdat de fout een oneven basis vereist, kan deze niet voorkomen op moderne machines (die $b=2^64$ gebruiken). Hoewel de correctie nodig was voor een volledig bewijs, werkte het algoritme in de praktijk correct.

Mogelijke use-cases voor oneven limbs zijn p-adische getallen voor oneven $p$, of radix-conversie. Bijvoorbeeld: het converteren van een getal van basis 3 naar hexadecimaal. De deling $(1200)3 / (121)3$ is precies het tegenvoorbeeld waarin de bug optreedt.

De llvm-"bug"

Bij het zoeken naar implementaties van Algoritme D merkte ik op dat de meeste implementaties een while-loop gebruiken voor stap D3, wat de bug onbewust oplost. Eén implementatie in de arbitrary-precision integer library van llvm (APInt.cpp) gebruikte echter twee if-statements.

Bij nadere inspectie bleek dat llvm in stap D3 controleerde op qp == b, terwijl de gecorrigeerde tekst van Knuth qp >= b voorschrijft. Dit was een typo in de oorspronkelijke TAOCP Vol. II, gecorrigeerd in 2005. De code van llvm behield de oude versie.

Omdat $b=2^{32}$ even is, is dit technisch gezien geen bug die tot foute resultaten leidt, omdat het geval $\hat{q}=b+1$ wordt opgevangen door de andere check ($\hat{q} \cdot v{n-2} > b \cdot \hat{r} + u{n-2}$). Toch heb ik een PR geopend om de check in lijn te brengen met de gecorrigeerde tekst.

AI vond het niet

Ik gaf Claude 3.5 een PDF van de relevante pagina's uit TAOCP en vroeg hem om een fout of bug te vinden. De AI besteedde veel tijd aan het verifiëren van de complexiteitsanalyse van Program D, maar vond de bug in het algoritme niet. Pas nadat ik hem specifiek wees op de grootte van $\hat{q}$ in de overgang van D3 naar D4, vond hij het, maar pas nadat hij de nieuwste errata van Knuth's website had opgehaald.

De cheque

Knuth beloont elke gevonden fout in TAOCP met één hexadecimale dollar (0x$1.00) gestort bij de Bank of San Serriffe. Ik ontving de cheque en een brief met handgeschreven notities van Knuth. Hij schreef: "Ik ben vooral blij met deze correctie, want ik denk dat lezers van TAOCP Vol 2 Algoritme 4.3.1 D vaker bekijken dan welk ander algoritme dan ook!"

De correctie is ingediend op 14 mei 2026 en gepubliceerd in de errata op 9 juni 2026.

Nog een klein beetje extra

Hier is een kort overzicht van drie methoden om staartdeling te versnellen:

Sterkere grenzen

Door de primitieve bewerking te veranderen van een $2/1$ deling naar een $(k+1)/k$ deling (voor $k \ge 2$), kunnen we sterkere grenzen verkrijgen. Voor $k \ge 2$ geldt $q \le \hat{q} \le q+1$, waardoor normalisatie niet meer nodig is.

Het verdubbelen van de quotiënt-limbs

We kunnen $4/2$ trial-delingen uitvoeren om twee limbs van het quotiënt per iteratie te verkrijgen. Dit halveert het aantal passes over $u$. Dit kan worden uitgebreid naar $8/4$, $16/8$, etc., wat uiteindelijk leidt tot het recursieve delingsalgoritme van Burnikel-Ziegler.

Deling door een constante

Aangezien de divisor $v$ constant blijft tijdens de gehele berekening, kan deling door $v$ worden vervangen door vermenigvuldiging met een vooraf berekende reciproke waarde, gevolgd door een shift. Dit idee, oorspronkelijk van Granlund en Montgomery, wordt gebruikt in GMP.

---

[^1]: In cryptografische bibliotheken wordt deling vaak vervangen door vermenigvuldiging met een modulair inversum. [^2]: Verschillende groottefuncties leiden tot verschillende delingsalgoritmen. [^3]: Schoolboek-vermenigvuldiging gaat van laag naar hoog; deling gaat van hoog naar laag. [^4]: Zie TAOCP Vol. II, Exercise 23. [^5]: MIX is een abstracte computer met een CISC-instructieset. [^6]: Op ARM geeft deling door nul 0 terug. [^7]: __udivti3 is de software-builtin die GCC/Clang/Rust gebruiken voor 128-bit deling. [^8]: x86 genereert een #DE (division error) exception. [^9]: In oudere edities waren vectoren 1-geïndexeerd. [^10]: MMIX is een 64-bit RISC-machine. [^11]: Setun was een Sovjet-computer die gebruikmaakte van gebalanceerd ternair. [^12]: De p-adische representatie $\mathbb{Z}p$ kan worden gezien als inverse limieten $\varprojlim \mathbb{Z}/p^n$. [^13]: Dit zou een nieuwe instructieset vereisen voor de oneven grootte $p$. [^14]: Een tegenvoorbeeld voor de $v{n-2}$ check is $u=(t,0,0,0)b, v=(t,0,1)b$ waar $t=\lfloor b/2 \rfloor$. [^15]: De check kan falen bij overflow als er geen derde register is voor $b^2$. [^16]: Bewijs: Als $\hat{q}=b+1$, dan volgt uit $\hat{r} < v{n-2}$ dat $\hat{q} v{n-2} > b \hat{r} + u_{n-2}$. [^17]: A. Svoboda, Stroje na Zpracování Informací 9 (1963), 25–32.