Faster modular exponentiation of integers

[fredrik-johansson]

Developed using Claude Fable 5.

Adds flint_mpn_powm and flint_mpn_powm_preinvn and uses these when applicable in fmpz_powm, fmpz_powm_ui, fmpz_mod_pow_fmpz, fmpz_mod_pow_ui, and adds mpn_mod_pow_*.

For large moduli and exponents this uses Montgomery reduction with several FFT tricks; a detailed writeup by Claude is attached below.

The speedup is ~1.3x asymptotically and up to 1.8x in a narrow band (moduli around 65K bits).

Speedup for fmpz_powm for a modulus with pb bits and exponent with eb bits:

    pb \ eb     1     2     4     8    16    32    64   128   256   512  1024
       32   2.500 1.264 1.014 1.028 1.271 1.465 1.000 1.000 1.000 1.000 1.026
       64   0.925 0.986 1.000 1.033 1.032 1.023 1.000 1.017 1.000 1.000 1.000
      128   0.928 0.987 1.016 1.000 1.000 1.000 1.013 1.008 0.996 1.021 1.011
      256   0.923 0.966 1.000 1.030 1.019 1.000 1.000 0.974 1.000 1.000 1.000
      512   0.923 0.981 1.023 1.013 1.008 1.004 1.021 1.023 1.006 1.029 1.000
     1024   0.929 0.969 1.000 1.000 1.000 1.000 1.006 1.000 1.000 1.000 1.000
     2048   0.925 1.000 0.979 1.010 1.007 1.033 1.000 1.000 1.005 1.000 1.012
     4096   0.890 1.031 1.000 1.000 1.000 1.019 0.995 1.000 1.000 1.000 1.000
     8192   0.914 1.000 1.027 1.000 1.000 1.000 1.000 1.008 1.000 1.022 0.989
    16384   0.918 0.996 1.073 1.428 1.113 1.042 1.039 1.057 1.059 1.053 1.050
    32768   0.961 1.000 1.405 1.654 1.333 1.441 1.500 1.514 1.515 1.544 1.540
    65536   0.978 0.995 1.613 1.427 1.522 1.816 1.849 1.851 1.875 1.904 1.883
   131072   0.996 1.000 1.000 1.000 1.111 1.172 1.175 1.290 1.271 1.283 1.273
   262144   0.931 1.000 1.000 1.018 1.052 1.181 1.229 1.266 1.275 1.288 1.300
   524288   1.010 0.990 1.000 1.000 1.011 1.205 1.300 1.311 1.319 1.322 1.332
  1048576   1.000 1.000 1.008 0.992 1.175 1.193 1.271 1.302 1.310 1.337 1.345
  2097152   1.064 1.012 1.000 1.015 1.113 1.168 1.238 1.275 1.297 1.307 1.310
  4194304   0.989 1.000 1.000 1.007 0.959 1.160 1.216 1.243 1.257 1.268 1.264
  8388608   0.879 1.024 1.011 1.007 1.140 1.172 1.232 1.257 1.275 1.278 1.286
 16777216   0.907 0.978 1.000 1.009 1.012 1.206 1.261 1.288 1.296 1.319 1.315

