Een snellere manier om signed 64-bit floor division door constanten te berekenen

Wat is signed integer floor division?

Signed integer floor division is een delingsoperatie ($a/b$) waarbij het resultaat wordt afgerond naar $-\infty$, terwijl standaard ("truncated") signed division afrondt naar $0$.

Zoals vermeld in Hacker's Delight: "[Ik vermoed dat] floor division vaker het resultaat geeft dat je wilt". Dit artikel richt zich op full-range floor division waarbij de divisor (de deler) bekend is tijdens het compileren.

Ondanks het nut van deze functie is de wereld van integer floor division in een deplorabele staat:

  • Mainstream compilers en veel talen (bijv. C/C++/Go/Swift) bieden het niet native aan.
  • Talen en softwarebibliotheken die het wel implementeren, gebruiken vrijwel nooit de snelste algoritmen.
  • Ten minste één zeer prominente softwarebibliotheek had ruim vier jaar lang een bugge implementatie die leidde tot overflows.
  • Typische online zoekopdrachten en LLM-prompts brengen meestal alleen de trage algoritmen aan het licht.
  • Tot slot heeft Hacker's Delight (in feite de bijbel voor low-level computermathematica) bijna niets over dit onderwerp te zeggen.

Dit blogbericht heeft tot doel niet alleen licht te werpen op de snelle algoritmen, maar introduceert ook een nieuwe snelle aanpak die, bij ideale compilatie, de volgende resultaten behaalt:

  • Op ARM: 10–20% snellere throughput dan bestaande algoritmen voor de meeste 64-bit divisors (soms gelijk), met gelijke of betere latentie.
  • Op x86:
  • ~5-10% snellere throughput voor ongeveer 80% van de positieve 64-bit divisors (60% van de negatieve divisors).
  • Trager in de overige x86-gevallen, en meestal een slechtere (maar concurrerende) latentie.

Het werkt voor alle bit-breedtes en kan nuttig zijn voor sommige 32-bit of SIMD-applicaties, maar het grootste voordeel ligt bij 64-bit. Het kan ook worden gebruikt voor snellere positive-modulus en snellere "Euclidean" division (die verschilt van floor division bij negatieve divisors).

Benchmark voorbeeld: Deling door 86400

Floor division kan veel sneller zijn dan je waarschijnlijk realiseert. Hieronder staat een snapshot van de throughput-benchmark voor deling door 86400 (het aantal seconden in een dag). In de grafieken geldt: lager is sneller.

ARM64 Chips (MacBook Pro M4 Pro, Snapdragon X1-26-100)

  • Common (Rust, Zig, etc.): ~2×
  • Robison (2005): ~1.1×
  • Granlund-Montgomery (1994): 1×
  • New Algo (2026): ~0.8×

x86 Chips (AMD Ryzen 9 9950X3D, Intel Core i7 MacBook Pro 2016)

  • Common (Rust, Zig, etc.): ~2×
  • Granlund-Montgomery (1994): 1.05×
  • Robison (2005): ~1×
  • New Algo (2026): ~0.95×

\Rust biedt div_euclid, wat overeenkomt met floor-div voor positieve divisors.*

Voorbeeldimplementatie: Snelle Floor Division door 86400

Hieronder staat een implementatie voor 64-bit under MIT-licentie:

// Fast floor division by 86400.
// (rounds towards negative infinity).
// Accurate over the full signed 64-bit range.
// Sign-selected multiplier variant.
inline int64_t fast_i64_floordiv_86400(int64_t x) {
    const int64_t M = ((__int128(1) << 79) / 86400 + 1);
    const int64_t m = M + (x >> 63);
    return (__int128) x * m >> 79;
}

---

Beginnen met Truncated (Standaard) Division

Stel dat je een Unix-timestamp in seconden wilt converteren naar rata-die (een dagteller). Als je hiervoor standaard ("truncated") division gebruikt, deel je door 86400.

// C++ Standard Division (niet floor-div)
// Given: x — signed 64-bit integer
const int64_t rd = x / 86400;
return rd;

Deze functie geeft incorrecte resultaten voor de meeste negatieve waarden (timestamps van vóór 1970). Voor het tracken van herhalende patronen zoals timestamps, hoekrotaties of tegelposities is floor division essentieel. Hierbij is de bijbehorende % operatie een echte positieve modulus, in tegenstelling tot de "fake" modulus die % biedt in C++ of JavaScript.

Wat gebeurt er onder de motorkap?

De gecompileerde assembly voor standaard deling is in feite een mul-shift met een "bump-by-1" voor negatieve getallen. De pseudocode ziet er als volgt uit:

const i64 M = (i128(1) << 77) / 86400 + 1; // Precomputed
const i128 wide = i128(x) * M;
const i64 tbit = wide >> 77;      // high64(wide) >> 13
const i64 sign = u64(x) >> 63;    // x < 0 ? 1 : 0
return tbit + sign;

De stap u64(x) >> 63 verplaatst het hoogste bit (het teken-bit) naar de onderste positie, waardoor er een 1 overblijft voor negatieve getallen en een 0 voor positieve.

Uncooperative Divisors

