rhdl_ed25519_scalable_core/

programs.rs

//! The elliptic-curve programs the field machines run.
//!
//! | Program | Bundles | Multiplies | Use |
//! |---|---:|---:|---|
//! | [`fixed_base`] | 392 | 378 | `[r]B` for signing and key expansion |
//! | [`fixed_base_affine`] | 657 | 645 | `[r]B` plus its own inversion (fused mode) |
//! | [`batch_compress`]`(32)` | 423 | 422 | affine `x, y` of 32 points with one inversion |
//!
//! Each is plain Rust against [`Builder`],
//! compiled to bundles, and checked against `curve25519-dalek` through the
//! interpreter in this module's tests. The verification and key-loading
//! programs ([`verify_known_key`], [`decode_candidate`], [`table_bases`],
//! [`digit_table`]) are built the same way.
//!
//! # The fixed-base walk
//!
//! `[r]B = Σ dᵢ · (2^offsetᵢ · B)`: one table per digit, so no doublings. Each
//! digit costs one 7-multiply mixed addition, and each lookup costs one table
//! scan step per entry pair. Choosing the digit sizes balances the two:
//!
#![doc = concat!("<div style=\"overflow-x:auto\">", include_str!("../docs/images/slot-balance.svg"), "</div>")]
//!
//! The chosen recoding ([`host::recode`](crate::host::recode)) has 39 signed
//! radix-32 digits and 15 signed radix-16 digits (255 bits; every scalar here
//! is reduced mod `l` first, so the top digit stays small). How the program
//! went from 492 to 392 bundles:
//!
//! | Change | Bundles |
//! |---|---:|
//! | Dalek's radix-16 walk: 64 additions, 4 doublings, two interleaved chains | 492 |
//! | one table per digit (no doublings); chains start from their first lookup (3 multiplies, not 7) | 458 |
//! | mixed radix 32/16, digits in order | 435 |
//! | radix-16 lookups interleaved among radix-32 ones | 399 |
//! | lookup order from a search with the real scheduler (`LOOKUP_ORDER`) | **392** |
//!
//! Additions commute, so the order is free. In digit order the long scans
//! come first and starve the multiplier; interleaving keeps both busy. The
//! result fills the multiply, add and scan slots 96 %, 97 % and 95 %.
//!
//! # Batch compression
//!
//! Encoding needs `1/Z` (265 multiplies by Fermat's little theorem).
//! Montgomery's trick inverts the product of 32 `Z`s once and recovers each
//! inverse with three more multiplies: 13.2 per point instead of 265.
//!
//! # Key loading
//!
//! A verification key `A` arrives as a 32-byte encoding, and verification
//! needs the fixed-base table of `-A`. Three programs build it, none of them
//! branching: [`decode_candidate`] computes both
//! square-root candidates (the controller picks one), [`table_bases`] walks
//! the 54 bases `2^offset(i) · (-A)` by doubling, and [`digit_table`] turns
//! each base into its 8 or 16 affine Niels entries with one shared inversion.

use super::host::Fp;
use super::host::DIGITS;
use super::microcode::{Builder, Program, V};

/// An extended twisted-Edwards point `(X:Y:Z:T)`, `T = XY/Z`, as virtual
/// values.
#[derive(Clone, Copy)]
struct Point {
    /// `X`.
    x: V,
    /// `Y`.
    y: V,
    /// `Z`.
    z: V,
    /// `T`; some helpers leave a placeholder here when the caller never
    /// reads it.
    t: V,
}

