ic-ec 0.3.0

//! Arithmetic modulo the Curve25519 group order.
//!
//! `L = 2^252 + 27742317777372353535851937790883648493`
//!
//! Ed25519 signing needs `(r + h*a) mod L` where every operand is secret, so
//! reduction must be constant-time. This module uses a fixed 512-iteration
//! shift-and-conditional-subtract loop: slower than a Barrett or Montgomery
//! reduction, but with no data-dependent control flow and short enough to
//! audit by reading it.

/// The group order `L`, little-endian.
pub const L: [u8; 32] = [
    0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58, 0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10,
];

/// `L` as nine 32-bit little-endian limbs (the top limb is always zero).
const L_LIMBS: [u32; 9] = [
    0x5cf5_d3ed,
    0x5812_631a,
    0xa2f7_9cd6,
    0x14de_f9de,
    0x0000_0000,
    0x0000_0000,
    0x0000_0000,
    0x1000_0000,
    0x0000_0000,
];

/// Subtract `L` from `r` when `r >= L`, in constant time.
fn conditional_subtract_l(r: &mut [u32; 9]) {
    let mut diff = [0u32; 9];
    let mut borrow = 0u64;
    for i in 0..9 {
        let d = (r[i] as u64)
            .wrapping_sub(L_LIMBS[i] as u64)
            .wrapping_sub(borrow);
        diff[i] = d as u32;
        borrow = (d >> 63) & 1;
    }
    // `borrow == 0` means r >= L, so the difference is the value we want.
    //
    // Both barriers are load-bearing. Without them LLVM, under LTO, saw that
    // the mask is one of two values, recognised the masked merge below as a
    // select, and compiled it as a jump on the sign of the top limb's
    // difference: a branch on whether the secret nonce or S was at least L.
    // Valgrind found it (docs/CTGRIND.md); `ct_ed25519_reduce_wide` in
    // scripts/check-ct.py now holds the compiled code to having none. The
    // input barrier hides that the borrow is 0 or 1, and the output barrier
    // that the mask is 0 or all ones, as `ic_core::ct::Choice::mask` does.
    let keep = core::hint::black_box(borrow as u32) ^ 1;
    let mask = core::hint::black_box(keep.wrapping_neg());
    for i in 0..9 {
        r[i] ^= mask & (r[i] ^ diff[i]);
    }
}

/// The balanced base-2^21 digits of `-c`, where `L = 2^252 + c`.
///
/// `c` is 125 bits against `L`'s 253, and that is the whole reason this
/// reduction is cheap: `2^252 = L - c`, so `2^252 ≡ -c`, and folding a high
/// limb down multiplies it by something a fifth of the modulus rather than
/// something the size of it.
///
/// Derived, not transcribed: these are the balanced base-2^21 digits of
/// `-(L - 2^252)`, and `the_folding_constants_are_the_digits_of_minus_c`
/// recomputes them from `L` and checks. They agree with the constants ref10
/// publishes, which is the cross-check that the derivation is the right one.
const NEG_C_DIGITS: [i64; 6] = [666_643, 470_296, 654_183, -997_805, 136_657, -683_901];

/// Limbs in the working representation: 21 bits each, covering 512 bits.
const LIMBS: usize = 25;

/// Carry `limbs[i]` into `limbs[i + 1]`, keeping the representation balanced.
#[inline]
fn carry(limbs: &mut [i64; LIMBS], i: usize) {
    let c = (limbs[i] + (1 << 20)) >> 21;
    limbs[i] -= c << 21;
    limbs[i + 1] += c;
}

/// [`carry`] at limbs `0..n`, in order.
///
/// `n` is a constant at every call, and the calls are written out below it
/// rather than looped: as a loop, LLVM kept it rolled on the 32-bit targets,
/// and a rolled loop is a conditional jump the compiled-code gate cannot tell
/// from a secret one.
#[inline(always)]
fn carry_through(limbs: &mut [i64; LIMBS], n: usize) {
    macro_rules! carries {
        ($($k:literal)*) => {$(
            if $k < n {
                carry(limbs, $k);
            }
        )*};
    }
    carries!(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23);
}

