rhdl_ed25519_scalable_core/

host.rs

//! Host-side reference arithmetic for the scalable datapath.
//!
//! These functions run on the build/test host, never in hardware. They
//! generate ROM contents (the fixed-base table and constants) and give the
//! microcode interpreter bit-exact semantics: every field value is canonical,
//! strictly below `p = 2^255 - 19`, exactly as the RHDL multiplier and
//! add/subtract units produce.
//!
//! The module also holds the reference curve model ([`Ext`], [`Niels`]),
//! point encoding and strict decoding, the small-order encoding list, and the
//! mixed-radix scalar recoding ([`recode`]) with its digit layout
//! ([`DIGITS`], [`digit_offset`], [`table_row`]). None of it is constant
//! time; it only computes reference values.

use std::sync::OnceLock;

/// Canonical field element modulo `2^255 - 19`, four little-endian words.
///
/// Every operation returns a fully reduced value, so `==` on the words is
/// field equality and the low bit is the Ed25519 sign.
#[derive(Copy, Clone, Debug, PartialEq, Eq, Default)]
pub struct Fp(pub [u64; 4]);

/// The field prime `p = 2^255 - 19` as four little-endian words.
const P: [u64; 4] = [
    0xffff_ffff_ffff_ffed,
    0xffff_ffff_ffff_ffff,
    0xffff_ffff_ffff_ffff,
    0x7fff_ffff_ffff_ffff,
];

/// `a >= b` for four-word little-endian integers.
fn ge(a: &[u64; 4], b: &[u64; 4]) -> bool {
    // Compare from the most significant word down.
    for i in (0..4).rev() {
        if a[i] != b[i] {
            return a[i] > b[i];
        }
    }
    true
}

/// `a - b mod 2^256` for four-word little-endian integers.
fn sub_raw(a: &[u64; 4], b: &[u64; 4]) -> [u64; 4] {
    let mut out = [0u64; 4];
    let mut borrow = 0u64;
    for i in 0..4 {
        // Borrow out of this word if either subtraction wrapped.
        let (d, b1) = a[i].overflowing_sub(b[i]);
        let (d, b2) = d.overflowing_sub(borrow);
        out[i] = d;
        borrow = u64::from(b1 || b2);
    }
    out
}

impl Fp {
    /// Zero.
    pub const ZERO: Fp = Fp([0; 4]);
    /// One.
    pub const ONE: Fp = Fp([1, 0, 0, 0]);

    /// Reduce an arbitrary 512-bit little-endian value.
    fn reduce512(w: [u64; 8]) -> Fp {
        // 2^256 = 38 (mod p).
        // Fold the high 256 bits into the low ones: acc[i] = w[i] + 38 w[i+4]
        // (acc[4] stays zero).
        let mut acc = [0u128; 5];
        for i in 0..4 {
            acc[i] += u128::from(w[i]) + u128::from(w[i + 4]) * 38;
        }
        // Propagate carries into 64-bit limbs; limbs[4] is small (below 40).
        let mut limbs = [0u64; 5];
        let mut carry = 0u128;
        for i in 0..4 {
            let v = acc[i] + carry;
            limbs[i] = v as u64;
            carry = v >> 64;
        }
        limbs[4] = carry as u64;
        // Fold bit 255 and above: value = lo255 + top * 2^255, 2^255 = 19.
        let mut v = [limbs[0], limbs[1], limbs[2], limbs[3]];
        let top = (u128::from(limbs[4]) << 1) | u128::from(v[3] >> 63);
        v[3] &= 0x7fff_ffff_ffff_ffff;
        let mut c = top * 19;
        for limb in &mut v {
            let s = u128::from(*limb) + c;
            *limb = s as u64;
            c = s >> 64;
        }
        // The value is now below 2^255 + a small amount; subtract p until
        // canonical.
        let mut out = Fp(v);
        while ge(&out.0, &P) {
            out = Fp(sub_raw(&out.0, &P));
        }
        out
    }

