Merge pull request #87 from alxiong/lazy-mont-red

Expose deferred Product result for lazy reduction

Cargo.toml

@@ -71,4 +71,5 @@ gpu = ["alloc", "ec-gpu"]
 sqrt-table = ["alloc", "lazy_static"]
 repr-c = []
 uninline-portable = []
+deferred = []
 serde = ["hex", "serde_crate"]

src/deferred.rs

@@ -0,0 +1,284 @@
+//! Deferred normalization for field arithmetic.
+//!
+//! This module provides the [`DeferredField`] trait and a wide [`Product`]
+//! accumulator. Together they enable accumulating multiple unreduced
+//! Montgomery products before performing a single expensive reduction.
+//! This is useful for operations like inner products where many
+//! multiplications feed into a sum.
+
+use core::fmt::Debug;
+
+use crate::arithmetic::{adc, mac};
+
+/// A trait for fields that support deferred reduction of products.
+///
+/// Instead of reducing each multiplication result immediately, callers
+/// accumulate products into an [`Accumulator`](Self::Accumulator) via
+/// [`mul_accumulate`](Self::mul_accumulate) and
+/// [`square_accumulate`](Self::square_accumulate), then perform a single
+/// reduction at the end with [`reduce`](Self::reduce).
+pub trait DeferredField: ff::Field {
+    /// A wide accumulator for unreduced products.
+    type Accumulator: Copy + Clone + Debug + Default;
+
+    /// Multiplies `a` by `b` and adds the result into `acc`.
+    fn mul_accumulate(acc: &mut Self::Accumulator, a: &Self, b: &Self);
+
+    /// Squares `a` and adds the result into `acc`.
+    fn square_accumulate(acc: &mut Self::Accumulator, a: &Self);
+
+    /// Reduces the accumulator to a canonical field element.
+    fn reduce(acc: Self::Accumulator) -> Self;
+}
+
+/// A wide accumulator for unreduced Montgomery products over field `F`.
+///
+/// This stores a running sum of 512-bit products with a 64-bit carry for
+/// overflow beyond 512 bits. Products are added internally by
+/// [`DeferredField::mul_accumulate`] and [`DeferredField::square_accumulate`].
+///
+/// Call [`DeferredField::reduce`] to fold the carry back into range and
+/// perform Montgomery reduction.
+#[derive(Clone, Copy, Debug)]
+pub struct Product<F> {
+    limbs: [u64; 8],
+    carry: u64,
+    _marker: core::marker::PhantomData<F>,
+}
+
+impl<F> Default for Product<F> {
+    fn default() -> Self {
+        Self::ZERO
+    }
+}
+
+impl<F> Product<F> {
+    /// The zero (additive identity) accumulator.
+    pub const ZERO: Self = Product {
+        limbs: [0; 8],
+        carry: 0,
+        _marker: core::marker::PhantomData,
+    };
+
+    /// Adds a raw 512-bit product (8 limbs) into this accumulator.
+    ///
+    /// Each call contributes at most 1 to `carry`; overflow of the 64-bit
+    /// carry requires 2^64 accumulated products (~590 exabytes of input).
+    #[inline]
+    pub(crate) fn accumulate(&mut self, product: [u64; 8]) {
+        let (d0, c) = adc(self.limbs[0], product[0], 0);
+        let (d1, c) = adc(self.limbs[1], product[1], c);
+        let (d2, c) = adc(self.limbs[2], product[2], c);
+        let (d3, c) = adc(self.limbs[3], product[3], c);
+        let (d4, c) = adc(self.limbs[4], product[4], c);
+        let (d5, c) = adc(self.limbs[5], product[5], c);
+        let (d6, c) = adc(self.limbs[6], product[6], c);
+        let (d7, c) = adc(self.limbs[7], product[7], c);
+        self.limbs = [d0, d1, d2, d3, d4, d5, d6, d7];
+        let (carry, overflow) = self.carry.overflowing_add(c);
+        debug_assert!(!overflow, "carry overflow: too many accumulated products");
+        self.carry = carry;
+    }
+
+    /// Folds `carry` (bits 512+) and `limbs[7]` (bits 448–511) into the lower
+    /// 448 bits using precomputed residues of $2^{448}$ and $2^{512}$ modulo
+    /// the field prime.
+    ///
+    /// The result fits in 8 limbs with value $< 2^{449} < Rp$, safe for
+    /// Montgomery reduction.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) fn partial_reduce(&self, b448: &[u64; 4], r2: &[u64; 4]) -> [u64; 8] {
+        let b7 = self.limbs[7];
+        let b8 = self.carry;
+
+        // Compute b7 * b448 (5 limbs)
+        let (t0, c) = mac(0, b7, b448[0], 0);
+        let (t1, c) = mac(0, b7, b448[1], c);
+        let (t2, c) = mac(0, b7, b448[2], c);
+        let (t3, c) = mac(0, b7, b448[3], c);
+        let t4 = c;
+
+        // Accumulate b8 * r2
+        let (t0, c) = mac(t0, b8, r2[0], 0);
+        let (t1, c) = mac(t1, b8, r2[1], c);
+        let (t2, c) = mac(t2, b8, r2[2], c);
+        let (t3, c) = mac(t3, b8, r2[3], c);
+        let (t4, t5) = adc(t4, 0, c);
+        debug_assert!(
+            t5 == 0,
+            "folding term overflow: t4 + carry does not fit in 64 bits"
+        );
+
+        // Add to lower 7 limbs
+        let (d0, c) = adc(self.limbs[0], t0, 0);
+        let (d1, c) = adc(self.limbs[1], t1, c);
+        let (d2, c) = adc(self.limbs[2], t2, c);
+        let (d3, c) = adc(self.limbs[3], t3, c);
+        let (d4, c) = adc(self.limbs[4], t4, c);
+        let (d5, c) = adc(self.limbs[5], 0, c);
+        let (d6, c) = adc(self.limbs[6], 0, c);
+        let (d7, _) = adc(0, 0, c);
+
+        // B448 < 2^253 and r2 < 2^252, so the folding term
+        // b7 * B448 + b8 * r2 < 2^317 + 2^316 < 2^318.
+        // The full value is < 2^448 + 2^318 < 2^449, so d7 is at most 1.
+        debug_assert!(d7 <= 1);
+
+        [d0, d1, d2, d3, d4, d5, d6, d7]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::DeferredField;
+    use ff::Field;
+    use rand::SeedableRng;
+    use rand_xorshift::XorShiftRng;
+    use std::vec::Vec;
+
+    const SEED: [u8; 16] = [
+        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
+        0xe5,
+    ];
+
+    fn inner_product<F: DeferredField>(a: &[F], b: &[F]) -> F {
+        let mut acc = F::Accumulator::default();
+        for (x, y) in a.iter().zip(b.iter()) {
+            F::mul_accumulate(&mut acc, x, y);
+        }
+        F::reduce(acc)
+    }
+
+    macro_rules! deferred_field_tests {
+        ($F:ty, $mod:ident, $adversarial_a:expr, $adversarial_b:expr) => {
+            mod $mod {
+                use super::*;
+
+                #[test]
+                fn accumulate_roundtrip() {
+                    let mut rng = XorShiftRng::from_seed(SEED);
+                    for _ in 0..100 {
+                        let a = <$F>::random(&mut rng);
+                        let b = <$F>::random(&mut rng);
+                        let mut acc = <$F as DeferredField>::Accumulator::default();
+                        <$F>::mul_accumulate(&mut acc, &a, &b);
+                        assert_eq!(<$F>::reduce(acc), a * b);
+                    }
+                }
+
+                #[test]
+                fn square_accumulate_roundtrip() {
+                    let mut rng = XorShiftRng::from_seed(SEED);
+                    for _ in 0..100 {
+                        let a = <$F>::random(&mut rng);
+                        let mut acc = <$F as DeferredField>::Accumulator::default();
+                        <$F>::square_accumulate(&mut acc, &a);
+                        assert_eq!(<$F>::reduce(acc), a.square());
+                    }
+                }
+
+                #[test]
+                fn test_inner_product() {
+                    let mut rng = XorShiftRng::from_seed(SEED);
+                    for len in [0, 1, 2, 3, 4, 7, 8, 15, 16, 31, 32, 100, 255, 256, 10_000] {
+                        let a: Vec<$F> = (0..len).map(|_| <$F>::random(&mut rng)).collect();
+                        let b: Vec<$F> = (0..len).map(|_| <$F>::random(&mut rng)).collect();
+
+                        let eager: $F = a.iter().zip(b.iter()).map(|(x, y)| *x * *y).sum();
+                        let lazy = inner_product(&a, &b);
+
+                        assert_eq!(eager, lazy, "mismatch at len={len}");
+                    }
+                }
+
+                #[test]
+                fn reduce_zero() {
+                    assert_eq!(
+                        <$F>::reduce(<$F as DeferredField>::Accumulator::default()),
+                        <$F>::ZERO,
+                    );
+                }
+
+                #[test]
+                fn square_vs_mul() {
+                    let mut rng = XorShiftRng::from_seed(SEED);
+                    for _ in 0..100 {
+                        let a = <$F>::random(&mut rng);
+                        let mut sq_acc = <$F as DeferredField>::Accumulator::default();
+                        <$F>::square_accumulate(&mut sq_acc, &a);
+                        let mut mul_acc = <$F as DeferredField>::Accumulator::default();
+                        <$F>::mul_accumulate(&mut mul_acc, &a, &a);
+                        assert_eq!(
+                            <$F>::reduce(sq_acc),
+                            <$F>::reduce(mul_acc),
+                            "square_accumulate and mul_accumulate(a, a) diverged",
+                        );
+                    }
+                }
+
+                #[test]
+                fn mixed_accumulate() {
+                    let mut rng = XorShiftRng::from_seed(SEED);
+                    for _ in 0..20 {
+                        let a = <$F>::random(&mut rng);
+                        let b = <$F>::random(&mut rng);
+                        let c = <$F>::random(&mut rng);
+                        let mut acc = <$F as DeferredField>::Accumulator::default();
+                        <$F>::mul_accumulate(&mut acc, &a, &b);
+                        <$F>::square_accumulate(&mut acc, &c);
+                        assert_eq!(<$F>::reduce(acc), a * b + c.square());
+                    }
+                }
+
+                /// Regression: elements with top limb ~0x3F whose products have
+                /// limbs[7] ~0x0F. These adversarial elements exercise the
+                /// partial-reduction path in the lazy Product accumulator.
+                #[test]
+                fn regression_overflow() {
+                    let a = $adversarial_a;
+                    let b = $adversarial_b;
+                    let a_arr = [a; 100];
+                    let b_arr = [b; 100];
+
+                    let eager: $F = a_arr.iter().zip(b_arr.iter()).map(|(x, y)| *x * *y).sum();
+                    let lazy = inner_product(&a_arr, &b_arr);
+
+                    assert_eq!(eager, lazy, "inner_product returned non-canonical result");
+                }
+            }
+        };
+    }
+
+    deferred_field_tests!(
+        crate::Fp,
+        fp,
+        crate::Fp([
+            0x0361524c2cc0f859u64,
+            0xae68690a78bc7175,
+            0xe66cd36e68ef8f5f,
+            0x3fa6524a713b7e05,
+        ]),
+        crate::Fp([
+            0x7a1c5e3b9d204f61u64,
+            0xc48e0b71a2d5f389,
+            0xd9f247a0856c13be,
+            0x3d8a19f5e6c7b042,
+        ])
+    );
+    deferred_field_tests!(
+        crate::Fq,
+        fq,
+        crate::Fq([
+            0x31d0b6640589f877u64,
+            0xf87f43fdf6062541,
+            0xb7d6467b2f5a522a,
+            0x3eb025240950fd13,
+        ]),
+        crate::Fq([
+            0x5e9a3c71f8b20d46u64,
+            0xa3d1e6f504879c2b,
+            0xcb45a8d2e1f36790,
+            0x3c47d2a8b10e5f93,
+        ])
+    );
+}