/// Extended point plus affine Niels point from bank L (Dalek `+ AffineNielsPoint`).
///
/// 7 multiplies: 3 against the table entry, then 4 to turn the completed
/// point back into extended coordinates. `l` is `(y+x, y-x, 2dxy)` of the
/// table entry, whose `Z` is 1.
fn madd(b: &mut Builder, p: &Point, l: [V; 3]) -> Point {
    let ypx = b.add(p.y, p.x);
    let ymx = b.sub(p.y, p.x);
    // (Y+X)(y+x), (Y-X)(y-x), and T * 2dxy.
    let pp = b.mul(ypx, l[0]);
    let mm = b.mul(ymx, l[1]);
    let txy = b.mul(p.t, l[2]);
    // 2Z (the entry's Z is 1).
    let z2 = b.add(p.z, p.z);
    // Completed point: affine x = cx / cz, y = cy / ct.
    let cx = b.sub(pp, mm);
    let cy = b.add(pp, mm);
    let cz = b.add(z2, txy);
    let ct = b.sub(z2, txy);
    // Completed to extended.
    Point {
        x: b.mul(cx, ct),
        y: b.mul(cy, cz),
        z: b.mul(cz, ct),
        t: b.mul(cx, cy),
    }
}

/// Projective doubling (a = -1); `T` of the result only when requested.
///
/// Dalek's `ProjectivePoint::double` followed by the completed-to-extended
/// conversion. Doubling does not read `T`, so a chain of doublings computes
/// `T` only on its last step; when `want_t` is false the input's stale `T`
/// is passed through.
fn double(b: &mut Builder, p: &Point, want_t: bool) -> Point {
    let xx = b.sqr(p.x);
    let yy = b.sqr(p.y);
    let zz = b.sqr(p.z);
    // (X + Y)^2.
    let xpy = b.add(p.x, p.y);
    let xpy2 = b.sqr(xpy);
    let zz2 = b.add(zz, zz);
    let yypxx = b.add(yy, xx);
    let yymxx = b.sub(yy, xx);
    // Completed point: X = (X+Y)^2 - (Y^2+X^2) = 2XY, Y = Y^2+X^2,
    // Z = Y^2-X^2, T = 2Z^2 - (Y^2-X^2).
    let cx = b.sub(xpy2, yypxx);
    let ct = b.sub(zz2, yymxx);
    Point {
        x: b.mul(cx, ct),
        y: b.mul(yypxx, yymxx),
        z: b.mul(yymxx, ct),
        t: if want_t { b.mul(cx, yypxx) } else { p.t },
    }
}

/// Affine Niels point from bank L as an extended point with `Z = 1`
/// (bank L may only feed a multiplier's second factor).
///
/// `x = ((y+x) - (y-x)) / 2`, `y = ((y+x) + (y-x)) / 2` and
/// `t = xy = 2dxy / 2d`. Starting a chain this way costs 3 multiplies
/// instead of a 7-multiply addition to the identity. `half` and `inv_2d`
/// are the constants `1/2` and `1/(2d)`.
fn niels_to_ext(b: &mut Builder, l: [V; 3], half: V, inv_2d: V) -> Point {
    // The halving happens on the multiplier because L cannot feed the adder.
    let ypx_2 = b.mul(half, l[0]);
    let ymx_2 = b.mul(half, l[1]);
    let t = b.mul(inv_2d, l[2]);
    let one = b.constant(Fp::ONE);
    Point {
        x: b.sub(ypx_2, ymx_2),
        y: b.add(ypx_2, ymx_2),
        z: one,
        t,
    }
}

/// Full extended addition `P + Q` (Dalek `+ ProjectiveNielsPoint`); `T` of
/// the sum only when requested.
///
/// `d2` is the constant `2d`. When `want_t` is false the result's `t` holds
/// `Z1 * Z2`, a placeholder that saves a multiply where only `X, Y, Z` are
/// output.
fn add_ext(b: &mut Builder, p: &Point, q: &Point, d2: V, want_t: bool) -> Point {
    // Q in projective Niels form: (Y+X, Y-X, Z, 2dT).
    let q_ypx = b.add(q.y, q.x);
    let q_ymx = b.sub(q.y, q.x);
    let q_t2d = b.mul(q.t, d2);
    let ypx = b.add(p.y, p.x);
    let ymx = b.sub(p.y, p.x);
    let pp = b.mul(ypx, q_ypx);
    let mm = b.mul(ymx, q_ymx);
    let tt = b.mul(p.t, q_t2d);
    let zz = b.mul(p.z, q.z);
    let zz2 = b.add(zz, zz);
    // Completed point, then completed to extended as in `madd`.
    let cx = b.sub(pp, mm);
    let cy = b.add(pp, mm);
    let cz = b.add(zz2, tt);
    let ct = b.sub(zz2, tt);
    Point {
        x: b.mul(cx, ct),
        y: b.mul(cy, cz),
        z: b.mul(cz, ct),
        t: if want_t { b.mul(cx, cy) } else { zz },
    }
}

