rhdl_ed25519_fast_core/

scalar_parallel.rs

//! Initiation-interval-one scalar arithmetic for Ed25519.
//!
//! The group order is `l = 2^252 + q`, where `q` is only 125 bits. Reduction
//! therefore folds high bits through multiplication by the fixed eight-limb
//! radix-`2^16` constant `q`. The multiplier below is fully spatial: each of 24
//! coefficients owns eight registered `16 x 16` products, and normalization
//! advances two limbs per stage.

use rhdl::prelude::*;
use rhdl_fpga::core::dff::DFF;

use crate::{
    protocol::Words512,
    scalar::{ScalarWords, scalar_ge_kernel, scalar_l_kernel, scalar_sub_kernel},
};

/// Seventeen radix-`2^16` input limbs, least significant first.
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct Radix16X17 {
    /// Input limbs.
    pub limbs: [b16; 17],
}

/// Twenty-five radix-`2^16` product limbs, least significant first.
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct Radix16X25 {
    /// Product limbs.
    pub limbs: [b16; 25],
}

/// Tagged input to [`ConstQMulPipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ConstQMulIn {
    /// Marks this transaction valid.
    pub valid: bool,
    /// Multiplicand.
    pub x: Radix16X17,
    /// Caller-defined transaction tag.
    pub tag: b16,
}