    /// Canonical value from 32 little-endian bytes (reduced modulo `p`).
    ///
    /// Bit 255 is included in the value, so callers decoding a point must
    /// clear the sign bit first.
    pub fn from_bytes(bytes: &[u8; 32]) -> Fp {
        let mut w = [0u64; 8];
        for i in 0..4 {
            w[i] = u64::from_le_bytes(bytes[i * 8..i * 8 + 8].try_into().unwrap());
        }
        Fp::reduce512(w)
    }

    /// Canonical little-endian encoding.
    pub fn to_bytes(self) -> [u8; 32] {
        let mut out = [0u8; 32];
        for i in 0..4 {
            out[i * 8..i * 8 + 8].copy_from_slice(&self.0[i].to_le_bytes());
        }
        out
    }

    /// Sum modulo `p`.
    pub fn add(self, o: Fp) -> Fp {
        // 257-bit sum, reduced by the general path.
        let mut w = [0u64; 8];
        let mut carry = 0u128;
        for i in 0..4 {
            let s = u128::from(self.0[i]) + u128::from(o.0[i]) + carry;
            w[i] = s as u64;
            carry = s >> 64;
        }
        w[4] = carry as u64;
        Fp::reduce512(w)
    }

    /// Difference modulo `p`.
    pub fn sub(self, o: Fp) -> Fp {
        // -o as p - o, except that -0 is 0 (p itself is not canonical).
        let neg = if o == Fp::ZERO {
            Fp::ZERO
        } else {
            Fp(sub_raw(&P, &o.0))
        };
        self.add(neg)
    }

    /// Product modulo `p`.
    pub fn mul(self, o: Fp) -> Fp {
        // Schoolbook 256 x 256 -> 512-bit product, one 64-bit word per entry
        // of `w` (u128 entries hold intermediate sums).
        let mut w = [0u128; 8];
        for i in 0..4 {
            let mut carry = 0u128;
            for j in 0..4 {
                let t = w[i + j] + u128::from(self.0[i]) * u128::from(o.0[j]) + carry;
                w[i + j] = t & u128::from(u64::MAX);
                carry = t >> 64;
            }
            w[i + 4] += carry;
        }
        Fp::reduce512(w.map(|x| x as u64))
    }

    /// Square modulo `p`.
    pub fn square(self) -> Fp {
        self.mul(self)
    }

    /// `self^(2^k)`.
    pub fn pow2k(self, k: u32) -> Fp {
        let mut x = self;
        for _ in 0..k {
            x = x.square();
        }
        x
    }

    /// Multiplicative inverse (zero maps to zero).
    ///
    /// Fermat's little theorem, `self^(p-2)`, by plain square-and-multiply.
    pub fn invert(self) -> Fp {
        // p - 2 = 2^255 - 21.
        let mut result = Fp::ONE;
        let e: [u64; 4] = sub_raw(&P, &[2, 0, 0, 0]);
        // Exponent bits from the top (bit 254) down.
        for i in (0..255).rev() {
            result = result.square();
            if (e[i / 64] >> (i % 64)) & 1 == 1 {
                result = result.mul(self);
            }
        }
        result
    }

    /// Negation modulo `p`.
    pub fn neg(self) -> Fp {
        Fp::ZERO.sub(self)
    }

    /// Least-significant bit of the canonical value (Ed25519 "negative").
    pub fn is_negative(self) -> bool {
        self.0[0] & 1 == 1
    }

    /// Fifteen 17-bit digits of the canonical value, least significant first.
    ///
    /// This is the digit layout of the RHDL field units (15 x 17 = 255 bits).
    pub fn digits17(self) -> [u32; 15] {
        let mut out = [0u32; 15];
        for (k, digit) in out.iter_mut().enumerate() {
            let mut v = 0u32;
            for b in 0..17 {
                // Absolute bit index of bit `b` of digit `k`.
                let bit = 17 * k + b;
                v |= (((self.0[bit / 64] >> (bit % 64)) & 1) as u32) << b;
            }
            *digit = v;
        }
        out
    }

