Chapter 7: Modular Arithmetic

7.1 Basics — Why Avoid Division?

Modular remainder $a \bmod m$ is by definition derived from the division $\lfloor a / m \rfloor$:

As the previous chapter showed, multiprecision division is 2-5 times slower than multiplication. In cryptography and NTT, however, we often perform hundreds of millions of modular multiplications with the same modulus $m$. By pre-computing something about $m$, we can get the remainder with just 1-2 multiplications per call.

7.2 Barrett Reduction

Barrett's idea (1986) was to use a pre-computed "pseudo-reciprocal" $\mu = \lfloor 2^k / m \rfloor$ to replace division with two multiplications.

The algorithm

If $m$ has $n$ words, take $k = 2n \cdot 64$:

$\hat{q}$ is an approximation of the true quotient $\lfloor a/m \rfloor$ with error at most $1$. Hence $r < 2m$, correctable by a single subtraction. Writing $\mu = 2^k/m - \delta_1$ with $0 \le \delta_1 < 1$, the premise $a < 2^k$ gives $a\delta_1/2^k < 1$ and therefore $0 \le a/m - \hat{q} < 2$; together with $\hat{q} \le q$ (which follows from $\mu \le 2^k/m$) this forces $q - \hat{q} \in \{0,1\}$. The weaker bound of 2 quoted in most references belongs to the classical formulation (HAC 14.42), which uses only the top words of $a\mu$ instead of the full product used here.

Why is it fast?

With $m$ occupying $n$ words we have $2^k = \beta^{2n}$ (where $\beta = 2^{64}$), so $\mu = \lfloor \beta^{2n} / m \rfloor \ge \beta^{n}$: in general $\mu$ takes $n+1$ words, not $n$. The product $a \cdot \mu$ is therefore a $2n$-word by $(n{+}1)$-word multiplication.

What the quotient estimate needs, however, is $\lfloor a\mu / \beta^{2n} \rfloor$ — only the part of the product above $\beta^{2n}$. A truncated (high-product) multiplication can supply that without forming the full product, though how much this actually saves depends on the truncated-multiplication algorithm and its implementation. Either way it is far lighter than $O(n^2)$ division.

7.3 Montgomery Multiplication

Montgomery's idea (1985) was to keep integers in a transformed Montgomery form $\bar{a} = a R \bmod m$ (with $R = 2^{64n}$) and perform multiplication + reduction in this form "without any division."

Montgomery form and REDC

The REDC algorithm

This is the magic of Montgomery. Given $T$ with $0 \le T < Rm$ (at most $2n$ words):

1. Compute $u = T \cdot (-m^{-1}) \bmod R$ (only the bottom $n$ words are needed).
2. Compute $T' = T + u \cdot m$. Then $T' \equiv T \pmod m$ and $T' \equiv 0 \pmod R$.
3. $T'' = T' / R$ is just a right shift (drop the bottom $n$ words).
4. If needed, subtract $m$ once from $T''$ to canonicalize.

The computation is two multiplications ($u$ and $u \cdot m$), an add, a shift, and at most one subtraction. Division hides in the pre-computed $m^{-1} \bmod R$. With $m$ fixed, $m^{-1}$ is computed only once.

Cost of input/output conversion

With $R^2 \bmod m$ precomputed, converting an ordinary integer $a$ to Montgomery form costs a single Montgomery multiplication: $\text{REDC}(a \cdot (R^2 \bmod m)) = aR \bmod m$. The inverse conversion is REDC($\bar{a} \cdot 1$). For $k$ multiplications in Montgomery space, the overhead is 2 conversions plus $k$ Montgomery multiplications. Around $k \ge 3$ this starts to beat ordinary reduction, though the actual break-even point depends on the implementation and the size of the modulus.

7.4 The CIOS Algorithm and sangi Implementation

When implementing REDC, one can choose between SOS (Separated Operand Scanning — compute $T = a \cdot b$ fully, then reduce) and CIOS (Coarsely Integrated Operand Scanning — interleave multiplication and reduction). CIOS saves intermediate storage and is among the fastest at typical crypto sizes (n=4-64) in Koç et al.'s comparison; the actual ranking depends on the implementation and the CPU.