Speedup fmpz_powm vs mpz_powm:

    pb \ eb     1     2     4     8    16    32    64   128   256   512  1024
       32   1.242 3.937 2.939 2.583 2.169 1.951 0.943 0.951 0.973 0.990 1.026
       64   0.587 0.684 0.803 0.744 0.809 0.867 0.970 0.983 0.991 0.990 1.026
      128   0.565 1.012 1.422 0.894 0.924 0.953 0.988 0.993 0.992 1.000 1.011
      256   0.551 1.081 0.995 0.943 0.982 0.981 0.995 1.000 1.014 1.014 1.000
      512   0.571 1.185 0.955 0.974 0.983 0.987 1.000 1.011 0.994 1.000 1.000
     1024   0.614 2.273 1.154 0.986 1.022 1.000 1.006 1.000 1.016 1.000 1.000
     2048   0.717 1.634 1.125 1.010 1.013 1.000 1.017 1.000 1.009 1.000 1.012
     4096   0.863 2.833 1.333 1.000 1.000 1.000 1.005 1.026 1.000 1.000 1.000
     8192   1.028 2.053 1.421 1.000 1.000 1.000 1.000 0.992 0.996 1.000 1.011
    16384   1.222 3.298 1.477 1.423 1.093 1.021 1.033 1.057 1.059 1.037 1.038
    32768   1.615 2.134 1.632 1.673 1.308 1.442 1.459 1.507 1.500 1.521 1.510
    65536   1.804 3.364 1.855 1.414 1.496 1.769 1.851 1.866 1.891 1.887 1.894
   131072   1.920 2.682 2.189 2.032 1.889 2.136 2.209 2.226 2.333 2.321 2.321
   262144   2.000 3.155 2.405 2.500 2.068 2.460 2.543 2.635 2.655 2.678 2.668
   524288   2.019 5.445 2.492 2.540 2.164 2.605 2.714 2.761 2.779 2.777 2.793
  1048576   2.115 3.451 2.595 2.571 2.255 2.628 2.778 2.848 2.889 2.907 2.908
  2097152   1.960 6.190 2.697 2.242 2.381 2.807 2.949 3.032 3.073 3.092 3.112
  4194304   1.950 6.190 2.595 2.302 2.320 2.729 2.870 2.926 2.961 2.975 2.974
  8388608   2.020 3.941 2.841 2.533 2.610 3.096 3.251 3.317 3.342 3.350 3.362
 16777216   1.975 7.170 2.951 2.468 2.791 3.432 3.599 3.678 3.719 3.723 3.703

The algorithms behind flint_mpn_powm (AI-written)

This document describes the algorithms implemented in
src/mpn_extras/powm.c for computing

r = b^e mod m,        0 <= r < m,

with b, e, m given as mpn limb vectors (limb base B = 2^64,
modulus of mn limbs, exponent of en limbs). It states the precise
formulas used at each size tier and explains the underlying techniques.
A final section lists optimizations that remain on the table.

Throughout, "one mul" means one full mn x mn-limb product at the
relevant size; costs of partial products are quoted as fractions of it.

---

1. Overall structure and size tiers

Exponentiation is performed by a sliding-window square-and-multiply
ladder (section 2) whose inner operation is a modular multiplication
or squaring. Everything else in the file is about making that inner
operation cheap and about amortizing per-call setup. The tiers, keyed
by the modulus size mn and the exponent bit length:

condition	inner multiplication
exponent below 25 bits	Barrett mulmod_preinvn ladder, minimal setup (§4)
mn < 110 (large exponents)	delegated to mpz_powm (GMP's redc_1/redc_2)
110 <= mn < 120	Barrett mulmod_preinvn ladder (§4)
120 <= mn < 480	Barrett ladder with wraparound remainder (§4.3)
mn >= 480	Montgomery ladder with folded REDC (§5), quotient in the transform domain (§5.2)

Inside the large tier, the q*m mod (B^rlen - 1) product at the heart
of the folded REDC is itself computed by one of three engines chosen at
setup time (§5.3): a cyclic FFT plan, the recursive
"chain" (§7), or mpn_mulmod_bnm1.

Even moduli are handled by the 2-adic splitting m = 2^t * m_odd:
the odd part runs through the machinery above, the power modulo 2^t
is computed by a dedicated truncated ladder (powm_2exp, all
arithmetic modulo a power of two is plain low products), and the two
residues are recombined by CRT. The rest of this document assumes m
odd where Montgomery arithmetic is involved.

At the fmpz level, two fast paths run before any of this: an exact
power when e <= (bits(m) - 1)/bits(b) guarantees |b|^e < m, and a
mpz_powm_ui call for short bases with exponents under 25 bits below
the FFT tier, where per-call setup cannot be amortized.

---

2. Sliding-window exponentiation

For an exponent of l bits and window parameter k, the ladder
precomputes the table of odd powers

T[j] = b^(2j+1) mod m,        j = 0, ..., 2^(k-1) - 1,

which costs one squaring plus 2^(k-1) - 1 multiplications. The
exponent is then consumed from the top: at each step the scanner
(next_window) finds the longest window of at most k bits that ends
in a set bit, so every table multiplication is by an odd power, and
runs of zeros between windows cost only squarings. For a random
exponent this gives

l  squarings  +  approximately l/(k+1)  table multiplications.

k grows with l (as in _gr_pow_mpn_sliding) but is capped so that
the table of 2^(k-1) full-size entries stays modest. When the base
fits in one or two limbs (small_base), table multiplications are
replaced by mul_scalar_mod: a full product by the 1–2-limb scalar
followed by one division — asymptotically negligible next to the
squarings — with base 2 specialized further to a shift and a
conditional subtract.

The measured cost split at mn = 1030 reflects this design: the
squarings and the reductions dominate, the table multiplications are
a few percent.

---

3. Squaring

At every tier the squarings go through flint_mpn_sqr. In the FFT
range this detects the equal operands and computes only one forward
transform, one pointwise pass, and one inverse transform, so that

cost(sqr) ~ 0.74 * cost(mul)

(measured at 1030 limbs: 78 vs 105 microseconds). Since half or more
of the ladder's work is squarings, this factor is a large share of the
overall advantage over a generic-multiplication ladder.

---

4. The Barrett basecase: precomputed inverse of the modulus

4.1 The precomputed inverse

flint_mpn_preinvn computes, for the modulus shifted into normalized
position (norm = clz(m[mn-1]), d = m << norm),

dinv  ~  floor( (B^(2n) - 1) / d )  -  B^n,

an n-limb approximation of the scaled reciprocal. The powering
ladder keeps residues shifted left by norm so that every reduction
sees a normalized divisor.

4.2 Barrett reduction (flint_mpn_mulmod_preinvn)

Given a 2n-limb product X < d * B^n, the quotient estimate is

q = X_hi + mulhigh(X_hi, dinv),         X_hi = floor(X / B^n),

which satisfies q <= floor(X/d) <= q + k for a small constant k.
The remainder candidate

r0 = X - q*d      (an n-limb low product: mullow(q, d), subtracted)

then lies in [0, k*d) and at most k conditional subtractions of
d produce the canonical remainder. One reduction therefore costs one
mulhigh plus one mullow, about one full multiplication, on top of
the product being reduced.

4.3 The wraparound variant (mulmod_preinvn_fold, 120–2000 limbs)

The mullow(q, d) above only ever feeds a subtraction whose result is
known to be tiny. It can therefore be replaced by a product in the
ring Z / (B^g - 1) for a fold length g slightly larger than the
bound on r0: with the same quotient estimate,

r0 = ( X - q*d ) mod (B^g - 1)

computed by mpn_mulmod_bnm1 recovers r0 exactly because
0 <= r0 < B^g - 1. Reducing an argument modulo B^g - 1 is just
folding: split it into g-limb chunks and add them with wraparound
carry, since B^g ≡ 1. The wraparound product costs roughly half of
a full product at these sizes, which is where the 120-limb tier gets
its advantage.

---

5. The large tier: Montgomery reduction with a folded quotient ledger

5.1 Montgomery representation and REDC

For odd m, fix R = B^mn. Residues are kept in Montgomery form
x~ = x*R mod m; the product of two Montgomery forms followed by

REDC(X) = X * R^(-1) mod m

yields the Montgomery form of the product, and the final result is
recovered by one more REDC. The inverse used is

minv = m^(-1) mod B^mn,

computed by _flint_mpn_binvert via Hensel lifting
v <- v*(2 - m*v) mod B^(2^i), doubling the precision each step.

The classical REDC formula is

q = (X mod B^mn) * minv  mod B^mn          (one mullow)
t = (X + q*m) / B^mn                        (exact division: the
                                             low mn limbs cancel)

with t < 2m when X < m*B^mn, so one conditional subtraction
canonicalizes. Note the signs: because q*m ≡ X (mod B^mn) — this
implementation uses q = X_lo * minv and adds q*m, making the low
half X_lo + q*m ≡ 2*X_lo... in fact the cancellation is arranged as
follows, which is exactly what the folded version exploits.

5.2 The quotient step, optionally in the transform domain

The mullow computing q is the one dense product of the reduction
that is not shared with anything else. Above
FLINT_MPN_POWM_REDC_QSTEP_FFT_THRESHOLD = 480 limbs it is computed
in the transform domain against a cached transform of minv
(Fminv, built once at setup):

q = export( ifft( fft(X_lo) . Fminv ) )  truncated to mn limbs,

replacing an mn x mn mullow by one forward transform, one pointwise
pass, one inverse transform and one export — measured at 1030 limbs
this takes the whole reduction from 157 to 128 microseconds per call.
The threshold was tuned by same-process A/B on the merged pipeline;
the crossover sits between 420 and 500 limbs.

5.3 The folded ledger: recovering t from a wraparound residue

Writing X = X_hi * B^mn + X_lo, the quantity X_lo + q*m is
divisible by B^mn by construction; call the quotient word

c = (X_lo + q*m) / B^mn,        so that   t' = X_hi + c

is the REDC output before canonicalization (the code short-circuits
X_lo = 0, where q = 0 and c = 0, to t = X_hi directly). The
entire dense work is thus computing q*m — and since only c is
needed, q*m never has to be produced exactly: it suffices to know it
modulo B^rlen - 1 for any rlen > mn, because division by B^mn
in that ring is a rotation. Concretely, with k = rlen - mn:

S =  q*m                mod (B^rlen - 1)      (the fold engine)
U =  S - (B^mn - X_lo)  mod (B^rlen - 1)      (= q*m + X_lo - B^mn
                                                 = (c - 1) * B^mn)
H =  U * B^k            mod (B^rlen - 1)      (rotate left k limbs;
                                                 B^(mn+k) = B^rlen ≡ 1)
  =  c - 1.

Since 0 <= H <= m - 1 < B^rlen - 1, the only ambiguity of the
wraparound representative — H = 0 versus the all-ones vector — is
resolved by inspecting the top limbs, and t' = X_hi + H + 1 follows
by one addition. The reduction has become: one quotient product (§5.2),
one multiplication modulo B^rlen - 1, and linear work.

The product S = q*m mod (B^rlen - 1) is served by one of three
engines chosen at setup:

• method 0, cyclic plan: a cyclic FFT convolution of length
  rlen = nnc chosen by fft_small_plan_init_mpn_cyclic as the
  cheapest admissible wraparound length at or slightly above mn,
  multiplied pointwise against a cached transform of m (Fm).
  A cyclic convolution is multiplication modulo B^N - 1: indexing
  chunks by powers of B^bits, the convolution wraps index N back
  to 0. Only q is transformed per call. Cost ~ 0.55 mul.
• method 1, the chain (§7): a recursive CRT factorization of
  B^rn - 1, rn the smallest power of two >= max(mn, 128), with
  cached negacyclic transforms of the modulus residues at each level.
• method 2, mpn_mulmod_bnm1: GMP-style wraparound recursion, used
  in the mid band where FFT plans do not yet pay.

Methods 0 and 1 are compared at setup by the calibrated cost model
27 * rn_chain <= 33 * nnc (chain wins when its power-of-two padding
is mild): the constants were validated on two microarchitectures to
within a few percent of measured ratios.

---

6. The two-prime FFT layer underneath

The transforms used above come from FLINT's fft_small machinery.
Its relevant properties for powm:

Chunked representation. An integer is split into slots of bits
bits (bits around 40–50 chosen by the plan); multiplication becomes
polynomial multiplication of the slot sequences followed by carry
propagation.

Two-prime CRT (Garner). Slot products can exceed one word, so the
convolution is computed modulo two ~50-bit primes p1, p2 in
double-precision FFTs. Each slot is recovered by Garner's mixed-radix
formula:

z  =  r1  +  p1 * ( (r2 - r1) * p1^(-1)  mod p2 ),
r1 = z mod p1,   r2 = z mod p2,   0 <= z < p1*p2,

evaluated 8 lanes at a time in floating point, followed by a scalar
sweep that shifts each z into its bit position j*bits of the
output and accumulates with carries. The reconstruction is
destructive on the transform lanes (output conversion is destructive
by default throughout fft_small; non-destructive conversions copy
first), and runs as one full vectorized sweep followed by one scalar
sweep — the separated shape that store-forwarding rewards.

Truncated, low, high and window products. The transforms support
bit-granular truncation: only itr slots are transformed and only an
output window [zl, zh) of limbs is exported, with slot bounds
guaranteeing exactness of the window. mulhigh and mullow are the
window specializations used by the Barrett and Montgomery quotient
steps. For signed accumulations the exports use centered
representatives in (-P/2, P/2] with the sign resolved during
recomposition.

Cyclic and negacyclic products. A length-N cyclic convolution
computes multiplication modulo B^(N*bits/64) - 1 with no zero
padding (used by the fold engine); a negacyclic convolution — the same
transform with a weight w^j applied per slot, w^(2N) = 1,
w^N = -1 — computes multiplication modulo B^h + 1 (used at every
chain level). Both avoid the 2x padding of a plain product, which is
the entire point of the fold: the ladder's reductions run in rings
where the transform length matches the modulus size instead of
doubling it.

---

7. The chain: recursive CRT over B^s - 1 = (B^h - 1)(B^h + 1)

powm_chain_init fixes rn, the smallest power of two
>= max(mn, 128), and builds the tower s = rn, rn/2, ..., 64.
At each level, multiplication modulo B^s - 1 splits by CRT
(gcd(B^h - 1, B^h + 1) = 1 for even B, using h = s/2):

x1 = x mod (B^h - 1)      -> recurse to the next level
x2 = x mod (B^h + 1)      -> one negacyclic product against the
                             cached transform of  m mod (B^h + 1)

The negacyclic products use sd_fft_mpn_mulmod_2expp1 with the
modulus-residue transform Fm[lev] precomputed once per level at
setup when h >= 128 (CHAIN_NEG_H), and a basecase
mulmod_2expp1 below that. Residue extraction is linear: modulo
B^h - 1 fold-and-add, modulo B^h + 1 alternate-and-subtract.

Recombination uses the explicit CRT for this coprime pair, in the
same form as flint_mpn_mulmod_bnm1: with u ≡ y (mod B^h - 1) and
v ≡ y (mod B^h + 1),

s  =  (u - v)/2   mod (B^h + 1),
y  =  u + s*B^h - s          mod (B^s - 1),

where s*B^h - s = s*(B^h - 1) vanishes modulo B^h - 1 (preserving
u) and equals -2s ≡ v - u modulo B^h + 1 (correcting to v);
the division by 2 is a shift after a parity fix using
B^h + 1 ≡ 0. The recursion bottoms out at 64 limbs with a plain
product against the cached m mod (B^64 - 1).

The chain's appeal is that all of its per-level modulus transforms
are cached — the per-reduction work is one forward transform, one
pointwise pass and one inverse transform per level on the q side
only, at geometrically decreasing sizes. Its weakness is the
power-of-two rn: at mn just above a power of two the padding is
nearly 2x, which is when the cyclic plan (whose admissible lengths are
much denser) wins — hence the 27:33 selector.

---

8. Where the constant factors come from: an accounting

Per exponent bit at mn = 1030 limbs (measured):

squaring (transform-reuse)         ~ 0.74 mul     81 us
fold reduction:
    quotient step (FFT, cached Fminv)  ~ 0.28 mul
    q*m cyclic (cached Fm)             ~ 0.52 mul
    folds, rotation, carries           ~ 0.10 mul
                                   ---------------
                                   ~ 1.96 mul    209 us

against roughly 3 mul-equivalents per bit for a classical
Barrett ladder (1 squaring + ~2 for the reduction), matching the
observed ~1.3x asymptotic speedup over the previous implementation
and ~3x over mpz_powm at large sizes.

---

9. Remaining optimizations

Ordered roughly by expected value per unit of implementation effort.

1. Cache the transforms of the sliding-window table. Table
   multiplications currently run as full products; transforming each
   T[j] once at table-build time and multiplying pointwise against
   the (already transformed) accumulator would remove one forward
   transform per table multiplication, ~30% of its cost. At large
   exponents table muls are only a few percent of the total, but at
   moderate exponent lengths (a few hundred bits) the table build plus
   its multiplications are a visibly larger share, and the same cached
   transforms would accelerate the build itself
   (T[j+1] = T[j] * T[0]^2 reuses the transform of T[0]^2).

2. Cache the wraparound transforms in the 120–2000 band. The
   mulmod_preinvn_fold basecase calls mpn_mulmod_bnm1 afresh per
   reduction, re-transforming the modulus every time; likewise dinv
   is re-consumed by a plain mulhigh per reduction. A persistent
   B^g - 1 context holding the transforms of d (and, where the
   quotient product is large enough to transform, of dinv) across
   the whole ladder would mirror what the large tier already does with
   Fm/Fminv. This band covers four octaves of sizes and currently
   pays a per-reduction setup the large tier has eliminated.

3. Fuse the squaring with the reduction in the transform domain.
   The fold path today exports the squaring X = acc^2 to limbs, then
   immediately transforms X_lo for the quotient step and q for the
   ledger. The squaring's own transforms are thrown away at export.
   A window-export design could keep acc transformed across the
   square-reduce pair: export only X_lo (a truncated window) for the
   quotient, and feed the ledger from the retained transform where
   lengths allow. This is precisely the chain path's philosophy
   extended into the fold; even a partial fusion (sharing the
   fft(X_lo) between the quotient step and the -X_lo correction)
   removes one transform per bit.

4. Bit-granular chain lengths. The chain's 2x padding cliff at
   mn slightly above a power of two comes from rn being a power of
   two. The negacyclic machinery supports weighted lengths at finer
   granularity (as the cyclic plans already exploit); a chain over
   B^rn - 1 with rn from a denser admissible set would move its
   crossover into territory the cyclic plan currently owns, and make
   the 27:33 selector nearly moot.

5. Signed-window (NAF) recoding. With centered signed exports
   already supported by the FFT layer, a signed sliding window would
   shrink the table by half for the same window width or lengthen the
   effective window for the same table, trading table muls for
   essentially free negations mod m.

6. A negacyclic fold. The ledger ring B^rlen - 1 could equally
   be B^rlen + 1 (the rotation trick works with a sign), and the
   negacyclic transform of the same length covers twice the integer
   size; where admissible cyclic lengths near mn are sparse, the
   negacyclic ring of half the transform length may be cheaper.

7. Base-2 at large sizes. With a base-2 exponent the table entries
   are powers 2^(2j+1) and every table multiplication is a shift;
   the current scalar special covers this, but the squaring chain
   itself could exploit the sparsity of the initial segments (GMP does
   at some sizes; the timing grid shows a remaining 6x gap at
   pb = 32768, eb = 16 for base 2, above the fmpz-level fallback's
   size cutoff).

8. Per-architecture constants. The tier thresholds and the 27:33
   selector are calibrated on two x86-64 microarchitectures and agree
   between them; they live in powm.c pending a flint-mparam.h
   migration if a future architecture (or the zen chain-selector
   observation at sizes moderately below a power of two, where the
   selector oversells the chain by up to ~4%) demands per-arch values.