Add small-value delayed reduction with Barrett algorithm

[wu-s-john]

PR: Add small-value delayed reduction with Barrett algorithm

NOTE: This PR is a small chunk of #98
NOTE: This PR is necessary for adding small-value sumcheck which is in #112
NOTE: A following PR is expected to create more optimized reduction algorithms for Pasta fields since they are large primes in the nice form of 2^l + c where c is small, like two limbs.

Summary

This PR adds infrastructure for delayed modular reduction when accumulating field × small-integer products. This is a prerequisite for the upcoming small-value sumcheck optimization, which is a significantly larger change.

Motivation

In the sumcheck protocol, we frequently accumulate sums of products: Σ(field × value). Previously, we supported delayed reduction only for field × field products:

field × field → Both operands in Montgomery form (factor R²)
             → Use Montgomery REDC (9-limb → 4-limb)
             → Final result automatically in Montgomery form

However, small-value sumcheck uses coefficients that fit in native integers (i32, i64, i128). These are not in Montgomery form:

field × small_int → Field in Montgomery form, small_int is raw
                  → Product has factor R (not R²)
                  → Montgomery REDC would divide by R, giving wrong result

Why Barrett Reduction?

Montgomery reduction computes x × R⁻¹ mod p. When both operands are in Montgomery form, (a·R) × (b·R) = ab·R², and REDC gives ab·R — exactly what we want.

Barrett reduction computes x mod p directly, without any R factor. When accumulating (field·R) × small_int, the product is field·R·small_int. Barrett gives us field·small_int mod p in Montgomery form — exactly right.

The algorithm replaces expensive division by multiplication with a precomputed reciprocal:

μ = ⌊2^512 / p⌋           // precomputed constant

To compute x mod p:
1. q = ⌊(x × μ) / 2^512⌋  // quotient estimate via multiply + shift
2. r = x − q × p           // remainder
3. if r ≥ p: r −= p        // conditional correction

Implementation Efficiency

The implementation is optimal with at most 1 conditional subtract for finalization:

• For BN254 (where 2p < 2^256): The remainder fits in 4 limbs, enabling a fast path with mul_3x4_lo4 and a single if r ≥ p: r -= p
• For fields like T256 (where 2p ≈ 2^256): Uses 5-limb arithmetic with one correction

The bound is tight because:

• Barrett's quotient estimate q satisfies q ≤ ⌊x/p⌋ ≤ q + 2
• Therefore r = x - q·p satisfies 0 ≤ r < 3p
• After one subtraction: 0 ≤ r < 2p
• After at most one more: 0 ≤ r < p ✓

In practice for our field sizes, a single correction suffices (verified by debug_assert!).

Verification of Barrett Constants

All precomputed constants are verified at test time against reference implementations using num-bigint:

Constant	Verification
BARRETT_MU	Computed as ⌊2^512 / p⌋ using BigUint, compared against hardcoded limbs
R384_MOD	Computed as 2^384 mod p using BigUint, compared against hardcoded limbs
USE_4_LIMB_BARRETT	Implicitly verified — if wrong, Barrett reduction would produce non-canonical results caught by debug_assert!

The test helpers in field_reduction_constants.rs implement this:

pub(crate) fn test_barrett_mu_impl<F: BarrettReductionConstants>() {
  let p = limbs_to_biguint(&F::MODULUS);
  let two_pow_512 = BigUint::from(1u64) << 512;
  let expected = &two_pow_512 / &p;
  let actual = limbs5_to_biguint(&F::BARRETT_MU);
  assert_eq!(actual, expected, "BARRETT_MU mismatch");
}

pub(crate) fn test_r384_mod_impl<F: BarrettReductionConstants>() {
  let p = limbs_to_biguint(&F::MODULUS);
  let two_pow_384 = BigUint::from(1u64) << 384;
  let expected = &two_pow_384 % &p;
  let actual = limbs_to_biguint(&F::R384_MOD);
  assert_eq!(actual, expected, "R384_MOD mismatch");
}

Each field provider (BN254, Pasta, P256, T256) runs these tests via macros:

crate::test_barrett_reduction_constants!(scalar_brc, Scalar);
crate::test_barrett_reduction!(scalar_br, Scalar, barrett_reduce_6::<Scalar>);