Niet elke divisor is zo eenvoudig. Bij deling door bijvoorbeeld 100 is de "magic multiplier" $M$ zo groot dat alle 64 bits worden benut. Dit maakt $M$ een negatief getal, waardoor de waarde van $M$ $2^{64}$ lager uitvalt dan nodig. Om dit te corrigeren, moet $x$ worden toegevoegd aan het bovenste register van het resultaat. Divisors die deze correcties vereisen, worden "uncooperative divisors" genoemd.

---

Het oplossen van de bug

Wanneer je merkt dat x / 86400 fouten geeft bij negatieve waarden, is de meest voor de hand liggende fix:

const int64_t div = x / 86400;
const int64_t rem = x % 86400;
return div - (rem < 0);

Echter, het berekenen van de modulus is bijna net zo duur als de deling zelf. Een veelvoorkomende "optimale" oplossing (zoals gebruikt in Java Math.floorDiv, Python numpy, Zig en Ruby) is:

// Common Solution for signed divisors
const int64_t div = x / y;
const int64_t rem = x % y;
return div - (rem != 0 && (x ^ y) < 0);

Wat zegt Hacker's Delight hierover?

In Hacker's Delight, 2e Editie, staat in de inhoudsopgave een sectie: "10–12 Applicability to Modulus and Floor Division". De volledige inhoud van deze sectie is echter slechts één alinea:

“Het zou kunnen lijken dat het omzetten van modulus of floor division door een constante in een vermenigvuldiging eenvoudiger zou zijn, omdat de stap 'voeg 1 toe als het dividend negatief is' zou kunnen worden weggelaten. Dit is niet het geval. De bovenstaande methoden zijn niet op een voor de hand liggende manier toepasbaar op modulus en floor division. Misschien kan er iets worden bedacht; dit zou kunnen inhouden dat de multiplier $m$ enigszins wordt gewijzigd, afhankelijk van het teken van het dividend.” — Henry S. Warren, Jr.

Hoewel Warren geen oplossing biedt, voorspelt hij exact de aanpak die later in dit artikel wordt geïntroduceerd.

---

Snellere Technieken

De methode van ufoguy

Op Stack Overflow stelde gebruiker ufoguy in 2017 een snellere variant voor voor positieve divisors ($b > 0$):

int floor_div(int a, int b) {
    return (a + (a < 0)) / b - (a < 0);
}

Benchmarks tonen aan dat dit de throughput met ongeveer 50% verbetert.

Granlund-Montgomery (1994)

Voor een algemene oplossing (inclusief negatieve divisors) kan men terugvallen op de techniek uit de paper Division by Invariant Integers using Multiplication (1994). Een C++ template-implementatie ziet er als volgt uit:

template <auto D>
constexpr decltype(D) floor_div(decltype(D) x) {
    static_assert(std::is_integral_v<decltype(D)> &&
                  std::is_signed_v<decltype(D)> && D != 0,
                  "divisor must be a signed non-zero integer");
    if constexpr (D > 0)
        return (x + (x < 0)) / D - (x < 0);
    else
        return (x - (x > 0)) / D - (x > 0);
}

De Positive Divisor XOR Trick

Granlund en Montgomery bieden een nog snellere route voor positieve divisors via een XOR-truc:

  1. s = x >> 63 (is -1 voor negatief, 0 voor positief)
  2. a = x ^ s (verplaatst negatieve inputs naar het non-negatieve domein: negatieve waarden worden $-1 - x$)
  3. b = (u128) a * M >> S
  4. c = b ^ s
  5. return c

Dit is effectief een reflectie rond $-0.5$. Het grote voordeel is dat $a$ niet alleen unsigned is, maar ook een beperkt bereik heeft ($0..2^{63}-1$), waardoor een mul-shift altijd werkt, zelfs voor uncooperative divisors.

---

Robison’s Multiply-Add Techniek

In de paper N-Bit Unsigned Division Via N-Bit Multiply-Add (2005) introduceert Arch D. Robison een methode die vaak reduceert tot een simpele mul-add-shift zonder branches of ternaire operaties.

Voorbeeld voor 32-bit deling door 100:

// Floor Div 100 (Robison)
i32 A = 1374389535;
i32 B = 0x28F50000; // 687144960
return (i64(x) * A + B) >> 37;

Op ARM compileert dit naar een enkele multiply-accumulate operatie gevolgd door een right-shift. Er zijn echter twee problemen:

  1. Wijdere vermenigvuldiging: Sommige divisors (zoals 101) vereisen een multiplier $A$ die net één bit te groot is voor de native breedte, waardoor er overgestapt moet worden naar een 64-bit type (voor 32-bit input).
  2. Complexiteit bij 64-bit: Bij 64-bit moet men werken met 128-bit emulatie. Dit introduceert een "add back" correctie voor negatieve $A$, wat op ARM64 resulteert in twee vermenigvuldigingsoperaties (mul en smulh), wat de winst beperkt.

---

De Nieuwe Aanpak: Sign-Selected Multiplier

