Montgomery reduction: for input x, returns y where y ≡ x*R^(-1) (mod q). This is a generic implementation parameterized by the modulus q, its inverse qInv, the Montgomery constant R, and the result bound.
For ML-DSA: R = 2^32, returns y < 2q For ML-KEM: R = 2^16, returns y in range (-q, q)
fn montgomeryReduce(
comptime InT: type,
comptime OutT: type,
comptime q: comptime_int,
comptime qInv: comptime_int,
comptime r_bits: comptime_int,
x: InT,
) OutT
fn montgomeryReduce(
comptime InT: type,
comptime OutT: type,
comptime q: comptime_int,
comptime qInv: comptime_int,
comptime r_bits: comptime_int,
x: InT,
) OutT {
const mask = (@as(InT, 1) << r_bits) - 1;
const m_full = (x *% qInv) & mask;
const m: OutT = @truncate(m_full);
const yR = x -% @as(InT, m) * @as(InT, q);
const y_shifted = @as(@Int(.unsigned, @typeInfo(InT).Int.bits), @bitCast(yR)) >> r_bits;
return @bitCast(@as(@Int(.unsigned, @typeInfo(OutT).Int.bits), @truncate(y_shifted)));
}