src/fields/fp.rs

@@ -12,6 +12,8 @@ use lazy_static::lazy_static;
 use ff::{FieldBits, PrimeFieldBits};
 
 use crate::arithmetic::{adc, mac, sbb, SqrtTableHelpers};
+#[cfg(feature = "deferred")]
+use crate::deferred::{DeferredField, Product};
 
 #[cfg(feature = "sqrt-table")]
 use crate::arithmetic::SqrtTables;
@@ -308,33 +310,8 @@ impl Fp {
     /// Squares this element.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn square(&self) -> Fp {
-        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let u = self.square_unreduced();
-        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        Fp::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
-        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
-        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
-
-        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
-
-        let r7 = r6 >> 63;
-        let r6 = (r6 << 1) | (r5 >> 63);
-        let r5 = (r5 << 1) | (r4 >> 63);
-        let r4 = (r4 << 1) | (r3 >> 63);
-        let r3 = (r3 << 1) | (r2 >> 63);
-        let r2 = (r2 << 1) | (r1 >> 63);
-        let r1 = r1 << 1;
-
-        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
-        let (r1, carry) = adc(0, r1, carry);
-        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
-        let (r3, carry) = adc(0, r3, carry);
-        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
-        let (r5, carry) = adc(0, r5, carry);
-        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
-        let (r7, _) = adc(0, r7, carry);
-
-        Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
     }
 
     #[allow(clippy::too_many_arguments)]