Structure of CIOS

The inner loop over each word $b_i$ of $b$ does:

1. Multiply step: $t \leftarrow t + a \cdot b_i$ (fixes one more word of $T$).
2. Reduce step: compute $q \leftarrow t_0 \cdot m_{\text{inv}} \bmod \beta$, then $t \leftarrow (t + q \cdot m) / \beta$ (zeroes out the bottom word).

After $n$ iterations, $t$ is the REDC result.

The generic REDC in sangi

Below is sangi's generic reduction routine. As the name says it is REDC itself: it does not contain the multiplication stage of the CIOS loop described above. The product $T = a \cdot b$ is formed by the caller and this function only reduces it — structurally the SOS arrangement. To make it CIOS, the accumulation of $a b_i$ and of $q_i m$ would have to alternate inside the same outer loop.

// IntModular.cpp: mont_redc (generic REDC; the caller supplies the product T)
// Requires: T has room for 2n+1 words and satisfies 0 <= T < Rm.
// Note: mpn::addmul_1 returns the single-word carry out of the n-word addmul.
//       When tn > i+n, that carry must be propagated up through the higher
//       words (the inner for loop).
inline void mont_redc(uint64_t* r, uint64_t* t, size_t tn,
                      const uint64_t* m, size_t n, uint64_t m_inv) {
    for (size_t i = 0; i < n; ++i) {
        uint64_t q = t[i] * m_inv;                     // q = T[i] * (-m^{-1}) mod beta
        uint64_t carry = mpn::addmul_1(t + i, m, n, q); // T += q*m*beta^i (local result)
        // Propagate the single carry out of addmul_1 upward (only while non-zero)
        for (size_t j = i + n; carry && j < tn; ++j) {
            uint64_t sum = t[j] + carry;
            carry = (sum < t[j]) ? 1 : 0;
            t[j] = sum;
        }
    }
    // The result is t[n..2n-1] TOGETHER WITH the carry word t[2n] above it.
    // When t[2n] != 0 the value is >= beta^n > m, so the subtraction is mandatory.
    const uint64_t top = (tn > 2 * n) ? t[2 * n] : 0;
    std::memcpy(r, t + n, n * sizeof(uint64_t));
    if (top != 0 || mpn::cmp(r, n, m, n) >= 0) {
        mpn::sub(r, r, n, m, n);                        // normalize (r < m)
    }
}

The generic path is independent of $n$. For specific sizes (n=4, 8, 16, 32) we dispatch to hand-written assembly kernels.

7.5 Montgomery n=8 Specialization — sangi Assembly

512 bits (8 words) is a representative benchmark size for multiprecision modular arithmetic. (Real cryptographic moduli differ: the RSA-4096 modulus $n$ is 4096 bits = 64 words, and even the CRT half-moduli $p, q$ are 2048 bits = 32 words; ECDSA P-521 needs $\lceil 521/64 \rceil = 9$ words.) An assembly kernel specialized to this one size runs 30-50% faster than the generic C++ path in our measurements; the margin depends on the CPU and on the baseline it is compared against. The sangi mpn_mont_mul_8 is organized into three phases.

Three-phase structure

Accumulator register residency

The 8-word accumulator $r_8, r_9, \ldots, r_{15}$ sits in GPRs (r8-r15), significantly reducing memory traffic:

; Phase 1 skeleton from mpn_x64_mont.asm
xor r8d, r8d          ; Zero 8 accumulators
xor r9d, r9d
...                   ; through r15
mov rdx, [rcx]        ; b[0] into rdx (implicit MULX operand)
ADDMUL_FIRST_8        ; Seed: r8 = a[0]*b[0]_lo, carry into r9
mov QWORD PTR [rdi], r8  ; r8 is final, store to product buffer
ADDMUL_REST_8         ; Accumulate remaining products into r9-r15
; Repeat for b[1]..b[7]

ADDMUL macros