/// Lookup order of the fixed-base walk (see the module docs).
///
/// Radix-16 digits (39..54, 4 scans each) lead, so the multiplier starts
/// sooner, and are then spread among the radix-32 digits (8 scans each) so
/// scans and multiplies stay balanced. Found by searching orders with the
/// real scheduler (`tools/hostmodel`): 392 bundles, against 435 in digit
/// order. In [`fixed_base`], positions alternate between its two chains;
/// [`verify_known_key`] runs the whole order as one chain.
pub const LOOKUP_ORDER: [u8; DIGITS] = [
    39, 40, 41, 42, 43, 44, 0, 1, 2, 3, 45, 4, 5, 6, 46, 7, 8, 9, 10, 47, 11, 12, 13, 48, 14, 15,
    16, 17, 49, 18, 19, 20, 50, 21, 22, 23, 24, 51, 25, 26, 27, 52, 28, 29, 30, 31, 53, 32, 33, 34,
    35, 36, 37, 38,
];

/// The fixed-base walk `[s]B` over the digits in `order`, as an extended
/// point. `order` must have even length; even positions feed one chain and
/// odd positions the other.
fn fixed_base_point(b: &mut Builder, order: &[u8]) -> Point {
    let d2 = b.constant(super::host::edwards_d2());
    // 1/2 and 1/(2d), for `niels_to_ext`.
    let half = b.constant(Fp::ONE.add(Fp::ONE).invert());
    let inv_2d = b.constant(super::host::edwards_d2().invert());
    // [s]B = sum_i d_i (2^offset(i) B). Lookups alternate between two
    // independent chains (so the scheduler can fill both slots), each started
    // from its first lookup instead of adding to the identity. Radix-16
    // lookups (4 scans) are spread among the radix-32 ones (8 scans) so scans
    // and multiplies stay balanced along the whole program.
    //
    // `odd` takes positions 0, 2, 4, ... of `order` (odd when counted from
    // 1) and `even` positions 1, 3, 5, ...
    let l = b.lookup(order[0]);
    let mut odd = niels_to_ext(b, l, half, inv_2d);
    let l = b.lookup(order[1]);
    let mut even = niels_to_ext(b, l, half, inv_2d);
    for pair in order[2..].chunks(2) {
        let l = b.lookup(pair[0]);
        odd = madd(b, &odd, l);
        let l = b.lookup(pair[1]);
        even = madd(b, &even, l);
    }
    // Join the chains; only X, Y, Z are needed.
    add_ext(b, &odd, &even, d2, false)
}

/// [`fixed_base`] with an explicit lookup order (for schedule exploration).
#[doc(hidden)]
pub fn fixed_base_ordered(order: &[u8]) -> Program {
    let mut b = Builder::new();
    let p = fixed_base_point(&mut b, order);
    b.out(p.x, 0);
    b.out(p.y, 1);
    b.out(p.z, 2);
    b.compile("fixed_base")
}