/// One round of `reduce_wide`'s folding: carry every limb, then push each
/// limb at or above `2^252` down twelve places, times the digits of `-c`.
#[inline(always)]
fn fold_round(limbs: &mut [i64; LIMBS]) {
    carry_through(limbs, LIMBS - 1);
    // Positions 24 down to 12, written out: as a loop LLVM kept it rolled.
    macro_rules! fold_down {
        ($($i:literal)*) => {$({
            let t = limbs[$i];
            limbs[$i] = 0;
            for (j, d) in NEG_C_DIGITS.iter().enumerate() {
                limbs[$i - 12 + j] += d * t;
            }
        })*};
    }
    fold_down!(24 23 22 21 20 19 18 17 16 15 14 13 12);
}

/// Reduce a 512-bit little-endian integer modulo `L`.
///
/// This used to be long division, one bit at a time: 512 rounds of shifting a
/// nine-limb accumulator and conditionally subtracting `L`, some fifteen
/// thousand operations. Signing does three of these -- two to turn a hash into
/// a scalar and one inside `mul_add` -- and they were about a third of its
/// time.
///
/// Folding instead. A limb at position `21i` for `i >= 12` carries a factor of
/// `2^252`, which is congruent to `-c`, so it can be pushed down twelve places
/// and multiplied by the six digits of `-c`. Thirteen such folds clear
/// everything above `2^252`, and the carries between them keep every limb
/// inside an `i64`.
///
/// Constant time: the loops are fixed, the folding is data-independent, and
/// the final subtractions are masked. The scalar being reduced is secret when
/// signing.
pub fn reduce_wide(input: &[u8; 64]) -> [u8; 32] {
    // Unpack into 21-bit limbs, eight bytes at a time.
    //
    // A limb spans at most four bytes, so one unaligned 64-bit read and a
    // shift produces it. Doing this a bit at a time -- 525 iterations with a
    // bounds test each -- cost more than the folding it feeds.
    //
    // The array is built limb by limb, each index a literal, rather than
    // filled by a loop from a padded copy of the input: on the 32-bit targets
    // the loop stayed rolled and the copy and the zeroing became `memcpy` and
    // `memset` calls, none of which the compiled-code gate can read. Reads
    // past byte 63 are of zero, which is what the padding used to supply.
    #[inline(always)]
    fn limb(input: &[u8; 64], i: usize) -> i64 {
        let bit = i * 21;
        let mut v = 0u64;
        for k in 0..8 {
            let byte = input.get(bit / 8 + k).copied().unwrap_or(0);
            v |= (byte as u64) << (8 * k);
        }
        ((v >> (bit % 8)) & ((1 << 21) - 1)) as i64
    }
    macro_rules! unpack {
        ($($i:literal)*) => { [$(limb(input, $i)),*] };
    }
    let mut limbs: [i64; LIMBS] =
        unpack!(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24);

    // Fold everything at or above 2^252 down, in rounds of carry-then-fold.
    //
    // The two steps have to alternate rather than interleave. A carry pass
    // writes into the limb above, so carrying inside the folding loop puts
    //a value back into a position the downward loop has already passed and will
    // never fold again. Carrying first bounds every limb below 2^21, which is
    // what keeps the products inside an i64: a digit is under 2^20, so a limb
    // receives at most six contributions under 2^41 each.
    //
    // Three rounds suffice -- a fold divides the excess above 2^252 by roughly
    // 2^127, so 2^512 comes down to 2^252 + 2^133 and then to nothing -- and
    // four are run because the count is fixed rather than tested, this being
    // constant-time code. `folding_agrees_with_long_division` covers the
    // largest input there is, which is where too few rounds would show.
    //
    // The four rounds are written out rather than looped. The loop count is
    // public, but LLVM kept the nest rolled, and a compiled loop is a
    // conditional jump that `scripts/check-ct.py` cannot tell from a secret
    // one; written out, `ct_ed25519_reduce_wide` has no jump at all.
    fold_round(&mut limbs);
    fold_round(&mut limbs);
    fold_round(&mut limbs);
    fold_round(&mut limbs);
    // Folding subtracts `c` times the high part, so the residue is congruent
    // to the input but may be negative -- it lies in roughly `(-L, L)`. Two
    // copies of `L` are added to put it safely above zero; the conditional
    // subtractions at the end take them off again.
    //
    // `L = 2^252 + c`, and `2^252` is exactly limb 12, so adding `L` is one
    // increment there and the digits of `c` at the bottom. `c`'s digits are
    // the negation of the folding constants, which is where they come from.
    for _ in 0..2 {
        limbs[12] += 1;
        for (j, d) in NEG_C_DIGITS.iter().enumerate() {
            limbs[j] -= d;
        }
    }
    carry_through(&mut limbs, 13);

    // Now non-negative. Settle into 21-bit limbs.
    let mut u = [0u64; 14];
    let mut borrow: i64 = 0;
    for (slot, &l) in u.iter_mut().zip(limbs.iter().take(14)) {
        let v = l + borrow;
        let m = v & ((1 << 21) - 1);
        borrow = (v - m) >> 21;
        *slot = m as u64;
    }
    debug_assert!(borrow >= 0, "reduction left a negative value");

    // Pack the 21-bit limbs into eight 32-bit words, each from the at most
    // three limbs it overlaps, with every index a literal. This used to merge
    // them into a zeroed byte buffer, which on the 32-bit targets compiled to
    // a `memset` call the compiled-code gate cannot read. Bits from 256 up are
    // dropped, as the buffer's truncation to 32 bytes dropped them: at most a
    // couple of multiples of L remain, well below 2^256.
    #[inline(always)]
    fn word(u: &[u64; 14], w: usize) -> u32 {
        let bit = 32 * w;
        let (i, s) = (bit / 21, bit % 21);
        ((u[i] >> s) | (u[i + 1] << (21 - s)) | (u[i + 2] << (42 - s))) as u32
    }
    macro_rules! pack {
        ($($w:literal)*) => { [$(word(&u, $w)),*, 0] };
    }
    let mut r: [u32; 9] = pack!(0 1 2 3 4 5 6 7);

    // At most a couple of multiples of L remain.
    conditional_subtract_l(&mut r);
    conditional_subtract_l(&mut r);
    conditional_subtract_l(&mut r);
    conditional_subtract_l(&mut r);
    let mut out = [0u8; 32];
    for i in 0..8 {
        out[i * 4..i * 4 + 4].copy_from_slice(&r[i].to_le_bytes());
    }
    out
}