; ADDMUL_FIRST_8: seed multiplication for b[j] (a[0]*rdx into rbp:rbx)
ADDMUL_FIRST_8 MACRO
    xor eax, eax
    mulx rbx, rbp, [rsi]    ; rbp:rbx = a[0] * rdx (b[j])
    adox r8, rbp            ; r8 += rbp (preserves OF)
ENDM

; ADDMUL_REST_8: accumulate remaining a[1]..a[7] into r9..r15 via two carry chains
ADDMUL_REST_8 MACRO
    mulx rbp, rax, [rsi+8]   ; a[1]*rdx low half
    adcx r8, rbx             ; prior rbx into r8 via CF chain
    adox r9, rax             ; low half into r9 via OF chain
    mulx rbx, rax, [rsi+16]  ; a[2]*rdx
    adcx r9, rbp
    adox r10, rax
    ; ... same pattern unrolled for a[3]..a[7] into r11..r15
    ; final: adcx r15, rbx; adox r15, 0 to merge both carries
ENDM

; REDC_ITER: one reduction round (iter_idx = 0..7)
REDC_ITER MACRO iter_idx
    mov rax, r8
    imul rax, QWORD PTR [rsp+168]   ; rax = r8 * m_inv
    mov rdx, rax                     ; q = rax (implicit MULX operand)
    ADDMUL_FIRST_8                   ; T += q * m (first word)
    ADDMUL_REST_8                    ; T += q * m (remaining 7 words)
    add r15, QWORD PTR [rdi + (8 + iter_idx)*8]  ; merge upper product buffer
    jnc ri_nc
ENDM

MULX / ADCX / ADOX dual carry chains

BMI2's MULX does not affect CF/OF. ADX's ADCX / ADOX treat CF and OF as independent carry flags. This lets two carry chains run simultaneously, halving dependency chains. Standard technique for modern Intel/AMD Montgomery code.

Dispatch

// IntModular.cpp: n==8 branch
if (n == 8 && mpn::detail::has_bmi2_adx()) {
    mpn_mont_mul_8(r, a, b, m, m_inv);
    return;
}

7.6 Chinese Remainder Theorem

Given pairwise coprime moduli $m_1, m_2, \ldots, m_k$, the Chinese Remainder Theorem gives a one-to-one correspondence:

This correspondence is a ring isomorphism $\mathbb{Z}/M\mathbb{Z} \cong \prod_i \mathbb{Z}/m_i\mathbb{Z}$: addition, subtraction and multiplication can all be carried out independently modulo each $m_i$, so large-modulus computations can be split into $k$ smaller-modulus ones. sangi uses CRT both for NTT prime combining and for stand-alone integer CRT used in RSA key generation and polynomial interpolation.

Garner's algorithm

An efficient way to reconstruct $a$ from residues $r_i = a \bmod m_i$:

Finally $a = v_1 + v_2 m_1 + v_3 m_1 m_2 + \cdots$. The key advantage of Garner is that it does not build the full product $M = \prod m_i$ explicitly; each step accumulates $m_1 m_2 \cdots m_{i-1}$ incrementally, so only modular inverses with respect to the small $m_i$ are needed.

// Garner's algorithm (incremental version, pseudocode)
// std::integral is a C++20 concept; on C++17 or earlier, replace with template <typename T>.
// Note: the mixed-radix coefficients v[i] are below m_i and fit in T, but the final
//       reconstruction a = v[0] + v[1]*m_0 + ... is as large as M = prod m_i and will
//       overflow a built-in integer type. Use a multiprecision type for the return
//       value and for the reconstruction loop.
template <std::integral T>
T garner_crt(const std::vector<T>& remainders,
             const std::vector<T>& moduli) {
    std::vector<T> v(moduli.size());
    v[0] = remainders[0];
    for (size_t i = 1; i < moduli.size(); ++i) {
        T t = remainders[i];
        for (size_t j = 0; j < i; ++j) {
            t = ((t - v[j]) * mod_inverse(moduli[j], moduli[i])) % moduli[i];
            if (t < 0) t += moduli[i];
        }
        v[i] = t;
    }
    // Finally assemble a = v[0] + v[1]*m_0 + v[2]*m_0*m_1 + ... in multiprecision
    T a = v.back();
    for (size_t i = moduli.size() - 1; i-- > 0;) a = a * moduli[i] + v[i];
    return a;
}