/// Known-key verification core: `[S]B + [k](-A)` as two fixed-base walks.
///
/// Digit slots `0..54` hold `S`, slots `64..118` hold `k = H(R ‖ A ‖ M) mod l`,
/// and the key table holds [`point_table`](super::host::point_table)`(-A)`.
/// Outputs projective `X, Y, Z` (indices 0..3); the controller encodes the
/// point and compares it with `R`. `key_order` is the lookup order of the key
/// chain (the base chain uses [`LOOKUP_ORDER`]).
///
/// The two scalars give two independent chains, so this program does not
/// split each walk into odd and even chains as [`fixed_base`] does.
pub fn verify_known_key_ordered(key_order: &[u8]) -> Program {
    let mut b = Builder::new();
    let d2 = b.constant(super::host::edwards_d2());
    let half = b.constant(Fp::ONE.add(Fp::ONE).invert());
    let inv_2d = b.constant(super::host::edwards_d2().invert());
    // Each chain starts from its first lookup.
    let l = b.lookup(LOOKUP_ORDER[0]);
    let mut base = niels_to_ext(&mut b, l, half, inv_2d);
    let l = b.lookup_key(key_order[0]);
    let mut key = niels_to_ext(&mut b, l, half, inv_2d);
    // Base and key lookups alternate, keeping both chains in flight.
    for k in 1..DIGITS {
        let l = b.lookup(LOOKUP_ORDER[k]);
        base = madd(&mut b, &base, l);
        let l = b.lookup_key(key_order[k]);
        key = madd(&mut b, &key, l);
    }
    let p = add_ext(&mut b, &base, &key, d2, false);
    b.out(p.x, 0);
    b.out(p.y, 1);
    b.out(p.z, 2);
    b.compile("verify_known_key")
}

/// [`verify_known_key_ordered`] with the key chain in [`LOOKUP_ORDER`].
pub fn verify_known_key() -> Program {
    verify_known_key_ordered(&LOOKUP_ORDER)
}

/// Key loading, step 1: square-root candidate for a point with `y` in
/// input 0. Outputs `x` (0), `x * sqrt(-1)` (1), `v * x^2` (2) and `u` (3)
/// with `u = y^2 - 1`, `v = d y^2 + 1`. The controller takes `x` if
/// `v x^2 = u`, `x sqrt(-1)` if `v x^2 = -u`, and rejects the point
/// otherwise, then applies the sign bit; the program never branches.
///
/// The candidate is the RFC 8032 one, `x = u v^3 (u v^7)^((p-5)/8)`, which
/// computes a root of `u / v` without a separate inversion.
pub fn decode_candidate() -> Program {
    let mut b = Builder::new();
    let y = V::Input(0);
    let one = b.constant(Fp::ONE);
    let d = b.constant(super::host::edwards_d());
    let sqrt_m1 = b.constant(super::host::sqrt_m1());
    let y2 = b.sqr(y);
    let u = b.sub(y2, one);
    let dy2 = b.mul(y2, d);
    let v = b.add(dy2, one);
    // v^3 and v^7.
    let v2 = b.sqr(v);
    let v3 = b.mul(v2, v);
    let v6 = b.sqr(v3);
    let v7 = b.mul(v6, v);
    let uv3 = b.mul(u, v3);
    let uv7 = b.mul(u, v7);
    let t = b.pow_p58(uv7);
    let x = b.mul(uv3, t);
    // The other candidate, for when x^2 = -u/v.
    let xi = b.mul(x, sqrt_m1);
    // v x^2, which the controller compares with u and -u.
    let x2 = b.sqr(x);
    let vx2 = b.mul(v, x2);
    b.out(x, 0);
    b.out(xi, 1);
    b.out(vx2, 2);
    b.out(u, 3);
    b.compile("decode_candidate")
}