/// Reduce a 256-bit little-endian integer modulo `L`.
#[cfg(test)]
pub fn reduce(input: &[u8; 32]) -> [u8; 32] {
    let mut wide = [0u8; 64];
    wide[..32].copy_from_slice(input);
    reduce_wide(&wide)
}

/// Compute `(a * b + c) mod L` for little-endian 32-byte scalars.
pub fn mul_add(a: &[u8; 32], b: &[u8; 32], c: &[u8; 32]) -> [u8; 32] {
    let al = to_limbs(a);
    let bl = to_limbs(b);

    // Schoolbook 8x8 -> 16 limbs of 32 bits, starting from `c` rather than
    // adding it afterwards. Every limb stays below 2^32 between steps, so
    // each step is at most (2^32 - 1) + (2^32 - 1)^2 + (2^32 - 1) = 2^64 - 1
    // and fits a u64; `c` in the bottom eight limbs is under that bound
    // already.
    //
    // The rows are written out, and `c` enters here, for the reason the
    // folding rounds in `reduce_wide` are written out: the row loop and a
    // byte-wise carry loop that used to add `c` stayed rolled, and a rolled
    // loop is a jump the compiled-code gate cannot tell from a secret one.
    let cl = to_limbs(c);
    let mut prod = [0u64; 16];
    for i in 0..8 {
        prod[i] = cl[i] as u64;
    }
    #[inline(always)]
    fn row(prod: &mut [u64; 16], i: usize, ai: u32, bl: &[u32; 8]) {
        let mut carry = 0u64;
        for j in 0..8 {
            let t = prod[i + j] + (ai as u64) * (bl[j] as u64) + carry;
            prod[i + j] = t & 0xFFFF_FFFF;
            carry = t >> 32;
        }
        prod[i + 8] = carry;
    }
    row(&mut prod, 0, al[0], &bl);
    row(&mut prod, 1, al[1], &bl);
    row(&mut prod, 2, al[2], &bl);
    row(&mut prod, 3, al[3], &bl);
    row(&mut prod, 4, al[4], &bl);
    row(&mut prod, 5, al[5], &bl);
    row(&mut prod, 6, al[6], &bl);
    row(&mut prod, 7, al[7], &bl);

    let mut wide = [0u8; 64];
    for i in 0..16 {
        wide[i * 4..i * 4 + 4].copy_from_slice(&(prod[i] as u32).to_le_bytes());
    }
    reduce_wide(&wide)
}