    /// Inverse of [`Self::digits17`]; the digits must encode a value below `p`.
    ///
    /// The value is not reduced, so out-of-range digits give a non-canonical
    /// `Fp`.
    pub fn from_digits17(d: &[u32; 15]) -> Fp {
        let mut w = [0u64; 4];
        for (k, digit) in d.iter().enumerate() {
            for b in 0..17 {
                if (digit >> b) & 1 == 1 {
                    let bit = 17 * k + b;
                    w[bit / 64] |= 1 << (bit % 64);
                }
            }
        }
        Fp(w)
    }

    /// A small integer as a field element (canonical because `v < p`).
    fn from_u64(v: u64) -> Fp {
        Fp([v, 0, 0, 0])
    }
}

/// Edwards curve constant `d = -121665/121666`.
pub fn edwards_d() -> Fp {
    Fp::from_u64(121665)
        .neg()
        .mul(Fp::from_u64(121666).invert())
}

/// `2d`, the constant the addition formulas and the Niels form use.
pub fn edwards_d2() -> Fp {
    let d = edwards_d();
    d.add(d)
}

/// `sqrt(-1)`, i.e. `2^((p-1)/4)`.
///
/// 2 is a non-square modulo `p`, so `2^((p-1)/2) = -1` and its square root
/// `2^((p-1)/4)` is a square root of `-1`.
pub fn sqrt_m1() -> Fp {
    // (p - 1)/4 = 2^253 - 5.
    let mut result = Fp::ONE;
    let two = Fp::from_u64(2);
    let e: [u64; 4] = [
        0xffff_ffff_ffff_fffb,
        0xffff_ffff_ffff_ffff,
        0xffff_ffff_ffff_ffff,
        0x1fff_ffff_ffff_ffff,
    ];
    // Square-and-multiply from bit 252 down.
    for i in (0..253).rev() {
        result = result.square();
        if (e[i / 64] >> (i % 64)) & 1 == 1 {
            result = result.mul(two);
        }
    }
    result
}

/// Extended twisted-Edwards point `(X:Y:Z:T)`.
///
/// Affine `x = X/Z`, `y = Y/Z`, on `-x^2 + y^2 = 1 + d x^2 y^2`.
#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub struct Ext {
    /// X.
    pub x: Fp,
    /// Y.
    pub y: Fp,
    /// Z.
    pub z: Fp,
    /// T = XY/Z.
    pub t: Fp,
}

/// Affine Niels coordinates `(y+x, y-x, 2dxy)`.
///
/// The form the tables store: a mixed addition with it needs no `Z` and no
/// multiply by `2d`, and negation only swaps two coordinates and negates the
/// third.
#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub struct Niels {
    /// `y + x`.
    pub ypx: Fp,
    /// `y - x`.
    pub ymx: Fp,
    /// `2dxy`.
    pub xy2d: Fp,
}

impl Niels {
    /// Niels form of the identity `(0, 1)`: `(1, 1, 0)`.
    pub fn identity() -> Self {
        Niels {
            ypx: Fp::ONE,
            ymx: Fp::ONE,
            xy2d: Fp::ZERO,
        }
    }
    /// Negated point: `-(x, y) = (-x, y)` swaps `y+x` and `y-x` and negates
    /// `2dxy`.
    pub fn neg(self) -> Self {
        Niels {
            ypx: self.ymx,
            ymx: self.ypx,
            xy2d: self.xy2d.neg(),
        }
    }
}

impl Ext {
    /// Negated point `(-X:Y:Z:-T)`.
    pub fn neg(self) -> Ext {
        Ext {
            x: self.x.neg(),
            y: self.y,
            z: self.z,
            t: self.t.neg(),
        }
    }

    /// Identity `(0:1:1:0)`.
    pub fn identity() -> Self {
        Ext {
            x: Fp::ZERO,
            y: Fp::ONE,
            z: Fp::ONE,
            t: Fp::ZERO,
        }
    }

    /// Ed25519 basepoint: `y = 4/5`, `x` non-negative (even).
    pub fn basepoint() -> Self {
        let y = Fp::from_u64(4).mul(Fp::from_u64(5).invert());
        let x = recover_x(y, false).expect("basepoint");
        Ext {
            x,
            y,
            z: Fp::ONE,
            t: x.mul(y),
        }
    }