/// Key loading, step 2: the 54 table bases `2^offset(i) * P` of an affine
/// point `P = (x, y)` (inputs 0, 1), as projective `X, Y, Z, T`, in digit
/// order (output `4i + c`, counted by the controller).
///
/// Every output uses indices 0..4; the controller tells bases apart by
/// counting outputs, which arrive in program order.
pub fn table_bases() -> Program {
    let mut b = Builder::new();
    let one = b.constant(Fp::ONE);
    let x = V::Input(0);
    let y = V::Input(1);
    // P in extended coordinates with Z = 1, T = xy.
    let mut p = Point { x, y, z: one, t: b.mul(x, y) };
    for i in 0..DIGITS {
        b.out(p.x, 0);
        b.out(p.y, 1);
        b.out(p.z, 2);
        b.out(p.t, 3);
        if i + 1 < DIGITS {
            // Next base: digit i's width in doublings, computing T only on
            // the last one because the output needs it.
            let w = super::host::digit_width(i);
            for step in 0..w {
                p = double(&mut b, &p, step + 1 == w);
            }
        }
    }
    b.compile("table_bases")
}

/// Key loading, step 3: the `entries` table entries `e * Q` (`e = 1..`) of
/// a projective base `Q` (inputs 0..4 = `X, Y, Z, T`) in affine Niels form,
/// normalized with one shared inversion. Output `3(e-1) + c` is coordinate
/// `c` (`y+x`, `y-x`, `2dxy`) of entry `e`.
///
/// `entries` is 16 for a radix-32 digit and 8 for a radix-16 one. The
/// multiples are `Q`, `2Q` (one doubling), then repeated additions of `Q`;
/// their `Z`s are inverted together with Montgomery's trick as in
/// [`batch_compress`].
pub fn digit_table(entries: usize) -> Program {
    // The program always builds Q and 2Q before the rest.
    assert!(entries >= 2, "a digit table has at least two entries");
    let mut b = Builder::new();
    let d2 = b.constant(super::host::edwards_d2());
    let q = Point { x: V::Input(0), y: V::Input(1), z: V::Input(2), t: V::Input(3) };
    // m[e - 1] = e * Q, with T (the next addition reads it).
    let mut m = vec![q, double(&mut b, &q, true)];
    while m.len() < entries {
        let next = add_ext(&mut b, m.last().unwrap(), &q, d2, true);
        m.push(next);
    }
    // prefix[e] = Z_0 * ... * Z_e.
    let mut prefix = vec![m[0].z];
    for e in 1..entries {
        let p = b.mul(prefix[e - 1], m[e].z);
        prefix.push(p);
    }
    // inv = 1 / (Z_0 * ... * Z_e) for the current e, walking down.
    let mut inv = b.invert(prefix[entries - 1]);
    // Placeholders, all overwritten below.
    let mut zinv = vec![inv; entries];
    for e in (0..entries).rev() {
        if e == 0 {
            // inv is now 1 / Z_0.
            zinv[0] = inv;
        } else {
            // 1/Z_e = (Z_0 ... Z_{e-1}) / (Z_0 ... Z_e); then drop Z_e from
            // inv.
            zinv[e] = b.mul(inv, prefix[e - 1]);
            inv = b.mul(inv, m[e].z);
        }
    }
    for e in 0..entries {
        // Affine x, y, then the Niels coordinates.
        let x = b.mul(m[e].x, zinv[e]);
        let y = b.mul(m[e].y, zinv[e]);
        let ypx = b.add(y, x);
        let ymx = b.sub(y, x);
        let xy = b.mul(x, y);
        let xy2d = b.mul(xy, d2);
        b.out(ypx, (3 * e) as u8);
        b.out(ymx, (3 * e + 1) as u8);
        b.out(xy2d, (3 * e + 2) as u8);
    }
    b.compile("digit_table")
}

/// Fixed-base multiplication followed by its own inversion: outputs affine
/// `x` (index 0) and `y` (index 1). Used when no batch-compress machine exists.
pub fn fixed_base_affine() -> Program {
    let mut b = Builder::new();
    let p = fixed_base_point(&mut b, &LOOKUP_ORDER);
    let zi = b.invert(p.z);
    let x = b.mul(p.x, zi);
    let y = b.mul(p.y, zi);
    b.out(x, 0);
    b.out(y, 1);
    b.compile("fixed_base_affine")
}