@@ -388,29 +365,8 @@ impl Fp {
     /// Multiplies `rhs` by `self`, returning the result.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn mul(&self, rhs: &Self) -> Self {
-        // Schoolbook multiplication
+        let u = self.mul_unreduced(rhs);
-
+        Fp::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
-        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
-        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
-        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
-
-        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
-        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
-        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
-        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
-
-        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
-        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
-        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
-        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
-        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
-        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
-        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
-
-        Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
     }
 
     /// Subtracts `rhs` from `self`, returning the result.
@@ -461,6 +417,96 @@ impl Fp {
 
         Fp([d0 & mask, d1 & mask, d2 & mask, d3 & mask])
     }
+
+    /// Multiplies `rhs` by `self`, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn mul_unreduced(&self, rhs: &Self) -> [u64; 8] {
+        // Schoolbook multiplication
+
+        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
+        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
+        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
+
+        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
+        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
+        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
+        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
+
+        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
+        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
+        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
+        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
+        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
+        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
+        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+
+    /// Squares this element, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn square_unreduced(&self) -> [u64; 8] {
+        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
+        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
+
+        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
+
+        let r7 = r6 >> 63;
+        let r6 = (r6 << 1) | (r5 >> 63);
+        let r5 = (r5 << 1) | (r4 >> 63);
+        let r4 = (r4 << 1) | (r3 >> 63);
+        let r3 = (r3 << 1) | (r2 >> 63);
+        let r2 = (r2 << 1) | (r1 >> 63);
+        let r1 = r1 << 1;
+
+        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
+        let (r1, carry) = adc(0, r1, carry);
+        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
+        let (r3, carry) = adc(0, r3, carry);
+        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
+        let (r5, carry) = adc(0, r5, carry);
+        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
+        let (r7, _) = adc(0, r7, carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+}
+
+#[cfg(feature = "deferred")]
+impl DeferredField for Fp {
+    type Accumulator = Product<Fp>;
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn mul_accumulate(acc: &mut Self::Accumulator, a: &Fp, b: &Fp) {
+        acc.accumulate(a.mul_unreduced(b));
+    }
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn square_accumulate(acc: &mut Self::Accumulator, a: &Fp) {
+        acc.accumulate(a.square_unreduced());
+    }
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn reduce(acc: Self::Accumulator) -> Fp {
+        /// 2^448 mod p (little-endian limbs).
+        const B448: [u64; 4] = [
+            0x9b9858f294cf91ba,
+            0x8635bd2c4252b065,
+            0x496d41af7b9cb714,
+            0x1b4b3c4bfffffffc,
+        ];
+        let limbs = acc.partial_reduce(&B448, &R2.0);
+        Fp::montgomery_reduce(
+            limbs[0], limbs[1], limbs[2], limbs[3], limbs[4], limbs[5], limbs[6], limbs[7],
+        )
+    }
 }
 
 impl From<Fp> for [u8; 32] {

src/fields/fq.rs

@@ -12,6 +12,8 @@ use lazy_static::lazy_static;
 use ff::{FieldBits, PrimeFieldBits};
 
 use crate::arithmetic::{adc, mac, sbb, SqrtTableHelpers};
+#[cfg(feature = "deferred")]
+use crate::deferred::{DeferredField, Product};
 
 #[cfg(feature = "sqrt-table")]
 use crate::arithmetic::SqrtTables;
@@ -308,33 +310,8 @@ impl Fq {
     /// Squares this element.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn square(&self) -> Fq {
-        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let u = self.square_unreduced();
-        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        Fq::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
-        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
-        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
-
-        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
-
-        let r7 = r6 >> 63;
-        let r6 = (r6 << 1) | (r5 >> 63);
-        let r5 = (r5 << 1) | (r4 >> 63);
-        let r4 = (r4 << 1) | (r3 >> 63);
-        let r3 = (r3 << 1) | (r2 >> 63);
-        let r2 = (r2 << 1) | (r1 >> 63);
-        let r1 = r1 << 1;
-
-        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
-        let (r1, carry) = adc(0, r1, carry);
-        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
-        let (r3, carry) = adc(0, r3, carry);
-        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
-        let (r5, carry) = adc(0, r5, carry);
-        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
-        let (r7, _) = adc(0, r7, carry);
-
-        Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
     }
 
     #[allow(clippy::too_many_arguments)]
@@ -388,29 +365,8 @@ impl Fq {
     /// Multiplies `rhs` by `self`, returning the result.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn mul(&self, rhs: &Self) -> Self {
-        // Schoolbook multiplication
+        let u = self.mul_unreduced(rhs);
-
+        Fq::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
-        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
-        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
-        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
-
-        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
-        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
-        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
-        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
-
-        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
-        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
-        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
-        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
-        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
-        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
-        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
-
-        Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
     }
 
     /// Subtracts `rhs` from `self`, returning the result.
@@ -461,6 +417,96 @@ impl Fq {
 
         Fq([d0 & mask, d1 & mask, d2 & mask, d3 & mask])
     }
+
+    /// Multiplies `rhs` by `self`, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn mul_unreduced(&self, rhs: &Self) -> [u64; 8] {
+        // Schoolbook multiplication
+
+        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
+        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
+        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
+
+        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
+        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
+        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
+        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
+
+        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
+        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
+        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
+        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
+        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
+        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
+        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+
+    /// Squares this element, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn square_unreduced(&self) -> [u64; 8] {
+        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
+        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
+
+        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
+
+        let r7 = r6 >> 63;
+        let r6 = (r6 << 1) | (r5 >> 63);
+        let r5 = (r5 << 1) | (r4 >> 63);
+        let r4 = (r4 << 1) | (r3 >> 63);
+        let r3 = (r3 << 1) | (r2 >> 63);
+        let r2 = (r2 << 1) | (r1 >> 63);
+        let r1 = r1 << 1;
+
+        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
+        let (r1, carry) = adc(0, r1, carry);
+        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
+        let (r3, carry) = adc(0, r3, carry);
+        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
+        let (r5, carry) = adc(0, r5, carry);
+        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
+        let (r7, _) = adc(0, r7, carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+}
+
+#[cfg(feature = "deferred")]
+impl DeferredField for Fq {
+    type Accumulator = Product<Fq>;
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn mul_accumulate(acc: &mut Self::Accumulator, a: &Fq, b: &Fq) {
+        acc.accumulate(a.mul_unreduced(b));
+    }
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn square_accumulate(acc: &mut Self::Accumulator, a: &Fq) {
+        acc.accumulate(a.square_unreduced());
+    }
+
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    fn reduce(acc: Self::Accumulator) -> Fq {
+        /// 2^448 mod q (little-endian limbs).
+        const B448: [u64; 4] = [
+            0xcc920bb9994a8dd9,
+            0x87a7dcbe1ff6e0d7,
+            0x496d41af7ccfdaa9,
+            0x0ee4537bfffffffc,
+        ];
+        let limbs = acc.partial_reduce(&B448, &R2.0);
+        Fq::montgomery_reduce(
+            limbs[0], limbs[1], limbs[2], limbs[3], limbs[4], limbs[5], limbs[6], limbs[7],
+        )
+    }
 }
 
 impl From<Fq> for [u8; 32] {

src/lib.rs

@@ -22,6 +22,9 @@ mod curves;
 mod fields;
 
 pub mod arithmetic;
+#[cfg(feature = "deferred")]
+#[cfg_attr(docsrs, doc(cfg(feature = "deferred")))]
+pub mod deferred;
 pub mod pallas;
 pub mod vesta;