A simpler Lagrange-style CRT that precomputes $M$ ($a = \sum r_i \cdot M_i \cdot (M_i^{-1} \bmod m_i) \bmod M$ with $M_i = M / m_i$) is also viable and performs comparably for small $k$. Choose based on whether building $M$ up front is cheaper than repeated small modular inverses.

7.7 Fast Modular Exponentiation — Sliding Window

The simplest way to compute $m^e \bmod n$ is to process the binary expansion of $e$ bit-by-bit (binary method): $O(\log e)$ multiplications.

$k$-bit window method

Processing $k$ consecutive bits at a time shortens the addition chain further. Pre-compute $m^1, m^3, m^5, \ldots, m^{2^k - 1}$ (odd powers), then at each window $w$:

The sangi window width selection:

// IntModular.cpp: choose_window_width (sangi, as of 2026-05)
inline int choose_window_width(size_t expBits) {
    if (expBits <= 24)   return 1;
    if (expBits <= 64)   return 3;
    if (expBits <= 256)  return 4;
    if (expBits <= 1024) return 5;
    if (expBits <= 4096) return 6;
    return 7;
}

Sliding window

A fixed window still multiplies when the window starts with 0. A sliding window skips runs of zeros (only squaring in that interval) and aligns each window to start with a 1. On average this saves on the order of 20-30% of the multiplications; the actual figure depends on the bit pattern of the exponent and on the window width.

// IntModular.cpp: powerMod (sliding window, sketch)
int oddTableSize = 1 << (w - 1);   // pre-compute only odd powers
uint64_t g[oddTableSize];
g[0] = baseR;                       // m^1 (Montgomery form)
if (oddTableSize > 1) {
    uint64_t base2R = mont_mul(baseR, baseR);    // m^2
    for (int i = 1; i < oddTableSize; ++i)
        g[i] = mont_mul(g[i-1], base2R);         // m^3, m^5, m^7, ...
}
// Main loop: scan exponent bits left-to-right, detect windows, multiply

Since everything happens in Montgomery form, each mont_mul is division-free. The result is converted back at the end.

7.8 Constant-Time-Oriented Implementation — powerModSec

In cryptographic use, if compute time depends on secret keys, a timing attack can recover them. Since Kocher's 1996 attack, real-world timing attacks have been reported periodically (e.g., CVE-2018-0737).

What affects timing?

• Branches: if branch prediction misses, cycles vary.
• Cache access patterns: if a lookup table index is secret, L1 hit/miss can vary by 10×.
• Data-dependent instructions: on some CPUs, division and multiplication cycles depend on operand values.

Constant-time modular exponentiation (sangi powerModSec)

The sangi powerModSec changes three things from the ordinary powerMod:

1. Fixed window (no sliding): always performs the same number of multiplications. If the window is 0, still multiplies by $m^0 = 1$.
2. Masked table lookup: reads every table entry and selects via bitmask, so the cache access pattern never depends on secrets.
3. Branch-free: uses arithmetic masks rather than ifs.

// IntModular.cpp: powerModSec constant-time table select
int tableSize = 1 << w;  // all 0..2^w-1 (not only odd)
std::vector<uint64_t> g_buf(tableSize * n);
// No sliding: all entries always ready

auto ct_select = [&](uint64_t* dst, int idx) {
    std::memset(dst, 0, n * sizeof(uint64_t));
    for (int i = 0; i < tableSize; ++i) {
        // i==idx yields 1, otherwise 0. Cast to int64_t and negate to get
        // 0xFFFF...FFFF (all ones) or 0x0000...0000 via two's-complement semantics.
        uint64_t mask = -(static_cast<int64_t>(i == idx));
        for (size_t j = 0; j < n; ++j)
            dst[j] |= g_buf[i * n + j] & mask;              // branch-free masked OR
    }
};

