Definitions for GF(2) circular convolution via NTT

The transform below uses the standard toolkit for number-theoretic transforms: a 2³⁰-smooth prime modulus, Montgomery arithmetic for modular multiplication, and a radix-4 decimation. See [Har14] for a general treatment of fast NTT arithmetic.

def mod64 : UInt64

The NTT-friendly prime modulus p = 3·2³⁰ + 1 = 3221225473, as a UInt64. Its 2³⁰-smooth factor p - 1 = 3·2³⁰ provides power-of-two roots of unity for the transform.

Equations

• mod64 = 3 * 2 ^ 30 + 1

def mod32 : UInt32

The same prime modulus p = 3·2³⁰ + 1 = 3221225473 as mod64, but as a UInt32.

Equations

• mod32 = 3 * 2 ^ 30 + 1

def primRoot : UInt64

5 is a primitive root of (ZMod p)ˣ, used to generate the NTT roots of unity.

Equations

• primRoot = 5

def montPprime : UInt32

Montgomery negated inverse p' = -p⁻¹ mod 2³² (= 3221225471), used by montMul.

Equations

• montPprime = 3221225471

def montR2 : UInt32

Montgomery R² mod p for radix R = 2³² (= 1789569709). montMul a montR2 puts a into Montgomery form.

Equations

• montR2 = 1789569709

def montR1 : UInt32

Montgomery R mod p = 2³² mod p (= 1073741823); the Montgomery representation of 1.

Equations

• montR1 = 1073741823

@[inline]

def addMod32 (a b : UInt32) : UInt32

Equations

• addMod32 a b = if a ≥ mod32 - b then a - (mod32 - b) else a + b

@[inline]

def subMod32 (a b : UInt32) : UInt32

Equations

• subMod32 a b = if a ≥ b then a - b else a - b + mod32

@[inline]

def montMul (a b : UInt32) : UInt32

Equations

• One or more equations did not get rendered due to their size.

@[inline]

def toMont (a : UInt32) : UInt32

Equations

• toMont a = montMul a montR2

def powModAuxU64 (mod_ : UInt64) : ℕ → UInt64 → UInt64 → UInt64 → UInt64

Fuel-based square-and-multiply accumulator for modular exponentiation: with fuel rounds left it folds the remaining exponent bits e into the running result r, squaring base b each round. Called by powModU64 with fuel = 64, enough for any UInt64 exponent.

Equations

• powModAuxU64 mod_ 0 x✝² x✝¹ x✝ = x✝
• powModAuxU64 mod_ f.succ x✝² x✝¹ x✝ = if (x✝¹ == 0) = true then x✝ else powModAuxU64 mod_ f (x✝² * x✝² % mod_) (x✝¹ >>> 1) (if (x✝¹ &&& 1 != 0) = true then x✝ * x✝² % mod_ else x✝)

def powModU64 (base exp mod_ : UInt64) : UInt64

base ^ exp mod mod_, computed by square-and-multiply via powModAuxU64.

Equations

• powModU64 base exp mod_ = powModAuxU64 mod_ 64 (base % mod_) exp 1

def montPow (seed wm : UInt32) : ℕ → UInt32

Iterated Montgomery multiplication: apply montMul · wm exactly j times to seed.

Equations

• montPow seed wm 0 = seed
• montPow seed wm j.succ = montMul (montPow seed wm j) wm

def rootsInner (wm : UInt32) (halfLen : ℕ) {n : ℕ} (k i : ℕ) : Vector UInt32 n → Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.
• rootsInner wm halfLen 0 x✝¹ x✝ = x✝

def ensureRoots (n : ℕ) : Vector UInt32 n

roots[2^k + j] = w_{2^(k+1)}^j · R mod mod64 (Montgomery domain)

Equations

• ensureRoots n = if hn0 : n = 0 then #v[] else ensureRoots.outer n (Vector.replicate n 0) 2 ⋯ 64

def ensureRoots.outer (n : ℕ) (v : Vector UInt32 n) (len : UInt64) (hn : 0 < n) : ℕ → Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.
• ensureRoots.outer n v len hn 0 = v

def bitRevNext : ℕ → ℕ → ℕ → ℕ