Additionally, the Barrett reduction itself is tested against field arithmetic:

• Single product: reduce(field × small) == field * F::from(small)
• Sum of products: reduce(Σ field_i × small_i) == Σ (field_i * F::from(small_i))
• Stress test: 2000 random products accumulated and reduced (release builds only)

Future Work: Specialized Pasta Reduction (Already Implemented)

The Pasta curves (Pallas and Vesta) have primes in a special pseudo-Mersenne form that enables significantly faster reduction. From src/provider/pasta.rs:

// Pallas scalar field (Fq):
"40000000000000000000000000000000224698fc0994a8dd8c46eb2100000001"

// Vesta scalar field (Fp):
"40000000000000000000000000000000224698fc094cf91b992d30ed00000001"

Breaking down the structure in 64-bit limbs (little-endian):

         limb[3]            limb[2]            limb[1]            limb[0]
    ┌──────────────────┬──────────────────┬──────────────────┬──────────────────┐
Fq: │ 4000000000000000 │ 0000000000000000 │ 224698fc0994a8dd │ 8c46eb2100000001 │
Fp: │ 4000000000000000 │ 0000000000000000 │ 224698fc094cf91b │ 992d30ed00000001 │
    └──────────────────┴──────────────────┴──────────────────┴──────────────────┘
           ↑                    ↑                         ↑
          2^254            always zero              ε (~125 bits)

Both primes have the form p = 2²⁵⁴ + ε where:

• limb[3] = 0x4000000000000000 contributes exactly 2^62 × 2^192 = 2^254
• limb[2] = 0 (the 128-191 bit range is empty)
• ε = limb[0] + limb[1] × 2^64 spans only ~125 bits

This structure enables Solinas reduction, which exploits the congruence 2^254 ≡ -ε (mod p):

To reduce x (a 5-6 limb value):
1. Split: x = x_lo + x_hi × 2^254       // x_lo is low 254 bits
2. Fold:  r = x_lo - x_hi × ε           // multiply by ~125-bit ε, not 256-bit p
3. Correct: if r < 0 or r ≥ p: adjust

This is faster than generic Barrett because:

Operation	Generic Barrett	Pasta Solinas
Main multiply	3×5 = 15 muls (q × μ)	4×2 = 8 muls (x_hi × ε)
Secondary multiply	3×4 = 12 muls (q × p)	— (folded inline)
Total muls	~27	~8

Note: The specialized pasta_reduce_6 implementation is already complete and tested. It will be included in a follow-up PR after this one is approved and merged, to keep the review scope manageable.

What's Included

Component	Purpose
BarrettReductionConstants	Trait with μ = ⌊2^512/p⌋, R384_MOD, and fast-path flag
barrett_reduce_6/7	Generic Barrett for 6-7 limb inputs
SignedWideLimbs<N>	Accumulator for signed products (separate pos/neg)
SmallValueField<V>	Trait for small integer ↔ field conversion
DelayedReduction<i32/i64/i128>	Accumulation + single reduction

Overflow Capacity

Accumulator	Product Type	Capacity
WideLimbs<6>	field × i32/i64	2^66 products
SignedWideLimbs<7>	field × i128	2^66 products

Sumcheck polynomials are bounded by practical sizes (≤2^40), so overflow is impossible.

Why Split This Out?

The small-value sumcheck PR that uses this infrastructure is substantially larger. Splitting the Barrett reduction foundation into its own PR makes review more manageable and establishes a clean abstraction boundary.

[srinathsetty]

@wu-s-john, here are some bugs found by Copilot.

BUG 1 (Critical) — barrett_reduce_7 drops carry, produces wrong result

File: src/big_num/barrett.rs:125

let (s5, _) = c[5].carrying_add(0, cy); // carry DISCARDED

The fold step computes c[0..6] + c[6] * R384_MOD but discards the carry out of position 5. I constructed a concrete input proving the result is wrong:

For T256 with a 7-limb value near 2 × 2^384 − 1 (c[6]=1, c[0..6]≈2^384−1):

• Correct result: 0xfffffffd000000010000000000000004...
• Code produces: 0xfffffffe000000010000000000000002...
• Error = R384_MOD (off by 2^384 mod p)