/// Fixed-base scalar multiplication. Outputs projective `X, Y, Z` (indices 0..3).
///
/// `[s]B = sum_i d_i T_i` with `T_i = 2^offset(i) B`, split into two sums by
/// position in [`LOOKUP_ORDER`]: digits at even positions form one chain and
/// digits at odd positions the other. The two chains are independent and
/// interleaved, so the scheduler can fill the bundle slots that one chain's
/// add/multiply recurrence would leave empty; one final addition joins them.
pub fn fixed_base() -> Program {
    let mut b = Builder::new();
    let p = fixed_base_point(&mut b, &LOOKUP_ORDER);
    b.out(p.x, 0);
    b.out(p.y, 1);
    b.out(p.z, 2);
    b.compile("fixed_base")
}

/// Batch affine conversion of `g` points held in bank I as `X, Y, Z` triples.
///
/// Output `2p` is affine `x` and `2p + 1` is affine `y` of point `p`. All
/// values are canonical, so `x`'s least-significant bit is its Ed25519 sign.
/// A zero `Z` would zero the whole batch; valid curve points never have one.
///
/// Montgomery's trick: `g - 1` multiplies build the prefix products, one
/// inversion inverts the last, and a backward walk recovers each `1/Z` with
/// two multiplies.
pub fn batch_compress(g: u8) -> Program {
    assert!(g >= 1);
    let mut b = Builder::new();
    // Point p's coordinates in bank I.
    let x = |p: u8| V::Input(3 * p);
    let y = |p: u8| V::Input(3 * p + 1);
    let z = |p: u8| V::Input(3 * p + 2);
    // prefix[p] = Z_0 * ... * Z_p.
    let mut prefix = vec![z(0)];
    for p in 1..g {
        let next = b.mul(prefix[p as usize - 1], z(p));
        prefix.push(next);
    }
    // inv = 1 / (Z_0 * ... * Z_p) for the current p, walking down.
    let mut inv = b.invert(prefix[g as usize - 1]);
    for p in (0..g).rev() {
        let zi = if p == 0 {
            inv
        } else {
            // 1/Z_p, then drop Z_p from inv.
            let zi = b.mul(inv, prefix[p as usize - 1]);
            inv = b.mul(inv, z(p));
            zi
        };
        let ax = b.mul(x(p), zi);
        let ay = b.mul(y(p), zi);
        b.out(ax, 2 * p);
        b.out(ay, 2 * p + 1);
    }
    // Adds zero as C[0], the `Out` dummy operand. The program uses no other
    // constant, and `compile` would add the same zero to an empty bank, so
    // this only makes the entry explicit.
    let _ = b.constant(Fp::ZERO);
    b.compile("batch_compress")
}

/// Interpreter tests of every program against Dalek and the host model.
#[cfg(test)]
mod tests {
    use super::super::host::{Ext, Fp, L_BYTES, recode};
    use super::super::microcode::{Context, interpret};
    use super::*;

    /// Reduce a 64-byte little-endian value modulo `l`.
    fn reduce_l(bytes: [u8; 64]) -> [u8; 32] {
        use num_bigint::BigUint;
        let l = BigUint::from_bytes_le(&L_BYTES);
        let v = BigUint::from_bytes_le(&bytes) % l;
        let mut out = [0u8; 32];
        let b = v.to_bytes_le();
        out[..b.len()].copy_from_slice(&b);
        out
    }

    /// `n` deterministic pseudo-random scalars below `l` (xorshift64 bytes,
    /// reduced from 64 bytes like a hash).
    fn scalars(n: usize) -> Vec<[u8; 32]> {
        let mut state = 0x1234_5678_9abc_def0u64;
        (0..n)
            .map(|_| {
                let mut wide = [0u8; 64];
                for byte in &mut wide {
                    state ^= state << 13;
                    state ^= state >> 7;
                    state ^= state << 17;
                    *byte = state as u8;
                }
                reduce_l(wide)
            })
            .collect()
    }