    /// Complete unified addition (Hisil et al., a = -1).
    ///
    /// Complete means it is correct for every pair of inputs, including
    /// doubling and the identity, which keeps this reference model simple.
    pub fn add(self, o: Ext) -> Ext {
        let d2 = edwards_d2();
        let a = self.y.sub(self.x).mul(o.y.sub(o.x));
        let b = self.y.add(self.x).mul(o.y.add(o.x));
        let c = self.t.mul(d2).mul(o.t);
        let d = self.z.add(self.z).mul(o.z);
        let e = b.sub(a);
        let f = d.sub(c);
        let g = d.add(c);
        let h = b.add(a);
        Ext {
            x: e.mul(f),
            y: g.mul(h),
            z: f.mul(g),
            t: e.mul(h),
        }
    }

    /// Doubling, through the complete addition.
    pub fn double(self) -> Ext {
        self.add(self)
    }

    /// Affine coordinates `(X/Z, Y/Z)`.
    pub fn affine(self) -> (Fp, Fp) {
        let zi = self.z.invert();
        (self.x.mul(zi), self.y.mul(zi))
    }

    /// Standard 32-byte compressed encoding: `y` little endian with the sign
    /// of `x` in bit 255.
    pub fn compress(self) -> [u8; 32] {
        let (x, y) = self.affine();
        let mut s = y.to_bytes();
        // Bit 255 of a canonical y is zero, so XOR sets it.
        s[31] ^= u8::from(x.is_negative()) << 7;
        s
    }

    /// Affine Niels form.
    pub fn niels(self) -> Niels {
        let (x, y) = self.affine();
        Niels {
            ypx: y.add(x),
            ymx: y.sub(x),
            xy2d: edwards_d2().mul(x).mul(y),
        }
    }

    /// Scalar multiplication by a little-endian 256-bit integer (host only).
    ///
    /// Double-and-add from the top bit, not constant time.
    pub fn mul_bytes(self, k: &[u8; 32]) -> Ext {
        let mut acc = Ext::identity();
        for i in (0..256).rev() {
            acc = acc.double();
            if (k[i / 8] >> (i % 8)) & 1 == 1 {
                acc = acc.add(self);
            }
        }
        acc
    }
}

/// The `x` with the given sign for which `(x, y)` is on the curve, or
/// `None` if `y` has no such `x`.
///
/// From the curve equation, `x^2 = (y^2 - 1) / (d y^2 + 1)`. Returns
/// `Some(0)` for `x = 0` even when `negative` is set; [`decode_point`]
/// rejects that encoding.
fn recover_x(y: Fp, negative: bool) -> Option<Fp> {
    let d = edwards_d();
    let yy = y.square();
    // x^2 = u / v.
    let u = yy.sub(Fp::ONE);
    let v = d.mul(yy).add(Fp::ONE);
    let x2 = u.mul(v.invert());
    // x = x2^((p+3)/8), (p+3)/8 = 2^252 - 2.
    // Because p = 5 (mod 8), this is a square root of x2 or of -x2.
    let mut x = Fp::ONE;
    let e: [u64; 4] = [
        0xffff_ffff_ffff_fffe,
        0xffff_ffff_ffff_ffff,
        0xffff_ffff_ffff_ffff,
        0x0fff_ffff_ffff_ffff,
    ];
    for i in (0..252).rev() {
        x = x.square();
        if (e[i / 64] >> (i % 64)) & 1 == 1 {
            x = x.mul(x2);
        }
    }
    // A root of -x2: multiply by sqrt(-1) to get a root of x2.
    if x.square() != x2 {
        x = x.mul(sqrt_m1());
    }
    // Still not a root: x2 is not a square, so no point has this y.
    if x.square() != x2 {
        return None;
    }
    // Pick the root whose low bit matches the requested sign.
    if x.is_negative() != negative {
        x = x.neg();
    }
    Some(x)
}