fn to_limbs(bytes: &[u8; 32]) -> [u32; 8] {
    let mut l = [0u32; 8];
    for i in 0..8 {
        l[i] = u32::from_le_bytes([
            bytes[i * 4],
            bytes[i * 4 + 1],
            bytes[i * 4 + 2],
            bytes[i * 4 + 3],
        ]);
    }
    l
}

/// Constant-time check that a little-endian scalar is canonical, i.e. `< L`.
///
/// RFC 8032 §5.1.7 requires verifiers to reject signatures whose `S` is not
/// canonical; skipping this is what makes an implementation malleable.
#[must_use = "a false return means the scalar encoding was non-canonical"]
pub fn is_canonical(s: &[u8; 32]) -> bool {
    let mut borrow = 0u16;
    for i in 0..32 {
        let d = (s[i] as u16).wrapping_sub(L[i] as u16).wrapping_sub(borrow);
        borrow = (d >> 8) & 1;
    }
    borrow == 1
}

#[cfg(test)]
mod tests {
    use super::*;

    fn from_u64(v: u64) -> [u8; 32] {
        let mut b = [0u8; 32];
        b[..8].copy_from_slice(&v.to_le_bytes());
        b
    }

    #[test]
    fn small_values_are_unchanged() {
        for v in [0u64, 1, 2, 1000, u32::MAX as u64] {
            assert_eq!(reduce(&from_u64(v)), from_u64(v), "{v}");
        }
    }

    #[test]
    fn l_reduces_to_zero() {
        assert_eq!(reduce(&L), [0u8; 32]);
    }

    #[test]
    fn l_plus_one_reduces_to_one() {
        let mut l1 = L;
        // L ends in 0xed, so adding one cannot carry.
        l1[0] += 1;
        assert_eq!(reduce(&l1), from_u64(1));
    }

    #[test]
    fn maximum_wide_value_reduces_into_range() {
        let r = reduce_wide(&[0xffu8; 64]);
        assert!(is_canonical(&r), "reduction must land below L");
    }

    #[test]
    fn mul_add_matches_small_arithmetic() {
        let a = from_u64(1_000_003);
        let b = from_u64(7_919);
        let c = from_u64(65_537);
        let expected = from_u64(1_000_003u64 * 7_919 + 65_537);
        assert_eq!(mul_add(&a, &b, &c), expected);
    }

    #[test]
    fn mul_add_is_zero_for_multiples_of_l() {
        // L * 1 + 0 == 0 (mod L)
        assert_eq!(mul_add(&L, &from_u64(1), &[0u8; 32]), [0u8; 32]);
        // 0 * x + L == 0 (mod L)
        assert_eq!(mul_add(&[0u8; 32], &from_u64(5), &L), [0u8; 32]);
    }

    #[test]
    fn mul_add_result_is_always_canonical() {
        let a = [0xAAu8; 32];
        let b = [0x55u8; 32];
        let c = [0xF0u8; 32];
        assert!(is_canonical(&mul_add(&a, &b, &c)));
    }