Further mitigations

A fully constant-time implementation also considers:

• Blinding: for RSA decryption $c^d \bmod n$, draw a random $r$, form the blinded input $c^{\prime} = c\,r^e \bmod n$, compute $m^{\prime} = (c^{\prime})^d \bmod n$, and recover $m = m^{\prime} r^{-1} \bmod n$. The secret exponent then acts on a fresh input on every invocation.
• A fixed memory access pattern: lay the table out regularly and walk the same group of cache lines in the same order regardless of the secret. At 512 bits one entry is exactly 64 bytes = one cache line, but for larger moduli an entry no longer fits in a single line; the point is not to fit in one line but to keep the access sequence independent of the secret.
• Power analysis (DPA) mitigations: for physical attacks, masked (secret-shared) computation may be required.

Since sangi is not a cryptographic library, we do not implement blinding or DPA countermeasures. For real cryptographic use, rely on established libraries such as OpenSSL or libsodium.

7.9 Summary

• Barrett: pre-compute $\mu = \lfloor 2^k / m \rfloor$ to replace division with two multiplications. Best for single-shot remainders.
• Montgomery: use Montgomery form $\bar{a} = a R \bmod m$ to chain multiplications without division. REDC is the key. Dominates cryptographic/NTT hot paths.
• CIOS: interleaved multiply-and-reduce pattern. sangi's generic path instead forms the product first and then reduces it (REDC alone), dispatching to specialized assembly at n=4, 8, 16, 32.
• mpn_mont_mul_8: 512-bit specialization. 8 accumulators resident in r8-r15. Dual carry chains with MULX/ADCX/ADOX.
• CRT: parallel computation over coprime moduli, combined by Garner. RSA-CRT is standard optimization.
• Sliding-window modular exponentiation: pre-compute odd powers, skip zero runs for roughly 20-30% fewer multiplications.
• Constant-time orientation: powerModSec uses fixed windows + masked selects to make the skeleton constant-time. Data-dependent work remains in the underlying Montgomery multiplication, so a whole-routine claim requires machine-code-level verification.

The next chapter, Chapter 8: Fast Computation of Number-Theoretic Functions, covers multiplicative functions — Euler's $\phi$, the Möbius $\mu$, the divisor sums $\sigma_k$ — and the fast evaluation of the prime-counting function $\pi(x)$. (An English translation of that chapter is not yet available; the Japanese edition is linked from the series index.)

References

• Barrett, P. "Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor", CRYPTO '86, LNCS 263, pp. 311-323, 1987.
• Montgomery, P.L. "Modular Multiplication without Trial Division", Mathematics of Computation, 44(170), pp. 519-521, 1985.
• Menezes, A.J., van Oorschot, P.C., Vanstone, S.A. Handbook of Applied Cryptography, CRC Press, 1996. Chapter 14.
• Brent, R.P., Zimmermann, P. Modern Computer Arithmetic, Cambridge University Press, 2010. Chapter 2.
• Koç, Ç.K., Acar, T., Kaliski, B.S. "Analyzing and Comparing Montgomery Multiplication Algorithms", IEEE Micro, 16(3), pp. 26-33, 1996. Comparison of CIOS/SOS/FIOS.
• Kocher, P.C. "Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems", CRYPTO '96, LNCS 1109, pp. 104-113.
• Bernstein, D.J. "Curve25519: New Diffie-Hellman Speed Records", PKC '06, LNCS 3958, pp. 207-228. Constant-time implementation in practice.

Related Chapters

• Chapter 6: Fast Division Algorithms — the naive implementation of modular reduction
• Chapter 4: Integer GCD Algorithms — extended Euclidean and the computation of $m^{-1}$
• Chapter 10: NTT Implementation Comparison — Montgomery multiplication in NTT
• Chapter 11: Primality Testing — Miller-Rabin exponentiation
• Chapter 12: Integer Factorization — modular arithmetic in Pollard rho and ECM