Als we kijken naar een simpele mul-shift voor deling door 3, zien we dat een "round-up" multiplier perfect werkt voor positieve getallen, maar een discontinuïteit vertoont bij negatieve getallen. Een "round-down" multiplier doet precies het tegenovergestelde.

De oplossing is om een verschillende multiplier te gebruiken voor positieve en negatieve getallen:

// Floor Division by 3 (Sign Selected Multiplier)
i64 M = MAX_U64 / 3 + 1;
i64 m = M + (x >> 63);  // = x < 0 ? M - 1 : M
i64 y = (i128)x * m >> 64;

Classificatie van Divisors

Niet alle divisors zijn even flexibel. Ze kunnen worden ingedeeld in:

  • Round-Up + Negative Bump: Gebruikt een round-up multiplier en past $x$ aan als $x < -1$.
  • Round-Down + Positive Bump: Gebruikt een round-down multiplier en past $x$ aan in het bereik $0 \le x < \text{MAX\_INT64}$.
  • Sign-Selected Multiplier: Kan beide methoden gebruiken en kiest de multiplier op basis van het teken.

Negatieve Divisors

Dit principe werkt ook voor negatieve divisors. Men moet alleen voorzichtig zijn met de tekens en vergelijkingen. De implementatie past de grenzen aan op basis van het teken van de divisor.

---

Floor-Modulus en Euclidean Division

Floor modulus vloeit direct voort uit floor division: posmod = x - fastfloor_div<Divisor>(x) * Divisor

Euclidean division is identiek aan floor division voor positieve divisors. Voor divisors kleiner dan -1 kan de volgende identiteit worden gebruikt: eucdiv = -fastfloor_div<-Divisor>(x)

---

Implementatie en Keuze van Algoritme

Het kiezen van de juiste variant gebeurt tijdens het compileren door de worst-case scenario's te testen voor zowel de round-up als de round-down versies.

Welk algoritme moet ik gebruiken?

  • 8, 16 en 32-bit: Gebruik Robison. Het is simpel en zeer efficiënt op deze breedtes.
  • 64-bit ARM: Gebruik het Nieuwe Algoritme (Sign-Selected). Het is superieur of gelijk op alle metrics.
  • 64-bit x86: De keuze is complexer.
  • Voor maximale throughput: Het Nieuwe Algoritme wint in de meeste gevallen.
  • Voor minimale latentie: Robison wint voor bijna alle divisors.
  • Voor extreme optimalisatie: Gebruik een hybride aanpak waarbij je per divisor kiest tussen New en Robison.

---

Waarom je niet je eigen Floor Division algoritme moet schrijven

Het is zeer eenvoudig om fouten te maken. Een voorbeeld is de bug in Google's TensorFlow, die ruim vier jaar (2017-2021) in de codebase zat.

De foutieve TensorFlow implementatie (pseudocode):

if ((x < 0) != (y < 0)) {
    const i64 abs_x = abs(x);
    const i64 abs_y = abs(y);
    return -(abs_x + abs_y - 1) / abs_y; // OVERFLOW BUG HIER
} else {
    return x / y;
}

De term absx + absy - 1 veroorzaakt een overflow voor extreme waarden in de signed 64-bit range. Dit onderstreept het belang van het gebruik van bewezen algoritmen of native taalondersteuning.

---

Bijlage A: Bewijs (Proof)

Positieve divisor, non-negatieve x

Voor $x \ge 0$ is $x$ een unsigned $(W-1)$-bit getal. De resultaten van Granlund-Montgomery en ridiculousfish zijn direct toepasbaar. Zowel de round-up multiplier $M$ als de round-down multiplier $M'$ (met increment) zorgen ervoor dat het resultaat strikt tussen $q$ en $q+1$ ligt wanneer $n > 0$.

Negatieve x

Laat $u := \sim x = -1 - x$. De bitwise NOT mapt negatieve getallen één-op-één naar $0 \le u \le 2^N-1$. Er geldt: $\lfloor x/D \rfloor = -1 - \lfloor u/D \rfloor$ (voor $D > 0$).

  1. Round-up, $x \le -2$: Hier geldt $1 \le u \le 2^N-1$. De identiteit $\lfloor -y \rfloor = -1 - \lfloor y \rfloor$ voor non-integers $y$ zorgt ervoor dat de bump correct werkt.
  2. Round-up, $x = -1$: Deze waarde wordt overgeslagen door de bump omdat $u=0$. Omdat $D \ge 3$, valt de ratio strikt tussen $(-1, 0)$, wat resulteert in $\lfloor -1/D \rfloor = -1$.
  3. Round-down, $x \le -1$: Hier is $0 \le u \le 2^N-1$. De increment variant werkt voor elk $u \ge 0$.

Sign-selected

Dit combineert de on-bumped helften van zowel round-up als round-down. Het werkt bij elke shift waar beide individuele roundings werken.

Negatieve divisors

$\lfloor x/-|D| \rfloor = \lfloor -x/|D| \rfloor$. Dit reduceert het probleem tot de positieve-divisor case toegepast op $-x$, met kleine aanpassingen in het bereik vanwege de asymmetrie van signed integers (INT\_MIN heeft geen positieve tegenhanger).