    /// `[S]B + [k](-A)` equals `R` exactly when `R = [S]B - [k]A`; checked
    /// for random keys and scalars, plus one corrupted case.
    #[test]
    fn verify_known_key_computes_the_equation() {
        use super::super::host::point_table;
        let p = verify_known_key();
        println!(
            "verify_known_key: {} bundles, {} mul, {} add, {} scan, regs M{} A{}",
            p.bundles.len(), p.mul_ops, p.add_ops, p.scan_ops, p.regs_m, p.regs_a
        );
        let all = scalars(12);
        for t in 0..4 {
            // Secret key scalar a, signature scalar s, challenge k.
            let (a, s, k) = (all[3 * t], all[3 * t + 1], all[3 * t + 2]);
            let public = Ext::basepoint().mul_bytes(&a);
            // The R that makes (R, s) verify for challenge k.
            let r = Ext::basepoint().mul_bytes(&s).add(public.mul_bytes(&k).neg());
            let table = point_table(public.neg());
            let mut ctx = Context::with_key(&p, &[], recode(&s), recode(&k), table.clone());
            let out = interpret(&p, &mut ctx);
            // `compress` does not read T.
            let point = Ext { x: out[0].1, y: out[1].1, z: out[2].1, t: Fp::ZERO };
            assert_eq!(point.compress(), r.compress(), "case {t}");
            // A different k must not verify.
            let mut ctx = Context::with_key(&p, &[], recode(&s), recode(&a), table);
            let out = interpret(&p, &mut ctx);
            let point = Ext { x: out[0].1, y: out[1].1, z: out[2].1, t: Fp::ZERO };
            assert_ne!(point.compress(), r.compress(), "case {t} with the wrong k");
        }
    }

    /// Key loading end to end on the interpreter (with the controller's
    /// root choice done here): the built table equals `point_table(-A)`.
    #[test]
    fn key_load_programs_build_the_key_table() {
        use super::super::host::{decode_point, digit_entries, point_table};
        let dec = decode_candidate();
        let bases = table_bases();
        let wide = digit_table(16);
        let narrow = digit_table(8);
        for (name, p) in [("decode_candidate", &dec), ("table_bases", &bases), ("digit_table(16)", &wide), ("digit_table(8)", &narrow)] {
            println!("{name}: {} bundles, {} mul, regs M{} A{} I{}", p.bundles.len(), p.mul_ops, p.regs_m, p.regs_a, p.regs_i);
        }
        for a in scalars(2) {
            let public = Ext::basepoint().mul_bytes(&a);
            let enc = public.compress();
            // Split the encoding into the sign of x and y.
            let sign = enc[31] >> 7 == 1;
            let mut yb = enc;
            yb[31] &= 0x7f;
            let y = Fp::from_bytes(&yb);
            // Step 1 + controller root choice.
            let mut ctx = Context::new(&dec, &[y], [0; 64]);
            let out = interpret(&dec, &mut ctx);
            let (x, xi, vx2, u) = (out[0].1, out[1].1, out[2].1, out[3].1);
            let mut root = if vx2 == u { x } else if vx2 == u.neg() { xi } else { panic!("not on the curve") };
            if root.is_negative() != sign {
                root = root.neg();
            }
            assert_eq!(Some(root), decode_point(&enc).map(|p| p.x));
            // -A = (-x, y).
            let neg_x = root.neg();
            let mut ctx = Context::new(&bases, &[neg_x, y], [0; 64]);
            let out = interpret(&bases, &mut ctx);
            let want = point_table(public.neg());
            for i in 0..DIGITS {
                // Base i is outputs 4i..4i+4 (X, Y, Z, T).
                let q: Vec<Fp> = (0..4).map(|c| out[4 * i + c].1).collect();
                let prog = if digit_entries(i) == 16 { &wide } else { &narrow };
                let mut ctx = Context::new(prog, &q, [0; 64]);
                let entries = interpret(prog, &mut ctx);
                for e in 0..digit_entries(i) {
                    let n = super::super::host::Niels { ypx: entries[3 * e].1, ymx: entries[3 * e + 1].1, xy2d: entries[3 * e + 2].1 };
                    assert_eq!(n, want[i][e], "digit {i} entry {e}");
                }
            }
        }
    }