/// Digits of the fixed-base recoding (one table per digit).
pub const DIGITS: usize = 54;
/// The low `WIDE_DIGITS` digits are signed radix 32, the rest radix 16.
///
/// 39 * 5 + 15 * 4 = 255 bits. The split balances the machine's slots: a
/// radix-32 lookup scans 8 entry pairs and a radix-16 one 4, every digit
/// costs one 7-multiply mixed addition, and with this split scans and
/// multiplies both total about 375 per fixed-base multiplication.
pub const WIDE_DIGITS: usize = 39;
/// Table rows (entry pairs) over all digits: 39 * 8 + 15 * 4 = 372.
pub const TABLE_ROWS: usize = WIDE_DIGITS * 8 + (DIGITS - WIDE_DIGITS) * 4;

/// Bit width of digit `i`: 5 for radix 32, 4 for radix 16.
pub const fn digit_width(i: usize) -> u32 {
    if i < WIDE_DIGITS { 5 } else { 4 }
}

/// Bit offset of digit `i`: the sum of the widths of digits `0..i`.
pub const fn digit_offset(i: usize) -> u32 {
    if i < WIDE_DIGITS {
        5 * i as u32
    } else {
        5 * WIDE_DIGITS as u32 + 4 * (i - WIDE_DIGITS) as u32
    }
}

/// Table entries of digit `i` (magnitudes `1..=entries`).
///
/// A signed width-`w` digit has magnitude at most `2^(w-1)`: 16 entries for
/// radix 32, 8 for radix 16. Zero needs no entry (it selects the identity).
pub const fn digit_entries(i: usize) -> usize {
    1 << (digit_width(i) - 1)
}

/// First table row (entry pair) of digit `i`.
///
/// Rows are packed digit after digit: 8 per radix-32 digit, 4 per radix-16
/// digit.
pub const fn table_row(i: usize) -> usize {
    if i < WIDE_DIGITS {
        8 * i
    } else {
        8 * WIDE_DIGITS + 4 * (i - WIDE_DIGITS)
    }
}

/// `table[i][e-1] = e * 2^offset(i) * B` in affine Niels form.
///
/// Built once on first use (it takes many host-side inversions) and shared.
pub fn basepoint_table() -> &'static [Vec<Niels>] {
    static TABLE: OnceLock<Vec<Vec<Niels>>> = OnceLock::new();
    TABLE.get_or_init(|| point_table(Ext::basepoint()))
}

/// `table[i][e-1] = e * 2^offset(i) * P` in affine Niels form: the
/// fixed-base table of any point (the base point, or a verification key).
pub fn point_table(p: Ext) -> Vec<Vec<Niels>> {
    // 2^offset(i) * P for the current digit.
    let mut base = p;
    let mut table = Vec::with_capacity(DIGITS);
    for i in 0..DIGITS {
        let mut row = Vec::with_capacity(digit_entries(i));
        // e * base, for e = 1, 2, ...
        let mut multiple = base;
        for _ in 0..digit_entries(i) {
            row.push(multiple.niels());
            multiple = multiple.add(base);
        }
        table.push(row);
        // Advance to the next digit's base.
        for _ in 0..digit_width(i) {
            base = base.double();
        }
    }
    table
}

/// Whether a 32-byte encoding has a canonical `y` (below `p`, sign bit
/// ignored).
pub fn canonical_y(bytes: &[u8; 32]) -> bool {
    let mut y = *bytes;
    y[31] &= 0x7f;
    // Compare with p = 2^255 - 19, most significant byte first.
    let mut p = [0xffu8; 32];
    p[0] = 0xed;
    p[31] = 0x7f;
    for i in (0..32).rev() {
        if y[i] != p[i] {
            return y[i] < p[i];
        }
    }
    // y == p is not canonical.
    false
}

/// Strict point decoding: canonical `y`, a square root exists, and not the
/// `x = 0` encoding with the sign bit set. Same rules as the strict
/// verifier of `rhdl_ed25519_fast_core`.
pub fn decode_point(bytes: &[u8; 32]) -> Option<Ext> {
    if !canonical_y(bytes) {
        return None;
    }
    // Bit 255 is the sign of x; the rest is y.
    let negative = bytes[31] >> 7 == 1;
    let mut yb = *bytes;
    yb[31] &= 0x7f;
    let y = Fp::from_bytes(&yb);
    let x = recover_x(y, negative)?;
    // "Negative zero" has no point of its own; reject it.
    if x == Fp::ZERO && negative {
        return None;
    }
    Some(Ext {
        x,
        y,
        z: Fp::ONE,
        t: x.mul(y),
    })
}