/// Tagged output from [`ConstQMulPipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ConstQMulOut {
    /// Marks this transaction valid.
    pub valid: bool,
    /// Exact product `x * q`.
    pub product: Radix16X25,
    /// Delayed caller tag.
    pub tag: b16,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct ConstQCellIn {
    x: [b16; 8],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct ConstQCellOut {
    coefficient: b35,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct ConstQCellRegs {
    x: [b16; 8],
    products: [b32; 8],
    pairs: [b33; 4],
    quads: [b34; 2],
    coefficient: b35,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[doc(hidden)]
pub struct ConstQCell {
    regs: DFF<ConstQCellRegs>,
}

impl Default for ConstQCell {
    fn default() -> Self {
        Self {
            regs: DFF::new(ConstQCellRegs {
                x: [b16(0); 8],
                products: [b32(0); 8],
                pairs: [b33(0); 4],
                quads: [b34(0); 2],
                coefficient: b35(0),
            }),
        }
    }
}

impl SynchronousIO for ConstQCell {
    type I = ConstQCellIn;
    type O = ConstQCellOut;
    type Kernel = const_q_cell_kernel;
}

/// Pipeline transition for one constant-product coefficient.
#[kernel]
#[doc(hidden)]
pub fn const_q_cell_kernel(
    _cr: ClockReset,
    i: ConstQCellIn,
    q: ConstQCellQ,
) -> (ConstQCellOut, ConstQCellD) {
    let constants = [
        b16(0xd3ed),
        b16(0x5cf5),
        b16(0x631a),
        b16(0x5812),
        b16(0x9cd6),
        b16(0xa2f7),
        b16(0xf9de),
        b16(0x14de),
    ];
    let r = q.regs;
    let mut n = r;
    n.x = i.x;
    for index in 0..8 {
        n.products[index] = r.x[index].resize::<32>() * constants[index].resize::<32>();
    }
    for index in 0..4 {
        n.pairs[index] =
            r.products[index * 2].resize::<33>() + r.products[index * 2 + 1].resize::<33>();
    }
    n.quads[0] = r.pairs[0].resize::<34>() + r.pairs[1].resize::<34>();
    n.quads[1] = r.pairs[2].resize::<34>() + r.pairs[3].resize::<34>();
    n.coefficient = r.quads[0].resize::<35>() + r.quads[1].resize::<35>();
    (
        ConstQCellOut {
            coefficient: r.coefficient,
        },
        ConstQCellD { regs: n },
    )
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct ConstQBankIn {
    windows: [[b16; 8]; 4],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct ConstQBankOut {
    coefficients: [b35; 4],
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
#[doc(hidden)]
pub struct ConstQBank {
    cell0: ConstQCell,
    cell1: ConstQCell,
    cell2: ConstQCell,
    cell3: ConstQCell,
}

impl SynchronousIO for ConstQBank {
    type I = ConstQBankIn;
    type O = ConstQBankOut;
    type Kernel = const_q_bank_kernel;
}

/// Wiring transition for four coefficient cells.
#[kernel]
#[doc(hidden)]
pub fn const_q_bank_kernel(
    _cr: ClockReset,
    i: ConstQBankIn,
    q: ConstQBankQ,
) -> (ConstQBankOut, ConstQBankD) {
    let output = ConstQBankOut {
        coefficients: [
            q.cell0.coefficient,
            q.cell1.coefficient,
            q.cell2.coefficient,
            q.cell3.coefficient,
        ],
    };
    let mut d = ConstQBankD::dont_care();
    d.cell0 = ConstQCellIn { x: i.windows[0] };
    d.cell1 = ConstQCellIn { x: i.windows[1] };
    d.cell2 = ConstQCellIn { x: i.windows[2] };
    d.cell3 = ConstQCellIn { x: i.windows[3] };
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct NormalizeState {
    coefficients: [b35; 24],
    limbs: [b16; 24],
    carry: b20,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct NormalizePairIn {
    state: NormalizeState,
    start: b6,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct NormalizePairOut {
    state: NormalizeState,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[doc(hidden)]
pub struct NormalizePair {
    state: DFF<NormalizeState>,
}

impl Default for NormalizePair {
    fn default() -> Self {
        Self {
            state: DFF::new(NormalizeState {
                coefficients: [b35(0); 24],
                limbs: [b16(0); 24],
                carry: b20(0),
            }),
        }
    }
}

impl SynchronousIO for NormalizePair {
    type I = NormalizePairIn;
    type O = NormalizePairOut;
    type Kernel = normalize_pair_kernel;
}

/// Two-limb carry-normalization transition.
#[kernel]
#[doc(hidden)]
pub fn normalize_pair_kernel(
    _cr: ClockReset,
    i: NormalizePairIn,
    q: NormalizePairQ,
) -> (NormalizePairOut, NormalizePairD) {
    let mut next = i.state;
    let first = i.state.coefficients[i.start].resize::<36>() + i.state.carry.resize::<36>();
    let second_index = i.start + b6(1);
    let second = i.state.coefficients[second_index].resize::<36>() + (first >> 16);
    next.limbs[i.start] = first.resize::<16>();
    next.limbs[second_index] = second.resize::<16>();
    next.carry = (second >> 16).resize::<20>();
    (
        NormalizePairOut { state: q.state },
        NormalizePairD { state: next },
    )
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct NormalizeGroupIn {
    state: NormalizeState,
    start: b6,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct NormalizeGroupOut {
    state: NormalizeState,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
#[doc(hidden)]
pub struct NormalizeGroup {
    pair0: NormalizePair,
    pair1: NormalizePair,
    pair2: NormalizePair,
    pair3: NormalizePair,
}

impl SynchronousIO for NormalizeGroup {
    type I = NormalizeGroupIn;
    type O = NormalizeGroupOut;
    type Kernel = normalize_group_kernel;
}

/// Wiring transition for four registered normalization steps.
#[kernel]
#[doc(hidden)]
pub fn normalize_group_kernel(
    _cr: ClockReset,
    i: NormalizeGroupIn,
    q: NormalizeGroupQ,
) -> (NormalizeGroupOut, NormalizeGroupD) {
    let mut d = NormalizeGroupD::dont_care();
    d.pair0 = NormalizePairIn {
        state: i.state,
        start: i.start,
    };
    d.pair1 = NormalizePairIn {
        state: q.pair0.state,
        start: i.start + b6(2),
    };
    d.pair2 = NormalizePairIn {
        state: q.pair1.state,
        start: i.start + b6(4),
    };
    d.pair3 = NormalizePairIn {
        state: q.pair2.state,
        start: i.start + b6(6),
    };
    (
        NormalizeGroupOut {
            state: q.pair3.state,
        },
        d,
    )
}

/// Seventeen-stage multiplication by the Ed25519 reduction constant `q`.
///
/// The pipeline accepts one multiplication every cycle and contains 192
/// independent `16 x 16` multiply sites.
#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[rhdl(dq_no_prefix)]
pub struct ConstQMulPipe {
    bank0: ConstQBank,
    bank1: ConstQBank,
    bank2: ConstQBank,
    bank3: ConstQBank,
    bank4: ConstQBank,
    bank5: ConstQBank,
    normalize0: NormalizeGroup,
    normalize1: NormalizeGroup,
    normalize2: NormalizeGroup,
    valid: DFF<[bool; 17]>,
    tags: DFF<[b16; 17]>,
}

impl Default for ConstQMulPipe {
    fn default() -> Self {
        Self {
            bank0: ConstQBank::default(),
            bank1: ConstQBank::default(),
            bank2: ConstQBank::default(),
            bank3: ConstQBank::default(),
            bank4: ConstQBank::default(),
            bank5: ConstQBank::default(),
            normalize0: NormalizeGroup::default(),
            normalize1: NormalizeGroup::default(),
            normalize2: NormalizeGroup::default(),
            valid: DFF::new([false; 17]),
            tags: DFF::new([b16(0); 17]),
        }
    }
}

impl SynchronousIO for ConstQMulPipe {
    type I = ConstQMulIn;
    type O = ConstQMulOut;
    type Kernel = const_q_mul_kernel;
}

/// State transition and hierarchy wiring for [`ConstQMulPipe`].
#[kernel]
pub fn const_q_mul_kernel(_cr: ClockReset, i: ConstQMulIn, q: Q) -> (ConstQMulOut, D) {
    let bank_coefficients = [
        q.bank0.coefficients,
        q.bank1.coefficients,
        q.bank2.coefficients,
        q.bank3.coefficients,
        q.bank4.coefficients,
        q.bank5.coefficients,
    ];
    let mut coefficients = [b35(0); 24];
    let mut windows = [[[b16(0); 8]; 4]; 6];
    for bank in 0..6 {
        for slot in 0..4 {
            let coefficient = bank * 4 + slot;
            coefficients[coefficient] = bank_coefficients[bank][slot];
            for x_limb in 0..17 {
                for q_limb in 0..8 {
                    if coefficient == x_limb + q_limb {
                        windows[bank][slot][q_limb] = i.x.limbs[x_limb];
                    }
                }
            }
        }
    }
    let initial = NormalizeState {
        coefficients,
        limbs: [b16(0); 24],
        carry: b20(0),
    };
    let product_limbs = q.normalize2.state.limbs;
    let product = Radix16X25 {
        limbs: [
            product_limbs[0],
            product_limbs[1],
            product_limbs[2],
            product_limbs[3],
            product_limbs[4],
            product_limbs[5],
            product_limbs[6],
            product_limbs[7],
            product_limbs[8],
            product_limbs[9],
            product_limbs[10],
            product_limbs[11],
            product_limbs[12],
            product_limbs[13],
            product_limbs[14],
            product_limbs[15],
            product_limbs[16],
            product_limbs[17],
            product_limbs[18],
            product_limbs[19],
            product_limbs[20],
            product_limbs[21],
            product_limbs[22],
            product_limbs[23],
            q.normalize2.state.carry.resize::<16>(),
        ],
    };
    let output = ConstQMulOut {
        valid: q.valid[16],
        product,
        tag: q.tags[16],
    };
    let mut d = D::dont_care();
    d.bank0 = ConstQBankIn {
        windows: windows[0],
    };
    d.bank1 = ConstQBankIn {
        windows: windows[1],
    };
    d.bank2 = ConstQBankIn {
        windows: windows[2],
    };
    d.bank3 = ConstQBankIn {
        windows: windows[3],
    };
    d.bank4 = ConstQBankIn {
        windows: windows[4],
    };
    d.bank5 = ConstQBankIn {
        windows: windows[5],
    };
    d.normalize0 = NormalizeGroupIn {
        state: initial,
        start: b6(0),
    };
    d.normalize1 = NormalizeGroupIn {
        state: q.normalize0.state,
        start: b6(8),
    };
    d.normalize2 = NormalizeGroupIn {
        state: q.normalize1.state,
        start: b6(16),
    };
    d.valid[0] = i.valid;
    d.tags[0] = i.tag;
    for stage in 1..17 {
        d.valid[stage] = q.valid[stage - 1];
        d.tags[stage] = q.tags[stage - 1];
    }
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct Scalar252 {
    limbs: [b16; 16],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct Fold2Side {
    x0: Scalar252,
    p1_low: Scalar252,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct Fold3Side {
    x0: Scalar252,
    p1_low: Scalar252,
    p2_low: Scalar252,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct FoldSideDelayIn {
    value: Fold3Side,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct FoldSideDelayOut {
    value: Fold3Side,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[doc(hidden)]
pub struct FoldSideDelay {
    x0: DFF<[Scalar252; 17]>,
    p1_low: DFF<[Scalar252; 17]>,
    p2_low: DFF<[Scalar252; 17]>,
}

impl Default for FoldSideDelay {
    fn default() -> Self {
        let zero = zero_scalar252();
        Self {
            x0: DFF::new([zero; 17]),
            p1_low: DFF::new([zero; 17]),
            p2_low: DFF::new([zero; 17]),
        }
    }
}

impl SynchronousIO for FoldSideDelay {
    type I = FoldSideDelayIn;
    type O = FoldSideDelayOut;
    type Kernel = fold_side_delay_kernel;
}

/// Seventeen-cycle sideband delay transition.
#[kernel]
#[doc(hidden)]
pub fn fold_side_delay_kernel(
    _cr: ClockReset,
    i: FoldSideDelayIn,
    q: FoldSideDelayQ,
) -> (FoldSideDelayOut, FoldSideDelayD) {
    let mut d = FoldSideDelayD::dont_care();
    d.x0[0] = i.value.x0;
    d.p1_low[0] = i.value.p1_low;
    d.p2_low[0] = i.value.p2_low;
    for stage in 1..17 {
        d.x0[stage] = q.x0[stage - 1];
        d.p1_low[stage] = q.p1_low[stage - 1];
        d.p2_low[stage] = q.p2_low[stage - 1];
    }
    (
        FoldSideDelayOut {
            value: Fold3Side {
                x0: q.x0[16],
                p1_low: q.p1_low[16],
                p2_low: q.p2_low[16],
            },
        },
        d,
    )
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct FinalReduceRegs {
    valid: [bool; 6],
    tags: [b16; 6],
    positive: ScalarWords,
    negative: ScalarWords,
    biased: ScalarWords,
    reduced1: ScalarWords,
    reduced2: ScalarWords,
    low: [b64; 2],
    low_borrow: bool,
    high_operand: [b64; 2],
    original: ScalarWords,
    result: ScalarWords,
}

/// Tagged 512-bit reduction request for [`ScalarReducePipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarReducePipeIn {
    /// Marks this digest valid.
    pub valid: bool,
    /// Wide little-endian input.
    pub wide: Words512,
    /// Caller-defined transaction tag.
    pub tag: b16,
}

/// Tagged canonical scalar from [`ScalarReducePipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarReducePipeOut {
    /// Marks this result valid.
    pub valid: bool,
    /// Canonical scalar strictly below the Ed25519 group order.
    pub scalar: ScalarWords,
    /// Delayed transaction tag.
    pub tag: b16,
}

fn zero_scalar252() -> Scalar252 {
    Scalar252 {
        limbs: [b16(0); 16],
    }
}

fn final_reset() -> FinalReduceRegs {
    let zero = ScalarWords { words: [b64(0); 4] };
    FinalReduceRegs {
        valid: [false; 6],
        tags: [b16(0); 6],
        positive: zero,
        negative: zero,
        biased: zero,
        reduced1: zero,
        reduced2: zero,
        low: [b64(0); 2],
        low_borrow: false,
        high_operand: [b64(0); 2],
        original: zero,
        result: zero,
    }
}

/// Fifty-seven-stage, initiation-interval-one reduction modulo `l`.
///
/// Three fixed-constant folds remove bits above bit 251. Six final registered
/// add/subtract stages canonicalize the result. Every loop count and data path
/// is independent of the input value.
#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
pub struct ScalarReducePipe {
    fold1: ConstQMulPipe,
    fold2: ConstQMulPipe,
    fold3: ConstQMulPipe,
    side1: FoldSideDelay,
    side2: FoldSideDelay,
    side3: FoldSideDelay,
    final_reduce: DFF<FinalReduceRegs>,
}

impl Default for ScalarReducePipe {
    fn default() -> Self {
        Self {
            fold1: ConstQMulPipe::default(),
            fold2: ConstQMulPipe::default(),
            fold3: ConstQMulPipe::default(),
            side1: FoldSideDelay::default(),
            side2: FoldSideDelay::default(),
            side3: FoldSideDelay::default(),
            final_reduce: DFF::new(final_reset()),
        }
    }
}

impl SynchronousIO for ScalarReducePipe {
    type I = ScalarReducePipeIn;
    type O = ScalarReducePipeOut;
    type Kernel = scalar_reduce_pipe_kernel;
}

#[kernel]
fn wide_radix16(value: Words512) -> [b16; 32] {
    let mut limbs = [b16(0); 32];
    for word in 0..8 {
        for part in 0..4 {
            limbs[word * 4 + part] = (value.words[word] >> ((part * 16) as u128)).resize::<16>();
        }
    }
    limbs
}

#[kernel]
fn low252(value: [b16; 32]) -> Scalar252 {
    let mut limbs = [b16(0); 16];
    for index in 0..16 {
        limbs[index] = value[index];
    }
    limbs[15] &= b16(0x0fff);
    Scalar252 { limbs }
}

#[kernel]
fn high252(value: [b16; 32]) -> Radix16X17 {
    let mut limbs = [b16(0); 17];
    for index in 0..16 {
        limbs[index] = (value[15 + index] >> 12) | (value[16 + index] << 4);
    }
    limbs[16] = value[31] >> 12;
    Radix16X17 { limbs }
}

#[kernel]
fn product_low252(value: Radix16X25) -> Scalar252 {
    let mut limbs = [b16(0); 16];
    for index in 0..16 {
        limbs[index] = value.limbs[index];
    }
    limbs[15] &= b16(0x0fff);
    Scalar252 { limbs }
}

#[kernel]
fn product_high252(value: Radix16X25) -> Radix16X17 {
    let mut limbs = [b16(0); 17];
    for index in 0..9 {
        limbs[index] = (value.limbs[15 + index] >> 12) | (value.limbs[16 + index] << 4);
    }
    limbs[9] = value.limbs[24] >> 12;
    Radix16X17 { limbs }
}

#[kernel]
fn scalar252_words(value: Scalar252) -> ScalarWords {
    let mut words = [b64(0); 4];
    for word in 0..4 {
        for part in 0..4 {
            words[word] |= value.limbs[word * 4 + part].resize::<64>() << ((part * 16) as u128);
        }
    }
    ScalarWords { words }
}

#[kernel]
fn add_words(a: ScalarWords, b: ScalarWords) -> ScalarWords {
    let mut output = ScalarWords { words: [b64(0); 4] };
    let mut carry = false;
    for index in 0..4 {
        let first = a.words[index] + b.words[index];
        let first_carry = first < a.words[index];
        let second = first + if carry { b64(1) } else { b64(0) };
        carry = first_carry || second < first;
        output.words[index] = second;
    }
    output
}

#[kernel]
fn two_l() -> ScalarWords {
    ScalarWords {
        words: [
            b64(0xb024c634b9eba7da),
            b64(0x29bdf3bd45ef39ac),
            b64(0),
            b64(0x2000000000000000),
        ],
    }
}

/// State transition for [`ScalarReducePipe`].
#[kernel]
pub fn scalar_reduce_pipe_kernel(
    _cr: ClockReset,
    i: ScalarReducePipeIn,
    q: ScalarReducePipeQ,
) -> (ScalarReducePipeOut, ScalarReducePipeD) {
    let wide = wide_radix16(i.wide);
    let p1_low = product_low252(q.fold1.product);
    let p2_low = product_low252(q.fold2.product);
    let p3_low = product_low252(q.fold3.product);
    let mut final_next = q.final_reduce;
    final_next.valid = [
        q.fold3.valid,
        q.final_reduce.valid[0],
        q.final_reduce.valid[1],
        q.final_reduce.valid[2],
        q.final_reduce.valid[3],
        q.final_reduce.valid[4],
    ];
    final_next.tags = [
        q.fold3.tag,
        q.final_reduce.tags[0],
        q.final_reduce.tags[1],
        q.final_reduce.tags[2],
        q.final_reduce.tags[3],
        q.final_reduce.tags[4],
    ];
    final_next.positive = add_words(
        scalar252_words(q.side3.value.x0),
        scalar252_words(q.side3.value.p2_low),
    );
    final_next.negative = add_words(
        scalar252_words(q.side3.value.p1_low),
        scalar252_words(p3_low),
    );
    final_next.biased = scalar_sub_kernel(
        add_words(q.final_reduce.positive, two_l()),
        q.final_reduce.negative,
    );
    let modulus = scalar_l_kernel(false);
    final_next.reduced1 = if scalar_ge_kernel(q.final_reduce.biased, modulus) {
        scalar_sub_kernel(q.final_reduce.biased, modulus)
    } else {
        q.final_reduce.biased
    };
    final_next.reduced2 = if scalar_ge_kernel(q.final_reduce.reduced1, modulus) {
        scalar_sub_kernel(q.final_reduce.reduced1, modulus)
    } else {
        q.final_reduce.reduced1
    };
    let low0 = q.final_reduce.reduced2.words[0] - modulus.words[0];
    let borrow0 = q.final_reduce.reduced2.words[0] < modulus.words[0];
    let low1_first = q.final_reduce.reduced2.words[1] - modulus.words[1];
    let borrow1_first = q.final_reduce.reduced2.words[1] < modulus.words[1];
    let low1 = low1_first - if borrow0 { b64(1) } else { b64(0) };
    final_next.low = [low0, low1];
    final_next.low_borrow = borrow1_first || (borrow0 && low1_first == b64(0));
    final_next.high_operand = [
        q.final_reduce.reduced2.words[2],
        q.final_reduce.reduced2.words[3],
    ];
    final_next.original = q.final_reduce.reduced2;
    let high0 = q.final_reduce.high_operand[0]
        - modulus.words[2]
        - if q.final_reduce.low_borrow {
            b64(1)
        } else {
            b64(0)
        };
    let borrow_high0 = q.final_reduce.high_operand[0] < modulus.words[2]
        || (q.final_reduce.low_borrow && q.final_reduce.high_operand[0] == modulus.words[2]);
    let high1_first = q.final_reduce.high_operand[1] - modulus.words[3];
    let borrow_high1_first = q.final_reduce.high_operand[1] < modulus.words[3];
    let high1 = high1_first - if borrow_high0 { b64(1) } else { b64(0) };
    let underflow = borrow_high1_first || (borrow_high0 && high1_first == b64(0));
    final_next.result = if underflow {
        q.final_reduce.original
    } else {
        ScalarWords {
            words: [q.final_reduce.low[0], q.final_reduce.low[1], high0, high1],
        }
    };

    let output = ScalarReducePipeOut {
        valid: q.final_reduce.valid[5],
        scalar: q.final_reduce.result,
        tag: q.final_reduce.tags[5],
    };
    let mut d = ScalarReducePipeD::dont_care();
    d.fold1 = ConstQMulIn {
        valid: i.valid,
        x: high252(wide),
        tag: i.tag,
    };
    d.fold2 = ConstQMulIn {
        valid: q.fold1.valid,
        x: product_high252(q.fold1.product),
        tag: q.fold1.tag,
    };
    d.fold3 = ConstQMulIn {
        valid: q.fold2.valid,
        x: product_high252(q.fold2.product),
        tag: q.fold2.tag,
    };
    let zero = Scalar252 {
        limbs: [b16(0); 16],
    };
    d.side1 = FoldSideDelayIn {
        value: Fold3Side {
            x0: low252(wide),
            p1_low: zero,
            p2_low: zero,
        },
    };
    d.side2 = FoldSideDelayIn {
        value: Fold3Side {
            x0: q.side1.value.x0,
            p1_low,
            p2_low: zero,
        },
    };
    d.side3 = FoldSideDelayIn {
        value: Fold3Side {
            x0: q.side2.value.x0,
            p1_low: q.side2.value.p1_low,
            p2_low,
        },
    };
    d.final_reduce = final_next;
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffCellIn {
    x: [b16; 16],
    y: [b16; 16],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffCellOut {
    coefficient: b36,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffCellRegs {
    x: [b16; 16],
    y: [b16; 16],
    products: [b32; 16],
    pairs: [b33; 8],
    quads: [b34; 4],
    octets: [b35; 2],
    coefficient: b36,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[doc(hidden)]
pub struct WideCoeffCell {
    regs: DFF<WideCoeffCellRegs>,
}

impl Default for WideCoeffCell {
    fn default() -> Self {
        Self {
            regs: DFF::new(WideCoeffCellRegs {
                x: [b16(0); 16],
                y: [b16(0); 16],
                products: [b32(0); 16],
                pairs: [b33(0); 8],
                quads: [b34(0); 4],
                octets: [b35(0); 2],
                coefficient: b36(0),
            }),
        }
    }
}

impl SynchronousIO for WideCoeffCell {
    type I = WideCoeffCellIn;
    type O = WideCoeffCellOut;
    type Kernel = wide_coeff_cell_kernel;
}

/// Pipeline transition for one coefficient of a full scalar product.
#[kernel]
#[doc(hidden)]
pub fn wide_coeff_cell_kernel(
    _cr: ClockReset,
    i: WideCoeffCellIn,
    q: WideCoeffCellQ,
) -> (WideCoeffCellOut, WideCoeffCellD) {
    let r = q.regs;
    let mut n = r;
    n.x = i.x;
    n.y = i.y;
    for index in 0..16 {
        n.products[index] = r.x[index].resize::<32>() * r.y[index].resize::<32>();
    }
    for index in 0..8 {
        n.pairs[index] =
            r.products[index * 2].resize::<33>() + r.products[index * 2 + 1].resize::<33>();
    }
    for index in 0..4 {
        n.quads[index] = r.pairs[index * 2].resize::<34>() + r.pairs[index * 2 + 1].resize::<34>();
    }
    n.octets[0] = r.quads[0].resize::<35>() + r.quads[1].resize::<35>();
    n.octets[1] = r.quads[2].resize::<35>() + r.quads[3].resize::<35>();
    n.coefficient = r.octets[0].resize::<36>() + r.octets[1].resize::<36>();
    (
        WideCoeffCellOut {
            coefficient: r.coefficient,
        },
        WideCoeffCellD { regs: n },
    )
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffBankIn {
    windows: [[b16; 16]; 4],
    y: [b16; 16],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffBankOut {
    coefficients: [b36; 4],
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
#[doc(hidden)]
pub struct WideCoeffBank {
    cell0: WideCoeffCell,
    cell1: WideCoeffCell,
    cell2: WideCoeffCell,
    cell3: WideCoeffCell,
}

impl SynchronousIO for WideCoeffBank {
    type I = WideCoeffBankIn;
    type O = WideCoeffBankOut;
    type Kernel = wide_coeff_bank_kernel;
}

/// Wiring transition for four full-product coefficients.
#[kernel]
#[doc(hidden)]
pub fn wide_coeff_bank_kernel(
    _cr: ClockReset,
    i: WideCoeffBankIn,
    q: WideCoeffBankQ,
) -> (WideCoeffBankOut, WideCoeffBankD) {
    let output = WideCoeffBankOut {
        coefficients: [
            q.cell0.coefficient,
            q.cell1.coefficient,
            q.cell2.coefficient,
            q.cell3.coefficient,
        ],
    };
    let mut d = WideCoeffBankD::dont_care();
    d.cell0 = WideCoeffCellIn {
        x: i.windows[0],
        y: i.y,
    };
    d.cell1 = WideCoeffCellIn {
        x: i.windows[1],
        y: i.y,
    };
    d.cell2 = WideCoeffCellIn {
        x: i.windows[2],
        y: i.y,
    };
    d.cell3 = WideCoeffCellIn {
        x: i.windows[3],
        y: i.y,
    };
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffBankPairIn {
    windows: [[b16; 16]; 8],
    y: [b16; 16],
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideCoeffBankPairOut {
    coefficients: [b36; 8],
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
#[doc(hidden)]
pub struct WideCoeffBankPair {
    bank0: WideCoeffBank,
    bank1: WideCoeffBank,
}

impl SynchronousIO for WideCoeffBankPair {
    type I = WideCoeffBankPairIn;
    type O = WideCoeffBankPairOut;
    type Kernel = wide_coeff_bank_pair_kernel;
}

/// Wiring transition for eight full-product coefficients.
#[kernel]
#[doc(hidden)]
pub fn wide_coeff_bank_pair_kernel(
    _cr: ClockReset,
    i: WideCoeffBankPairIn,
    q: WideCoeffBankPairQ,
) -> (WideCoeffBankPairOut, WideCoeffBankPairD) {
    let output = WideCoeffBankPairOut {
        coefficients: [
            q.bank0.coefficients[0],
            q.bank0.coefficients[1],
            q.bank0.coefficients[2],
            q.bank0.coefficients[3],
            q.bank1.coefficients[0],
            q.bank1.coefficients[1],
            q.bank1.coefficients[2],
            q.bank1.coefficients[3],
        ],
    };
    let mut d = WideCoeffBankPairD::dont_care();
    d.bank0 = WideCoeffBankIn {
        windows: [i.windows[0], i.windows[1], i.windows[2], i.windows[3]],
        y: i.y,
    };
    d.bank1 = WideCoeffBankIn {
        windows: [i.windows[4], i.windows[5], i.windows[6], i.windows[7]],
        y: i.y,
    };
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideNormalizeState {
    coefficients: [b36; 32],
    limbs: [b16; 32],
    carry: b21,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideNormalizePairIn {
    state: WideNormalizeState,
    start: b7,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideNormalizePairOut {
    state: WideNormalizeState,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
#[doc(hidden)]
pub struct WideNormalizePair {
    state: DFF<WideNormalizeState>,
}

impl Default for WideNormalizePair {
    fn default() -> Self {
        Self {
            state: DFF::new(WideNormalizeState {
                coefficients: [b36(0); 32],
                limbs: [b16(0); 32],
                carry: b21(0),
            }),
        }
    }
}

impl SynchronousIO for WideNormalizePair {
    type I = WideNormalizePairIn;
    type O = WideNormalizePairOut;
    type Kernel = wide_normalize_pair_kernel;
}

/// Two-limb full-product normalization transition.
#[kernel]
#[doc(hidden)]
pub fn wide_normalize_pair_kernel(
    _cr: ClockReset,
    i: WideNormalizePairIn,
    q: WideNormalizePairQ,
) -> (WideNormalizePairOut, WideNormalizePairD) {
    let mut next = i.state;
    let first = i.state.coefficients[i.start].resize::<37>() + i.state.carry.resize::<37>();
    let second_index = i.start + b7(1);
    let second = i.state.coefficients[second_index].resize::<37>() + (first >> 16);
    next.limbs[i.start] = first.resize::<16>();
    next.limbs[second_index] = second.resize::<16>();
    next.carry = (second >> 16).resize::<21>();
    (
        WideNormalizePairOut { state: q.state },
        WideNormalizePairD { state: next },
    )
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideNormalizeGroupIn {
    state: WideNormalizeState,
    start: b7,
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct WideNormalizeGroupOut {
    state: WideNormalizeState,
}

#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
#[doc(hidden)]
pub struct WideNormalizeGroup {
    pair0: WideNormalizePair,
    pair1: WideNormalizePair,
    pair2: WideNormalizePair,
    pair3: WideNormalizePair,
}

impl SynchronousIO for WideNormalizeGroup {
    type I = WideNormalizeGroupIn;
    type O = WideNormalizeGroupOut;
    type Kernel = wide_normalize_group_kernel;
}

/// Wiring transition for four full-product normalization steps.
#[kernel]
#[doc(hidden)]
pub fn wide_normalize_group_kernel(
    _cr: ClockReset,
    i: WideNormalizeGroupIn,
    q: WideNormalizeGroupQ,
) -> (WideNormalizeGroupOut, WideNormalizeGroupD) {
    let mut d = WideNormalizeGroupD::dont_care();
    d.pair0 = WideNormalizePairIn {
        state: i.state,
        start: i.start,
    };
    d.pair1 = WideNormalizePairIn {
        state: q.pair0.state,
        start: i.start + b7(2),
    };
    d.pair2 = WideNormalizePairIn {
        state: q.pair1.state,
        start: i.start + b7(4),
    };
    d.pair3 = WideNormalizePairIn {
        state: q.pair2.state,
        start: i.start + b7(6),
    };
    (
        WideNormalizeGroupOut {
            state: q.pair3.state,
        },
        d,
    )
}

/// Exact 528-bit output of [`ScalarWideMulPipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct Radix16X33 {
    /// Product limbs, least significant first.
    pub limbs: [b16; 33],
}

/// Tagged request for [`ScalarWideMulPipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarWideMulIn {
    /// Marks this request valid.
    pub valid: bool,
    /// First factor.
    pub x: ScalarWords,
    /// Second factor.
    pub y: ScalarWords,
    /// Sideband addend delayed with the product.
    pub addend: ScalarWords,
    /// Caller tag.
    pub tag: b16,
}

/// Tagged response from [`ScalarWideMulPipe`].
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarWideMulOut {
    /// Marks this response valid.
    pub valid: bool,
    /// Exact product.
    pub product: Radix16X33,
    /// Delayed sideband addend.
    pub addend: ScalarWords,
    /// Delayed caller tag.
    pub tag: b16,
}

/// Twenty-two-stage, 512-DSP full scalar multiplier.
#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
pub struct ScalarWideMulPipe {
    banks0: WideCoeffBankPair,
    banks1: WideCoeffBankPair,
    banks2: WideCoeffBankPair,
    banks3: WideCoeffBankPair,
    normalize0: WideNormalizeGroup,
    normalize1: WideNormalizeGroup,
    normalize2: WideNormalizeGroup,
    normalize3: WideNormalizeGroup,
    valid: DFF<[bool; 22]>,
    tags: DFF<[b16; 22]>,
    addends: DFF<[ScalarWords; 22]>,
}

impl Default for ScalarWideMulPipe {
    fn default() -> Self {
        let zero = ScalarWords { words: [b64(0); 4] };
        Self {
            banks0: WideCoeffBankPair::default(),
            banks1: WideCoeffBankPair::default(),
            banks2: WideCoeffBankPair::default(),
            banks3: WideCoeffBankPair::default(),
            normalize0: WideNormalizeGroup::default(),
            normalize1: WideNormalizeGroup::default(),
            normalize2: WideNormalizeGroup::default(),
            normalize3: WideNormalizeGroup::default(),
            valid: DFF::new([false; 22]),
            tags: DFF::new([b16(0); 22]),
            addends: DFF::new([zero; 22]),
        }
    }
}

impl SynchronousIO for ScalarWideMulPipe {
    type I = ScalarWideMulIn;
    type O = ScalarWideMulOut;
    type Kernel = scalar_wide_mul_kernel;
}

#[kernel]
fn scalar_radix16(value: ScalarWords) -> [b16; 16] {
    let mut limbs = [b16(0); 16];
    for word in 0..4 {
        for part in 0..4 {
            limbs[word * 4 + part] = (value.words[word] >> ((part * 16) as u128)).resize::<16>();
        }
    }
    limbs
}

/// State transition for [`ScalarWideMulPipe`].
#[kernel]
pub fn scalar_wide_mul_kernel(
    _cr: ClockReset,
    i: ScalarWideMulIn,
    q: ScalarWideMulPipeQ,
) -> (ScalarWideMulOut, ScalarWideMulPipeD) {
    let x = scalar_radix16(i.x);
    let y = scalar_radix16(i.y);
    let bank_coefficients = [
        q.banks0.coefficients,
        q.banks1.coefficients,
        q.banks2.coefficients,
        q.banks3.coefficients,
    ];
    let mut coefficients = [b36(0); 32];
    let mut windows = [[[b16(0); 16]; 8]; 4];
    for bank_pair in 0..4 {
        for slot in 0..8 {
            let coefficient = bank_pair * 8 + slot;
            coefficients[coefficient] = bank_coefficients[bank_pair][slot];
            for x_limb in 0..16 {
                for y_limb in 0..16 {
                    if coefficient == x_limb + y_limb {
                        windows[bank_pair][slot][y_limb] = x[x_limb];
                    }
                }
            }
        }
    }
    let initial = WideNormalizeState {
        coefficients,
        limbs: [b16(0); 32],
        carry: b21(0),
    };
    let normalized = q.normalize3.state;
    let product = Radix16X33 {
        limbs: [
            normalized.limbs[0],
            normalized.limbs[1],
            normalized.limbs[2],
            normalized.limbs[3],
            normalized.limbs[4],
            normalized.limbs[5],
            normalized.limbs[6],
            normalized.limbs[7],
            normalized.limbs[8],
            normalized.limbs[9],
            normalized.limbs[10],
            normalized.limbs[11],
            normalized.limbs[12],
            normalized.limbs[13],
            normalized.limbs[14],
            normalized.limbs[15],
            normalized.limbs[16],
            normalized.limbs[17],
            normalized.limbs[18],
            normalized.limbs[19],
            normalized.limbs[20],
            normalized.limbs[21],
            normalized.limbs[22],
            normalized.limbs[23],
            normalized.limbs[24],
            normalized.limbs[25],
            normalized.limbs[26],
            normalized.limbs[27],
            normalized.limbs[28],
            normalized.limbs[29],
            normalized.limbs[30],
            normalized.limbs[31],
            normalized.carry.resize::<16>(),
        ],
    };
    let output = ScalarWideMulOut {
        valid: q.valid[21],
        product,
        addend: q.addends[21],
        tag: q.tags[21],
    };
    let mut d = ScalarWideMulPipeD::dont_care();
    d.banks0 = WideCoeffBankPairIn {
        windows: windows[0],
        y,
    };
    d.banks1 = WideCoeffBankPairIn {
        windows: windows[1],
        y,
    };
    d.banks2 = WideCoeffBankPairIn {
        windows: windows[2],
        y,
    };
    d.banks3 = WideCoeffBankPairIn {
        windows: windows[3],
        y,
    };
    d.normalize0 = WideNormalizeGroupIn {
        state: initial,
        start: b7(0),
    };
    d.normalize1 = WideNormalizeGroupIn {
        state: q.normalize0.state,
        start: b7(8),
    };
    d.normalize2 = WideNormalizeGroupIn {
        state: q.normalize1.state,
        start: b7(16),
    };
    d.normalize3 = WideNormalizeGroupIn {
        state: q.normalize2.state,
        start: b7(24),
    };
    d.valid[0] = i.valid;
    d.tags[0] = i.tag;
    d.addends[0] = i.addend;
    for stage in 1..22 {
        d.valid[stage] = q.valid[stage - 1];
        d.tags[stage] = q.tags[stage - 1];
        d.addends[stage] = q.addends[stage - 1];
    }
    (output, d)
}

#[derive(Copy, Clone, Debug, PartialEq, Digital)]
#[doc(hidden)]
pub struct MulAddStage {
    valid: bool,
    wide: Words512,
    tag: b16,
}

/// Tagged scalar multiply-add request.
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarMulAddIn {
    /// Marks this operation valid.
    pub valid: bool,
    /// First multiplicand.
    pub a: ScalarWords,
    /// Second multiplicand.
    pub b: ScalarWords,
    /// Addend.
    pub c: ScalarWords,
    /// Caller tag.
    pub tag: b16,
}

/// Tagged canonical scalar multiply-add response.
#[derive(Copy, Clone, Debug, PartialEq, Digital)]
pub struct ScalarMulAddOut {
    /// Marks the reduced scalar valid.
    pub valid: bool,
    /// Canonical scalar result.
    pub scalar: ScalarWords,
    /// Delayed caller tag.
    pub tag: b16,
}

/// Fully pipelined `a * b + c mod l` engine.
#[derive(Clone, Debug, Synchronous, SynchronousDQ)]
pub struct ScalarMulAddPipe {
    multiplier: ScalarWideMulPipe,
    add: DFF<MulAddStage>,
    reducer: ScalarReducePipe,
}

#[cfg(test)]
#[derive(Clone, Debug, Synchronous, SynchronousDQ, Default)]
struct DualReducerHarness {
    first: ScalarReducePipe,
    second: ScalarReducePipe,
}

#[cfg(test)]
impl SynchronousIO for DualReducerHarness {
    type I = ScalarReducePipeIn;
    type O = ScalarReducePipeOut;
    type Kernel = dual_reducer_harness_kernel;
}

#[cfg(test)]
#[kernel]
fn dual_reducer_harness_kernel(
    _cr: ClockReset,
    i: ScalarReducePipeIn,
    q: DualReducerHarnessQ,
) -> (ScalarReducePipeOut, DualReducerHarnessD) {
    let mut d = DualReducerHarnessD::dont_care();
    d.first = ScalarReducePipeIn {
        valid: false,
        wide: Words512 { words: [b64(0); 8] },
        tag: b16(0),
    };
    d.second = i;
    (q.second, d)
}

impl Default for ScalarMulAddPipe {
    fn default() -> Self {
        Self {
            multiplier: ScalarWideMulPipe::default(),
            add: DFF::new(MulAddStage {
                valid: false,
                wide: Words512::zero(),
                tag: b16(0),
            }),
            reducer: ScalarReducePipe::default(),
        }
    }
}

impl SynchronousIO for ScalarMulAddPipe {
    type I = ScalarMulAddIn;
    type O = ScalarMulAddOut;
    type Kernel = scalar_muladd_kernel;
}

/// Add a canonical scalar addend to a wide radix-`2^16` product.
#[kernel]
pub fn scalar_product_addend_kernel(product: Radix16X33, addend: ScalarWords) -> Words512 {
    let addend_low = scalar_radix16(addend);
    let mut addend = [b16(0); 32];
    for index in 0..16 {
        addend[index] = addend_low[index];
    }
    let mut limbs = [b16(0); 32];
    let mut carry = b17(0);
    for index in 0..32 {
        let sum = product.limbs[index].resize::<17>() + addend[index].resize::<17>() + carry;
        limbs[index] = sum.resize::<16>();
        carry = sum >> 16;
    }
    let mut words = [b64(0); 8];
    for word in 0..8 {
        for part in 0..4 {
            words[word] |= limbs[word * 4 + part].resize::<64>() << ((part * 16) as u128);
        }
    }
    Words512 { words }
}

/// State transition for [`ScalarMulAddPipe`].
#[kernel]
pub fn scalar_muladd_kernel(
    _cr: ClockReset,
    i: ScalarMulAddIn,
    q: ScalarMulAddPipeQ,
) -> (ScalarMulAddOut, ScalarMulAddPipeD) {
    let mut d = ScalarMulAddPipeD::dont_care();
    d.multiplier = ScalarWideMulIn {
        valid: i.valid,
        x: i.a,
        y: i.b,
        addend: i.c,
        tag: i.tag,
    };
    d.add = MulAddStage {
        valid: q.multiplier.valid,
        wide: scalar_product_addend_kernel(q.multiplier.product, q.multiplier.addend),
        tag: q.multiplier.tag,
    };
    d.reducer = ScalarReducePipeIn {
        valid: q.add.valid,
        wide: q.add.wide,
        tag: q.add.tag,
    };
    (
        ScalarMulAddOut {
            valid: q.reducer.valid,
            scalar: q.reducer.scalar,
            tag: q.reducer.tag,
        },
        d,
    )
}

#[cfg(test)]
mod tests {
    use super::*;
    use curve25519_dalek::scalar::Scalar;
    use num_bigint::BigUint;
    use std::iter;

    fn to_x17(value: &BigUint) -> Radix16X17 {
        let mask = BigUint::from(0xffffu32);
        let mut limbs = [b16(0); 17];
        for (index, limb) in limbs.iter_mut().enumerate() {
            let digits = ((value >> (index * 16)) & &mask).to_u32_digits();
            *limb = b16(digits.first().copied().unwrap_or(0) as u128);
        }
        Radix16X17 { limbs }
    }

    fn from_x25(value: Radix16X25) -> BigUint {
        value
            .limbs
            .iter()
            .enumerate()
            .fold(BigUint::from(0u8), |sum, (index, limb)| {
                sum + (BigUint::from(limb.raw()) << (index * 16))
            })
    }

    fn digest_words(bytes: [u8; 64]) -> Words512 {
        let mut value = Words512::zero();
        for (index, chunk) in bytes.chunks_exact(8).enumerate() {
            value.words[index] = b64(u64::from_le_bytes(chunk.try_into().unwrap()) as u128);
        }
        value
    }

    fn scalar_bytes(value: ScalarWords) -> [u8; 32] {
        let mut bytes = [0u8; 32];
        for (index, word) in value.words.iter().enumerate() {
            bytes[index * 8..index * 8 + 8].copy_from_slice(&(word.raw() as u64).to_le_bytes());
        }
        bytes
    }

    fn scalar_words(bytes: [u8; 32]) -> ScalarWords {
        let mut words = [b64(0); 4];
        for (index, chunk) in bytes.chunks_exact(8).enumerate() {
            words[index] = b64(u64::from_le_bytes(chunk.try_into().unwrap()) as u128);
        }
        ScalarWords { words }
    }

    fn from_x33(value: Radix16X33) -> BigUint {
        value
            .limbs
            .iter()
            .enumerate()
            .fold(BigUint::from(0u8), |sum, (index, limb)| {
                sum + (BigUint::from(limb.raw()) << (index * 16))
            })
    }

    #[test]
    fn constant_multiplier_accepts_one_value_per_cycle() {
        let q = BigUint::parse_bytes(b"14def9dea2f79cd65812631a5cf5d3ed", 16).unwrap();
        let values: Vec<_> = (0u64..40)
            .map(|index| {
                (BigUint::from(0x9e3779b97f4a7c15u64) * BigUint::from(index + 1)
                    + (BigUint::from(index * 13 + 5) << 244))
                    % (BigUint::from(1u8) << 272)
            })
            .collect();
        let idle = ConstQMulIn {
            valid: false,
            x: Radix16X17 {
                limbs: [b16(0); 17],
            },
            tag: b16(0),
        };
        let inputs = values.iter().enumerate().map(|(index, value)| ConstQMulIn {
            valid: true,
            x: to_x17(value),
            tag: b16(index as u128),
        });
        let outputs: Vec<_> = ConstQMulPipe::default()
            .run(
                iter::once(idle)
                    .chain(inputs)
                    .chain(iter::repeat_n(idle, 24))
                    .with_reset(1)
                    .clock_pos_edge(100),
            )
            .synchronous_sample()
            .map(|sample| sample.output)
            .filter(|output| output.valid)
            .collect();
        assert_eq!(outputs.len(), values.len());
        for (index, (output, input)) in outputs.into_iter().zip(values).enumerate() {
            assert_eq!(output.tag, b16(index as u128));
            assert_eq!(from_x25(output.product), input * &q);
        }
    }

    #[test]
    fn reducer_matches_dalek_for_back_to_back_digests() {
        std::thread::Builder::new()
            .name("scalar-reducer-sim".into())
            .stack_size(64 * 1024 * 1024)
            .spawn(reducer_matches_dalek_body)
            .unwrap()
            .join()
            .unwrap();
    }

    fn reducer_matches_dalek_body() {
        let digests: Vec<_> = (0u8..32)
            .map(|case| {
                let mut bytes = [0u8; 64];
                for (index, byte) in bytes.iter_mut().enumerate() {
                    *byte = case
                        .wrapping_mul(37)
                        .wrapping_add(index as u8)
                        .rotate_left((index % 7) as u32);
                }
                bytes
            })
            .collect();
        let idle = ScalarReducePipeIn {
            valid: false,
            wide: Words512::zero(),
            tag: b16(0),
        };
        let inputs = digests
            .iter()
            .enumerate()
            .map(|(index, digest)| ScalarReducePipeIn {
                valid: true,
                wide: digest_words(*digest),
                tag: b16(index as u128),
            });
        let outputs: Vec<_> = ScalarReducePipe::default()
            .run(
                iter::once(idle)
                    .chain(inputs)
                    .chain(iter::repeat_n(idle, 70))
                    .with_reset(1)
                    .clock_pos_edge(100),
            )
            .synchronous_sample()
            .map(|sample| sample.output)
            .filter(|output| output.valid)
            .collect();
        assert_eq!(outputs.len(), digests.len());
        for (index, (output, digest)) in outputs.into_iter().zip(digests).enumerate() {
            assert_eq!(output.tag, b16(index as u128));
            assert_eq!(
                scalar_bytes(output.scalar),
                Scalar::from_bytes_mod_order_wide(&digest).to_bytes()
            );
        }
    }

    #[test]
    fn wide_multiplier_and_muladd_match_dalek() {
        std::thread::Builder::new()
            .name("scalar-muladd-sim".into())
            .stack_size(128 * 1024 * 1024)
            .spawn(wide_multiplier_and_muladd_body)
            .unwrap()
            .join()
            .unwrap();
    }

    fn wide_multiplier_and_muladd_body() {
        let cases: Vec<_> = (0u8..24)
            .map(|case| {
                let mut a_bytes = [0u8; 64];
                let mut b_bytes = [0u8; 64];
                let mut c_bytes = [0u8; 64];
                for index in 0..64 {
                    a_bytes[index] = case.wrapping_mul(17).wrapping_add(index as u8);
                    b_bytes[index] = case.wrapping_mul(29).wrapping_add((index * 3) as u8);
                    c_bytes[index] = case.wrapping_mul(43).wrapping_add((index * 5) as u8);
                }
                (
                    Scalar::from_bytes_mod_order_wide(&a_bytes),
                    Scalar::from_bytes_mod_order_wide(&b_bytes),
                    Scalar::from_bytes_mod_order_wide(&c_bytes),
                )
            })
            .collect();
        let zero = ScalarWords { words: [b64(0); 4] };
        let wide_idle = ScalarWideMulIn {
            valid: false,
            x: zero,
            y: zero,
            addend: zero,
            tag: b16(0),
        };
        let wide_inputs = cases
            .iter()
            .enumerate()
            .map(|(index, (a, b, c))| ScalarWideMulIn {
                valid: true,
                x: scalar_words(a.to_bytes()),
                y: scalar_words(b.to_bytes()),
                addend: scalar_words(c.to_bytes()),
                tag: b16(index as u128),
            });
        let wide_outputs: Vec<_> = ScalarWideMulPipe::default()
            .run(
                iter::once(wide_idle)
                    .chain(wide_inputs)
                    .chain(iter::repeat_n(wide_idle, 30))
                    .with_reset(1)
                    .clock_pos_edge(100),
            )
            .synchronous_sample()
            .map(|sample| sample.output)
            .filter(|output| output.valid)
            .collect();
        assert_eq!(wide_outputs.len(), cases.len());
        for (index, (output, (a, b, c))) in wide_outputs.into_iter().zip(cases.iter()).enumerate() {
            assert_eq!(output.tag, b16(index as u128));
            assert_eq!(output.addend, scalar_words(c.to_bytes()));
            assert_eq!(
                from_x33(output.product),
                BigUint::from_bytes_le(&a.to_bytes()) * BigUint::from_bytes_le(&b.to_bytes())
            );
        }

        let muladd_idle = ScalarMulAddIn {
            valid: false,
            a: zero,
            b: zero,
            c: zero,
            tag: b16(0),
        };
        let muladd_inputs = cases
            .iter()
            .enumerate()
            .map(|(index, (a, b, c))| ScalarMulAddIn {
                valid: true,
                a: scalar_words(a.to_bytes()),
                b: scalar_words(b.to_bytes()),
                c: scalar_words(c.to_bytes()),
                tag: b16(index as u128),
            });
        let outputs: Vec<_> = ScalarMulAddPipe::default()
            .run(
                iter::once(muladd_idle)
                    .chain(muladd_inputs)
                    .chain(iter::repeat_n(muladd_idle, 100))
                    .with_reset(1)
                    .clock_pos_edge(100),
            )
            .synchronous_sample()
            .map(|sample| sample.output)
            .filter(|output| output.valid)
            .collect();
        assert_eq!(outputs.len(), cases.len());
        for (index, (output, (a, b, c))) in outputs.into_iter().zip(cases).enumerate() {
            assert_eq!(output.tag, b16(index as u128));
            assert_eq!(scalar_bytes(output.scalar), (a * b + c).to_bytes());
        }
    }

    #[test]
    fn muladd_preserves_a_single_sparse_valid() {
        std::thread::Builder::new()
            .name("scalar-muladd-sparse".into())
            .stack_size(128 * 1024 * 1024)
            .spawn(muladd_sparse_body)
            .unwrap()
            .join()
            .unwrap();
    }

    fn muladd_sparse_body() {
        let zero = ScalarWords { words: [b64(0); 4] };
        let idle = ScalarMulAddIn {
            valid: false,
            a: zero,
            b: zero,
            c: zero,
            tag: b16(0),
        };
        let request = ScalarMulAddIn {
            valid: true,
            a: ScalarWords {
                words: [b64(3), b64(0), b64(0), b64(0)],
            },
            b: ScalarWords {
                words: [b64(5), b64(0), b64(0), b64(0)],
            },
            c: ScalarWords {
                words: [b64(7), b64(0), b64(0), b64(0)],
            },
            tag: b16(0x5151),
        };
        let outputs: Vec<_> = ScalarMulAddPipe::default()
            .run(
                iter::once(idle)
                    .chain(iter::once(request))
                    .chain(iter::repeat_n(idle, 120))
                    .with_reset(1)
                    .clock_pos_edge(100),
            )
            .synchronous_sample()
            .map(|sample| sample.output)
            .filter(|output| output.valid)
            .collect();
        assert_eq!(outputs.len(), 1);
        assert_eq!(outputs[0].tag, b16(0x5151));
        assert_eq!(outputs[0].scalar.words[0], b64(22));
    }

    #[test]
    fn second_reducer_instance_preserves_sparse_valid() {
        std::thread::Builder::new()
            .name("dual-reducer-sparse".into())
            .stack_size(128 * 1024 * 1024)
            .spawn(|| {
                let idle = ScalarReducePipeIn {
                    valid: false,
                    wide: Words512 { words: [b64(0); 8] },
                    tag: b16(0),
                };
                let request = ScalarReducePipeIn {
                    valid: true,
                    wide: Words512 {
                        words: [
                            b64(22),
                            b64(0),
                            b64(0),
                            b64(0),
                            b64(0),
                            b64(0),
                            b64(0),
                            b64(0),
                        ],
                    },
                    tag: b16(0x5252),
                };
                let outputs: Vec<_> = DualReducerHarness::default()
                    .run(
                        iter::once(idle)
                            .chain(iter::once(request))
                            .chain(iter::repeat_n(idle, 100))
                            .with_reset(1)
                            .clock_pos_edge(100),
                    )
                    .synchronous_sample()
                    .map(|sample| sample.output)
                    .filter(|output| output.valid)
                    .collect();
                assert_eq!(outputs.len(), 1);
                assert_eq!(outputs[0].tag, b16(0x5252));
                assert_eq!(outputs[0].scalar.words[0], b64(22));
            })
            .unwrap()
            .join()
            .unwrap();
    }

    #[test]
    fn compiles_to_rhdl_netlist() -> miette::Result<()> {
        ConstQMulPipe::default()
            .descriptor(ScopedName::top())?
            .hdl()?;
        ScalarReducePipe::default()
            .descriptor(ScopedName::top())?
            .hdl()?;
        Ok(())
    }
}