    #[test]
    fn canonical_test_matches_the_boundary() {
        assert!(!is_canonical(&L), "L itself is not canonical");
        let mut below = L;
        below[0] -= 1;
        assert!(is_canonical(&below));
        assert!(is_canonical(&[0u8; 32]));
        assert!(!is_canonical(&[0xffu8; 32]));
    }

    /// Reduction must be a ring homomorphism: reducing before or after
    /// multiplying gives the same answer.
    #[test]
    fn reduction_commutes_with_multiplication() {
        let a = [0x37u8; 32];
        let b = [0x91u8; 32];
        let direct = mul_add(&a, &b, &[0u8; 32]);
        let pre = mul_add(&reduce(&a), &reduce(&b), &[0u8; 32]);
        assert_eq!(direct, pre);
    }

    /// The long division the folding replaced, kept as an oracle.
    ///
    /// Transcribed unchanged from the previous revision. It is obviously
    /// correct -- it is school long division with a conditional subtract at
    /// every bit -- and unusably slow, which is the combination an oracle
    /// wants.
    fn reduce_wide_by_long_division(input: &[u8; 64]) -> [u8; 32] {
        let mut r = [0u32; 9];
        for bit_index in (0..512).rev() {
            let mut carry = 0u32;
            for limb in r.iter_mut() {
                let next = *limb >> 31;
                *limb = (*limb << 1) | carry;
                carry = next;
            }
            r[0] |= ((input[bit_index / 8] >> (bit_index % 8)) & 1) as u32;
            conditional_subtract_l(&mut r);
        }
        let mut out = [0u8; 32];
        for i in 0..8 {
            out[i * 4..i * 4 + 4].copy_from_slice(&r[i].to_le_bytes());
        }
        out
    }

    /// The folding constants are the digits of `-c`, recomputed from `L`.
    ///
    /// They are the one part of this that looks like magic numbers, so they
    /// are checked against the definition rather than against a reference.
    #[test]
    fn the_folding_constants_are_the_digits_of_minus_c() {
        // c = L - 2^252, from L's own bytes.
        let mut c = [0u8; 32];
        c.copy_from_slice(&L);
        // Clear bit 252, which is L's leading term.
        c[31] &= !0x10;
        // Balanced base-2^21 digits of -c, low to high.
        let mut v: i128 = 0;
        for (i, b) in c.iter().enumerate().take(16) {
            v |= (*b as i128) << (8 * i);
        }
        let mut v = -v;
        let mut got = [0i64; 6];
        for slot in got.iter_mut() {
            let mut d = (v & ((1 << 21) - 1)) as i64;
            if d >= 1 << 20 {
                d -= 1 << 21;
            }
            *slot = d;
            v = (v - d as i128) >> 21;
        }
        assert_eq!(v, 0, "c did not fit in six digits");
        assert_eq!(got, NEG_C_DIGITS);
    }

    /// The folding reduction against the long division, over random inputs.
    ///
    /// Every bit pattern the folding can be handed, including the ones that
    /// make a limb go negative and borrow: the balanced representation is
    /// where this would go wrong, and it goes wrong on particular inputs
    /// rather than on all of them.
    #[test]
    fn folding_agrees_with_long_division() {
        let mut state = 0x2f6d_9c1b_a473_e850u64;
        let mut next = || {
            state ^= state >> 12;
            state ^= state << 25;
            state ^= state >> 27;
            state.wrapping_mul(0x2545_f491_4f6c_dd1d)
        };
        for case in 0..3_000 {
            let mut input = [0u8; 64];
            match case {
                0 => {}
                1 => input = [0xff; 64],
                2 => input[0] = 1,
                3 => input[63] = 0x80,
                _ => {
                    for chunk in input.chunks_exact_mut(8) {
                        chunk.copy_from_slice(&next().to_le_bytes());
                    }
                    // Sometimes only the low half, so the fold has little to do.
                    if case % 5 == 0 {
                        input[32..].fill(0);
                    }
                }
            }
            assert_eq!(
                reduce_wide(&input),
                reduce_wide_by_long_division(&input),
                "case {case}, input {input:?}"
            );
        }
    }
}