This is reachable with ~5 max-magnitude field × i128 products on T256/P256 (where R384_MOD ≈ 2^256). Any sumcheck accumulation with a handful of large same-sign products can trigger this.

Fix: After the fold addition, if there's a carry, add R384_MOD to the 6-limb result (since the dropped carry represents 2^384 ≡ R384_MOD (mod p)):

let (s5, extra_carry) = c[5].carrying_add(0, cy);
let mut folded = [s0, s1, s2, s3, s4, s5];
if extra_carry {
// Dropped carry = 2^384 ≡ R384_MOD (mod p), add it back
let (f, cy2) = add::<6>(&folded, &[F::R384_MOD[0], F::R384_MOD[1], F::R384_MOD[2], F::R384_MOD[3], 0, 0]);
folded = f;
// cy2 can be handled by another R384_MOD addition or absorbed by barrett_reduce_6
}

Alternatively, fold into a 7-limb temporary and call a wider Barrett, or reduce before folding.

---

BUG 2 (Medium) — accumulate_field_times_i128 panics on i128::MIN

File: src/big_num/delayed_reduction.rs:98

(&mut acc.neg, (-*value) as u128)

-i128::MIN overflows in signed i128 arithmetic. In debug builds, this panics. In release, it wraps and happens to give the correct magnitude by coincidence (bit reinterpretation).

The i32 and i64 impls in macros.rs correctly use wrapping_neg(). This one does not.

Fix:

(&mut acc.neg, (*value).wrapping_neg() as u128)

---

BUG 3 (Medium) — small_to_field for i32 panics on i32::MIN, wrong in release

File: src/big_num/macros.rs:130

-Self::from((-val) as u64)

For val = i32::MIN:

• Debug: panics (-i32::MIN overflows)
• Release: -i32::MIN wraps to i32::MIN, then i32::MIN as u64 sign-extends to 0xFFFFFFFF80000000 = 2^64 − 2^31. The result is −(2^64 − 2^31) mod p, but the correct answer is −2^31 mod p. These differ. I verified: i32::MIN as u64 =
  18446744071562067968 (0xFFFFFFFF80000000)
  Correct magnitude = 2147483648 (0x80000000)

The i64 impl (i64_to_field) correctly uses val.wrapping_neg() as u64.

Fix:

-Self::from((val as i64).wrapping_neg() as u64)

[Copilot]

Pull request overview

This PR adds big-number infrastructure to support delayed modular reduction for field × small-integer accumulations using a generic μ-Barrett reduction path, enabling upcoming small-value sumcheck optimizations.

Changes:

• Add BarrettReductionConstants + compile-time constant generation and provider wiring for BN254/Pasta/P256/T256.
• Introduce generic barrett_reduce_6 / barrett_reduce_7 reducers and signed wide accumulators to support delayed reduction of signed small values.
• Add SmallValueField/WideMul helpers and corresponding test macros and provider test coverage.

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
src/provider/pt256.rs	Wires Barrett constants, small-value conversions, delayed reduction, and tests for P256/T256 scalars.
src/provider/pasta.rs	Wires Barrett constants, small-value conversions, delayed reduction, and tests for Pasta scalars.
src/provider/bn254.rs	Wires Barrett constants, small-value conversions, delayed reduction, and tests for BN254 scalar.
src/errors.rs	Adds SmallValueOverflow error for failed small-int conversions.
src/big_num/wide_mul.rs	Adds widening-multiplication trait for small integer products.
src/big_num/small_value_field.rs	Adds SmallValueField trait + helpers and vec_to_small conversion routine.
src/big_num/mod.rs	Registers new modules and re-exports new traits/types.
src/big_num/macros.rs	Adds macros for Barrett constants, small-value field conversions, and delayed reduction impls.
src/big_num/limbs.rs	Adds SignedWideLimbs, sub_mag, and additional limb ops supporting new reducers/accumulators.
src/big_num/field_reduction_constants.rs	Introduces BarrettReductionConstants trait and test helpers for Barrett constants.
src/big_num/delayed_reduction.rs	Adds accumulation helpers for field×u64/u128 and test macros for small-value delayed reduction.
src/big_num/barrett.rs	Implements Barrett reduction for 6- and 7-limb inputs and associated tests/macros.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.