/// The eight encodings of low-order points that strict verification
/// rejects (the list `rhdl_ed25519_fast_core` checks).
///
/// These are the canonical encodings of the eight points of the torsion
/// subgroup (order dividing 8), compared as four little-endian 64-bit
/// words. With a small-order key, `[k]A` takes at most eight values, so a
/// signature can be forged without the secret key; strict verifiers refuse
/// these encodings.
pub fn small_order_encoding(bytes: &[u8; 32]) -> bool {
    const WORDS: [[u64; 4]; 8] = [
        // y = 1: the identity (order 1).
        [1, 0, 0, 0],
        // An order-8 point, x non-negative.
        [0x4fd84d3d706a17c7, 0x0f67100d760b3cba, 0xc6cc392cfa53202a, 0x7a03ac9277fdc74e],
        // y = 0 with the sign bit set: the order-4 point whose x is odd.
        [0, 0, 0, 0x8000_0000_0000_0000],
        // An order-8 point, x non-negative.
        [0xb027b2c28f95e826, 0xf098eff289f4c345, 0x3933c6d305acdfd5, 0x05fc536d880238b1],
        // y = p - 1 = -1: the point (0, -1) (order 2).
        [0xffff_ffff_ffff_ffec, 0xffff_ffff_ffff_ffff, 0xffff_ffff_ffff_ffff, 0x7fff_ffff_ffff_ffff],
        // The previous order-8 y with the sign bit set.
        [0xb027b2c28f95e826, 0xf098eff289f4c345, 0x3933c6d305acdfd5, 0x85fc536d880238b1],
        // y = 0 with the sign bit clear: the order-4 point whose x is even.
        [0, 0, 0, 0],
        // The first order-8 y with the sign bit set.
        [0x4fd84d3d706a17c7, 0x0f67100d760b3cba, 0xc6cc392cfa53202a, 0xfa03ac9277fdc74e],
    ];
    let words: [u64; 4] = std::array::from_fn(|k| u64::from_le_bytes(bytes[8 * k..8 * k + 8].try_into().unwrap()));
    WORDS.contains(&words)
}

/// Signed mixed-radix recoding into [`DIGITS`] digits (the rest of the
/// array is zero): digit `i` lies in `[-2^(w-1), 2^(w-1)]` for its width `w`.
///
/// The scalar must be below `2^253` (every scalar the signer multiplies is
/// reduced modulo `l`), which keeps the top digit in `[0, 4]`.
///
/// Signed digits halve the table: magnitudes run to `2^(w-1)` instead of
/// `2^w - 1`, and the sign is applied by negating the selected entry. The
/// array has 64 slots to match the machine's per-scalar digit slots.
pub fn recode(scalar: &[u8; 32]) -> [i8; 64] {
    assert!(scalar[31] < 32, "scalar must be below 2^253");
    // Bit `k` of the scalar.
    let bit = |k: u32| -> i32 { ((scalar[(k / 8) as usize] >> (k % 8)) & 1) as i32 };
    let mut d = [0i8; 64];
    // Carry out of the previous digit (0 or 1).
    let mut carry = 0;
    for i in 0..DIGITS {
        let w = digit_width(i);
        // The digit's w bits plus the incoming carry, in 0..=2^w.
        let mut v = carry;
        for b in 0..w {
            v += bit(digit_offset(i) + b) << b;
        }
        // Values from 2^(w-1) up become v - 2^w (in [-2^(w-1), 0]) and carry
        // 1 into the next digit. The top digit keeps its value: the scalar
        // bound keeps it at most 4.
        if i + 1 < DIGITS && v >= 1 << (w - 1) {
            d[i] = (v - (1 << w)) as i8;
            carry = 1;
        } else {
            d[i] = v as i8;
            carry = 0;
        }
    }
    d
}

/// Ed25519 group order `l = 2^252 + 27742317777372353535851937790883648493`,
/// little endian.
pub const L_BYTES: [u8; 32] = [
    0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0x10,
];