    /// The fused program's affine output encodes to Dalek's `[s]B`.
    #[test]
    fn fused_fixed_base_matches_dalek() {
        use curve25519_dalek::{EdwardsPoint, Scalar};
        let p = fixed_base_affine();
        println!(
            "fixed_base_affine: {} bundles, {} mul, regs M{} A{}",
            p.bundles.len(),
            p.mul_ops,
            p.regs_m,
            p.regs_a
        );
        for s in scalars(6) {
            let mut ctx = Context::new(&p, &[], recode(&s));
            let out = interpret(&p, &mut ctx);
            let (x, y) = (out[0].1, out[1].1);
            // Ed25519 encoding: y with x's sign in bit 255.
            let mut enc = y.to_bytes();
            enc[31] |= u8::from(x.is_negative()) << 7;
            let expect = EdwardsPoint::mul_base(&Scalar::from_bytes_mod_order(s)).compress();
            assert_eq!(enc, expect.to_bytes());
        }
    }

    /// `fixed_base` followed by `batch_compress(4)` encodes to Dalek's
    /// `[s]B`, including the scalars 0, 1 and `l - 1`.
    #[test]
    fn fixed_base_and_batch_compress_match_dalek() {
        use curve25519_dalek::{EdwardsPoint, Scalar};
        let fb = fixed_base();
        let g = 4;
        let cp = batch_compress(g);
        println!(
            "fixed_base: {} bundles, {} mul, {} add, {} scan, regs M{} A{}",
            fb.bundles.len(),
            fb.mul_ops,
            fb.add_ops,
            fb.scan_ops,
            fb.regs_m,
            fb.regs_a
        );
        println!(
            "batch_compress({g}): {} bundles, {} mul, regs M{} I{}",
            cp.bundles.len(),
            cp.mul_ops,
            cp.regs_m,
            cp.regs_i
        );
        // Edge-case scalars: 0, 1 and l - 1.
        let mut special = vec![[0u8; 32], {
            let mut one = [0u8; 32];
            one[0] = 1;
            one
        }];
        let mut lm1 = L_BYTES;
        lm1[0] -= 1;
        special.push(lm1);
        // 11 scalars: two full batches of 4 and a short batch of 3, so the
        // identity padding below runs too.
        let all: Vec<[u8; 32]> = special.into_iter().chain(scalars(8)).collect();
        for chunk in all.chunks(g as usize) {
            let mut inputs = Vec::new();
            for s in chunk {
                let mut ctx = Context::new(&fb, &[], recode(s));
                let out = interpret(&fb, &mut ctx);
                assert_eq!(out.len(), 3);
                inputs.extend(out.iter().map(|(_, v)| *v));
            }
            // Pad a short batch with the identity.
            while inputs.len() < 3 * g as usize {
                inputs.extend([Fp::ZERO, Fp::ONE, Fp::ONE]);
            }
            let mut ctx = Context::new(&cp, &inputs, [0; 64]);
            let out = interpret(&cp, &mut ctx);
            for (p, s) in chunk.iter().enumerate() {
                let x = out.iter().find(|(i, _)| *i as usize == 2 * p).unwrap().1;
                let y = out.iter().find(|(i, _)| *i as usize == 2 * p + 1).unwrap().1;
                let mut enc = y.to_bytes();
                enc[31] |= u8::from(x.is_negative()) << 7;
                let expect = EdwardsPoint::mul_base(&Scalar::from_bytes_mod_order(*s)).compress();
                assert_eq!(enc, expect.to_bytes());
                assert_eq!(enc, Ext::basepoint().mul_bytes(s).compress());
            }
        }
    }
}