Compute next index in bit-reveral order ("add 1 in bit-reversed binary")

Equations

• bitRevNext 0 x✝¹ x✝ = x✝¹ ^^^ x✝
• bitRevNext f.succ x✝¹ x✝ = if (x✝¹ &&& x✝ == 0) = true then x✝¹ ^^^ x✝ else bitRevNext f (x✝¹ ^^^ x✝) (x✝ >>> 1)

def bitRevLoop {n : ℕ} (k i : ℕ) : Vector UInt32 n → ℕ → Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.
• bitRevLoop 0 x✝² x✝¹ x✝ = x✝¹

def radix2Pass {n : ℕ} (k i : ℕ) : Vector UInt32 n → Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.
• radix2Pass 0 x✝¹ x✝ = x✝

@[inline]

def butterfly4 {n : ℕ} (a : Vector UInt32 n) (inverse : Bool) (roots : Vector UInt32 n) (s len i2 j2 : ℕ) : Vector UInt32 n

── Radix-4 butterfly ────────────────────────────────────────────────────────

Two-stage DIF decomposition. Stage 1 applies the same twiddle t1 to both aB (at offset s) and aD (at offset len + s) — this is intentional: in this grouping both elements occupy the "odd" slot of their respective half-group and share the same first-stage twiddle. t2 and t3 are the distinct second-stage twiddles applied after the first-stage ±-combine.

Equations

• One or more equations did not get rendered due to their size.

def radix4Inner {n : ℕ} (inverse : Bool) (roots : Vector UInt32 n) (s len i2 : ℕ) : ℕ → ℕ → Vector UInt32 n → Vector UInt32 n

Equations

• radix4Inner inverse roots s len i2 0 x✝¹ x✝ = x✝
• radix4Inner inverse roots s len i2 k.succ x✝¹ x✝ = radix4Inner inverse roots s len i2 k (x✝¹ + 1) (butterfly4 x✝ inverse roots s len i2 x✝¹)

def radix4Middle {n : ℕ} (inverse : Bool) (roots : Vector UInt32 n) (s len : ℕ) : ℕ → ℕ → Vector UInt32 n → Vector UInt32 n

Equations

• radix4Middle inverse roots s len 0 x✝¹ x✝ = x✝
• radix4Middle inverse roots s len k.succ x✝¹ x✝ = radix4Middle inverse roots s len k (x✝¹ + 1) (radix4Inner inverse roots s len (x✝¹ * 2 * len) s 0 x✝)

def nttInplace {n : ℕ} (arr : Vector UInt32 n) (inverse : Bool) (roots : Vector UInt32 n) : Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.

def nttInplace.go : ℕ → ℕ → UInt64 → UInt64

Equations

• nttInplace.go 0 a✝¹ a✝ = a✝
• nttInplace.go f.succ a✝¹ a✝ = if a✝¹ ≤ 1 then a✝ else nttInplace.go f (a✝¹ >>> 1) (a✝ + 1)

def nttInplace.outerLoop {n : ℕ} (inverse : Bool) (roots a : Vector UInt32 n) (len : UInt64) : ℕ → Vector UInt32 n

Equations

• One or more equations did not get rendered due to their size.
• nttInplace.outerLoop inverse roots a len 0 = a

def circularConvolutionGf2 {n : ℕ} (a b : BitVec n) : BitVec n

GF(2) circular convolution: (a * b)[i] = ⊕ₖ (a[k] ∧ b[(i−k) mod n]). Equivalently, multiplication in GF(2)[X] / (Xⁿ − 1). Computed via radix-4 NTT with Montgomery arithmetic. Uses 2n-point NTT to avoid aliasing; folds result[i] + result[i+n] to recover the circular (period-n) convolution from the linear product.

This is only correct for n < 2 ^ 29 ≈ 2.6e8.

Equations

• One or more equations did not get rendered due to their size.

def BitVec.circConvolutionBruteforce {n : ℕ} (v1 v2 : BitVec n) : BitVec n

Specification of circular convolution on bit vectors. Note that for Bool, summing is xor, while multiplication is and.

Equations

• v1.circConvolutionBruteforce v2 = BitVec.ofFnLE fun (i : Fin n) => ∑ j : Fin n, v1[j] * v2[i - j]