x25519_nostd/

lib.rs

//! Pure-Rust X25519 Elliptic Curve Diffie-Hellman
//!
//! Implements X25519 key exchange (RFC 7748).
//! This is a faithful translation of curve25519-donna-c64 by Adam Langley.
//! Avoids LLVM SIMD issues on x86_64-unknown-none bare-metal targets.
//!
//! Properties:
//! - 128-bit security level
//! - Constant-time Montgomery ladder
//! - 32-byte keys
//! - ~100k operations/sec on modern CPUs
//!
//! Algorithm:
//! - Montgomery curve: y² = x³ + 486662x² + x
//! - Prime: p = 2^255 - 19
//! - Base point: u = 9
//! - Scalar clamping: s[0] &= 248, s[31] &= 127, s[31] |= 64

#![no_std]
#![forbid(unsafe_code)]

/// Field element: 5 limbs of 51 bits each (255 bits total)
/// Represented as u64 to allow for overflow during computation
type Limb = u64;
type FieldElement = [Limb; 5];

/// Load a little-endian 64-bit number
#[inline]
fn load_limb(input: &[u8]) -> Limb {
    (input[0] as Limb) |
    ((input[1] as Limb) << 8) |
    ((input[2] as Limb) << 16) |
    ((input[3] as Limb) << 24) |
    ((input[4] as Limb) << 32) |
    ((input[5] as Limb) << 40) |
    ((input[6] as Limb) << 48) |
    ((input[7] as Limb) << 56)
}

/// Store a 64-bit number as little-endian bytes
#[inline]
fn store_limb(output: &mut [u8], input: Limb) {
    output[0] = (input & 0xff) as u8;
    output[1] = ((input >> 8) & 0xff) as u8;
    output[2] = ((input >> 16) & 0xff) as u8;
    output[3] = ((input >> 24) & 0xff) as u8;
    output[4] = ((input >> 32) & 0xff) as u8;
    output[5] = ((input >> 40) & 0xff) as u8;
    output[6] = ((input >> 48) & 0xff) as u8;
    output[7] = ((input >> 56) & 0xff) as u8;
}

/// Expand a 32-byte little-endian number into field element form (5 limbs of 51 bits)
fn fexpand(output: &mut FieldElement, input: &[u8; 32]) {
    output[0] = load_limb(&input[0..]) & 0x7ffffffffffff;
    output[1] = (load_limb(&input[6..]) >> 3) & 0x7ffffffffffff;
    output[2] = (load_limb(&input[12..]) >> 6) & 0x7ffffffffffff;
    output[3] = (load_limb(&input[19..]) >> 1) & 0x7ffffffffffff;
    output[4] = (load_limb(&input[24..]) >> 12) & 0x7ffffffffffff;
}

/// Contract a field element into a 32-byte little-endian array
fn fcontract(output: &mut [u8; 32], input: &FieldElement) {
    let mut t: [u128; 5] = [
        input[0] as u128,
        input[1] as u128,
        input[2] as u128,
        input[3] as u128,
        input[4] as u128,
    ];

    // First carry pass
    t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
    t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
    t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
    t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
    t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;

    // Second carry pass
    t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
    t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
    t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
    t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
    t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;

    // Now t is between 0 and 2^255-1, properly carried
    // Add 19 to handle case where value is between 2^255-19 and 2^255-1
    t[0] += 19;

    t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
    t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
    t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
    t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
    t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;

    // Now between 19 and 2^255-1, offset by 19
    // Add 2^255 - 19 to each limb to ensure positive and in range
    t[0] += 0x8000000000000u128 - 19;
    t[1] += 0x8000000000000u128 - 1;
    t[2] += 0x8000000000000u128 - 1;
    t[3] += 0x8000000000000u128 - 1;
    t[4] += 0x8000000000000u128 - 1;

    // Now between 2^255 and 2^256-20, offset by 2^255
    t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
    t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
    t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
    t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
    t[4] &= 0x7ffffffffffff;

    // Pack into bytes
    store_limb(&mut output[0..], (t[0] | (t[1] << 51)) as Limb);
    store_limb(&mut output[8..], ((t[1] >> 13) | (t[2] << 38)) as Limb);
    store_limb(&mut output[16..], ((t[2] >> 26) | (t[3] << 25)) as Limb);
    store_limb(&mut output[24..], ((t[3] >> 39) | (t[4] << 12)) as Limb);
}

/// Sum two field elements: output += input
#[inline]
fn fsum(output: &mut FieldElement, input: &FieldElement) {
    output[0] += input[0];
    output[1] += input[1];
    output[2] += input[2];
    output[3] += input[3];
    output[4] += input[4];
}

/// Subtract field elements: output = input - output
/// (Note the order of arguments!)
/// Adds large constants to avoid underflow
#[inline]
fn fdifference_backwards(out: &mut FieldElement, input: &FieldElement) {
    // 152 is 19 << 3
    const TWO54M152: Limb = (1u64 << 54) - 152;
    const TWO54M8: Limb = (1u64 << 54) - 8;

    out[0] = input[0] + TWO54M152 - out[0];
    out[1] = input[1] + TWO54M8 - out[1];
    out[2] = input[2] + TWO54M8 - out[2];
    out[3] = input[3] + TWO54M8 - out[3];
    out[4] = input[4] + TWO54M8 - out[4];
}

/// Multiply field element by scalar
fn fscalar_product(output: &mut FieldElement, input: &FieldElement, scalar: Limb) {
    let mut a: u128;

    a = (input[0] as u128) * (scalar as u128);
    output[0] = (a as Limb) & 0x7ffffffffffff;

    a = (input[1] as u128) * (scalar as u128) + (a >> 51);
    output[1] = (a as Limb) & 0x7ffffffffffff;

    a = (input[2] as u128) * (scalar as u128) + (a >> 51);
    output[2] = (a as Limb) & 0x7ffffffffffff;

    a = (input[3] as u128) * (scalar as u128) + (a >> 51);
    output[3] = (a as Limb) & 0x7ffffffffffff;

    a = (input[4] as u128) * (scalar as u128) + (a >> 51);
    output[4] = (a as Limb) & 0x7ffffffffffff;

    output[0] += ((a >> 51) as Limb) * 19;
}

/// Multiply two field elements
fn fmul(output: &mut FieldElement, in2: &FieldElement, input: &FieldElement) {
    let mut t: [u128; 5] = [0; 5];

    let r0 = input[0];
    let r1 = input[1];
    let r2 = input[2];
    let r3 = input[3];
    let r4 = input[4];

    let s0 = in2[0];
    let s1 = in2[1];
    let s2 = in2[2];
    let s3 = in2[3];
    let s4 = in2[4];

    t[0] = (r0 as u128) * (s0 as u128);
    t[1] = (r0 as u128) * (s1 as u128) + (r1 as u128) * (s0 as u128);
    t[2] = (r0 as u128) * (s2 as u128) + (r2 as u128) * (s0 as u128) + (r1 as u128) * (s1 as u128);
    t[3] = (r0 as u128) * (s3 as u128) + (r3 as u128) * (s0 as u128) + (r1 as u128) * (s2 as u128) + (r2 as u128) * (s1 as u128);
    t[4] = (r0 as u128) * (s4 as u128) + (r4 as u128) * (s0 as u128) + (r3 as u128) * (s1 as u128) + (r1 as u128) * (s3 as u128) + (r2 as u128) * (s2 as u128);

    let r4 = r4 * 19;
    let r1 = r1 * 19;
    let r2 = r2 * 19;
    let r3 = r3 * 19;

    t[0] += (r4 as u128) * (s1 as u128) + (r1 as u128) * (s4 as u128) + (r2 as u128) * (s3 as u128) + (r3 as u128) * (s2 as u128);
    t[1] += (r4 as u128) * (s2 as u128) + (r2 as u128) * (s4 as u128) + (r3 as u128) * (s3 as u128);
    t[2] += (r4 as u128) * (s3 as u128) + (r3 as u128) * (s4 as u128);
    t[3] += (r4 as u128) * (s4 as u128);

    let mut r0: Limb;
    let mut r1: Limb;
    let mut r2: Limb;
    let r3: Limb;
    let r4: Limb;
    let mut c: Limb;

    r0 = (t[0] as Limb) & 0x7ffffffffffff; c = (t[0] >> 51) as Limb;
    t[1] += c as u128; r1 = (t[1] as Limb) & 0x7ffffffffffff; c = (t[1] >> 51) as Limb;
    t[2] += c as u128; r2 = (t[2] as Limb) & 0x7ffffffffffff; c = (t[2] >> 51) as Limb;
    t[3] += c as u128; r3 = (t[3] as Limb) & 0x7ffffffffffff; c = (t[3] >> 51) as Limb;
    t[4] += c as u128; r4 = (t[4] as Limb) & 0x7ffffffffffff; c = (t[4] >> 51) as Limb;
    r0 += c * 19; c = r0 >> 51; r0 &= 0x7ffffffffffff;
    r1 += c; c = r1 >> 51; r1 &= 0x7ffffffffffff;
    r2 += c;

    output[0] = r0;
    output[1] = r1;
    output[2] = r2;
    output[3] = r3;
    output[4] = r4;
}

/// Square a field element (multiple times)
fn fsquare_times(output: &mut FieldElement, input: &FieldElement, count: usize) {
    let mut t: [u128; 5] = [0; 5];
    let mut r0 = input[0];
    let mut r1 = input[1];
    let mut r2 = input[2];
    let mut r3 = input[3];
    let mut r4 = input[4];

    for _ in 0..count {
        let d0 = r0 * 2;
        let d1 = r1 * 2;
        let d2 = r2 * 2 * 19;
        let d419 = r4 * 19;
        let d4 = d419 * 2;

        t[0] = (r0 as u128) * (r0 as u128) + (d4 as u128) * (r1 as u128) + (d2 as u128) * (r3 as u128);
        t[1] = (d0 as u128) * (r1 as u128) + (d4 as u128) * (r2 as u128) + (r3 as u128) * ((r3 * 19) as u128);
        t[2] = (d0 as u128) * (r2 as u128) + (r1 as u128) * (r1 as u128) + (d4 as u128) * (r3 as u128);
        t[3] = (d0 as u128) * (r3 as u128) + (d1 as u128) * (r2 as u128) + (r4 as u128) * (d419 as u128);
        t[4] = (d0 as u128) * (r4 as u128) + (d1 as u128) * (r3 as u128) + (r2 as u128) * (r2 as u128);

        let mut c: Limb;
        r0 = (t[0] as Limb) & 0x7ffffffffffff; c = (t[0] >> 51) as Limb;
        t[1] += c as u128; r1 = (t[1] as Limb) & 0x7ffffffffffff; c = (t[1] >> 51) as Limb;
        t[2] += c as u128; r2 = (t[2] as Limb) & 0x7ffffffffffff; c = (t[2] >> 51) as Limb;
        t[3] += c as u128; r3 = (t[3] as Limb) & 0x7ffffffffffff; c = (t[3] >> 51) as Limb;
        t[4] += c as u128; r4 = (t[4] as Limb) & 0x7ffffffffffff; c = (t[4] >> 51) as Limb;
        r0 += c * 19; c = r0 >> 51; r0 &= 0x7ffffffffffff;
        r1 += c; c = r1 >> 51; r1 &= 0x7ffffffffffff;
        r2 += c;
    }

    output[0] = r0;
    output[1] = r1;
    output[2] = r2;
    output[3] = r3;
    output[4] = r4;
}

/// Compute multiplicative inverse using Fermat's little theorem
/// out = z^(p-2) = z^(2^255 - 21)
fn crecip(out: &mut FieldElement, z: &FieldElement) {
    let mut a = [0u64; 5];
    let mut t0 = [0u64; 5];
    let mut b = [0u64; 5];
    let mut c = [0u64; 5];

    // 2
    fsquare_times(&mut a, z, 1);
    // 8
    fsquare_times(&mut t0, &a, 2);
    // 9
    fmul(&mut b, &t0, z);
    // 11
    let a_copy = a;
    fmul(&mut a, &b, &a_copy);
    // 22
    fsquare_times(&mut t0, &a, 1);
    // 2^5 - 2^0 = 31
    let b_copy = b;
    fmul(&mut b, &t0, &b_copy);
    // 2^10 - 2^5
    fsquare_times(&mut t0, &b, 5);
    // 2^10 - 2^0
    let b_copy = b;
    fmul(&mut b, &t0, &b_copy);
    // 2^20 - 2^10
    fsquare_times(&mut t0, &b, 10);
    // 2^20 - 2^0
    fmul(&mut c, &t0, &b);
    // 2^40 - 2^20
    fsquare_times(&mut t0, &c, 20);
    // 2^40 - 2^0
    let t0_copy = t0;
    fmul(&mut t0, &t0_copy, &c);
    // 2^50 - 2^10
    let t0_copy = t0;
    fsquare_times(&mut t0, &t0_copy, 10);
    // 2^50 - 2^0
    let b_copy = b;
    fmul(&mut b, &t0, &b_copy);
    // 2^100 - 2^50
    fsquare_times(&mut t0, &b, 50);
    // 2^100 - 2^0
    fmul(&mut c, &t0, &b);
    // 2^200 - 2^100
    fsquare_times(&mut t0, &c, 100);
    // 2^200 - 2^0
    let t0_copy = t0;
    fmul(&mut t0, &t0_copy, &c);
    // 2^250 - 2^50
    let t0_copy = t0;
    fsquare_times(&mut t0, &t0_copy, 50);
    // 2^250 - 2^0
    let t0_copy = t0;
    fmul(&mut t0, &t0_copy, &b);
    // 2^255 - 2^5
    let t0_copy = t0;
    fsquare_times(&mut t0, &t0_copy, 5);
    // 2^255 - 21
    fmul(out, &t0, &a);
}

/// Constant-time conditional swap
fn swap_conditional(a: &mut FieldElement, b: &mut FieldElement, iswap: Limb) {
    let swap = iswap.wrapping_neg(); // 0 or 0xFFFFFFFFFFFFFFFF

    for i in 0..5 {
        let x = swap & (a[i] ^ b[i]);
        a[i] ^= x;
        b[i] ^= x;
    }
}

/// Montgomery ladder step
#[allow(clippy::too_many_arguments)]
fn fmonty(
    x2: &mut FieldElement, z2: &mut FieldElement,
    x3: &mut FieldElement, z3: &mut FieldElement,
    x: &mut FieldElement, z: &mut FieldElement,
    xprime: &mut FieldElement, zprime: &mut FieldElement,
    qmqp: &FieldElement,
) {
    let mut origx = [0u64; 5];
    let mut origxprime = [0u64; 5];
    let mut zzz = [0u64; 5];
    let mut xx = [0u64; 5];
    let mut zz = [0u64; 5];
    let mut xxprime = [0u64; 5];
    let mut zzprime = [0u64; 5];
    let mut zzzprime = [0u64; 5];

    origx.copy_from_slice(x);
    fsum(x, z);
    fdifference_backwards(z, &origx);

    origxprime.copy_from_slice(xprime);
    fsum(xprime, zprime);
    fdifference_backwards(zprime, &origxprime);
    fmul(&mut xxprime, xprime, z);
    fmul(&mut zzprime, x, zprime);
    origxprime.copy_from_slice(&xxprime);
    fsum(&mut xxprime, &zzprime);
    fdifference_backwards(&mut zzprime, &origxprime);
    fsquare_times(x3, &xxprime, 1);
    fsquare_times(&mut zzzprime, &zzprime, 1);
    fmul(z3, &zzzprime, qmqp);

    fsquare_times(&mut xx, x, 1);
    fsquare_times(&mut zz, z, 1);
    fmul(x2, &xx, &zz);
    fdifference_backwards(&mut zz, &xx);
    fscalar_product(&mut zzz, &zz, 121665);
    fsum(&mut zzz, &xx);
    fmul(z2, &zz, &zzz);
}

/// Scalar multiplication: resultx/resultz = n * q
fn cmult(resultx: &mut FieldElement, resultz: &mut FieldElement, n: &[u8; 32], q: &FieldElement) {
    let mut a = [0u64; 5];
    let mut b = [1u64, 0, 0, 0, 0];
    let mut c = [1u64, 0, 0, 0, 0];
    let mut d = [0u64; 5];
    let mut e = [0u64; 5];
    let mut f = [1u64, 0, 0, 0, 0];
    let mut g = [0u64; 5];
    let mut h = [1u64, 0, 0, 0, 0];

    let nqpqx = &mut a;
    let nqpqz = &mut b;
    let nqx = &mut c;
    let nqz = &mut d;
    let nqpqx2 = &mut e;
    let nqpqz2 = &mut f;
    let nqx2 = &mut g;
    let nqz2 = &mut h;

    nqpqx.copy_from_slice(q);

    for i in 0..32 {
        let mut byte = n[31 - i];
        for _ in 0..8 {
            let bit = (byte >> 7) as Limb;

            swap_conditional(nqx, nqpqx, bit);
            swap_conditional(nqz, nqpqz, bit);
            fmonty(nqx2, nqz2, nqpqx2, nqpqz2, nqx, nqz, nqpqx, nqpqz, q);
            swap_conditional(nqx2, nqpqx2, bit);
            swap_conditional(nqz2, nqpqz2, bit);

            core::mem::swap(&mut *nqx, &mut *nqx2);
            core::mem::swap(&mut *nqz, &mut *nqz2);
            core::mem::swap(&mut *nqpqx, &mut *nqpqx2);
            core::mem::swap(&mut *nqpqz, &mut *nqpqz2);

            byte <<= 1;
        }
    }

    resultx.copy_from_slice(nqx);
    resultz.copy_from_slice(nqz);
}

/// X25519 scalar multiplication
fn curve25519_donna(mypublic: &mut [u8; 32], secret: &[u8; 32], basepoint: &[u8; 32]) {
    let mut bp = [0u64; 5];
    let mut x = [0u64; 5];
    let mut z = [0u64; 5];
    let mut zmone = [0u64; 5];
    let mut e = *secret;

    // Clamp scalar
    e[0] &= 248;
    e[31] &= 127;
    e[31] |= 64;

    fexpand(&mut bp, basepoint);
    cmult(&mut x, &mut z, &e, &bp);
    crecip(&mut zmone, &z);
    fmul(&mut z, &x, &zmone);
    fcontract(mypublic, &z);
}

/// Base point for X25519 (u=9)
const BASEPOINT: [u8; 32] = [
    9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
];

/// Compute public key from secret key
pub fn public_key(secret: &[u8; 32]) -> [u8; 32] {
    let mut public = [0u8; 32];
    curve25519_donna(&mut public, secret, &BASEPOINT);
    public
}

/// Compute shared secret from our secret key and their public key
pub fn diffie_hellman(our_secret: &[u8; 32], their_public: &[u8; 32]) -> [u8; 32] {
    let mut shared = [0u8; 32];
    curve25519_donna(&mut shared, our_secret, their_public);
    shared
}

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

    #[test]
    fn test_rfc7748_vector1() {
        // RFC 7748 Section 5.2 Test Vector 1
        let scalar = [
            0xa5, 0x46, 0xe3, 0x6b, 0xf0, 0x52, 0x7c, 0x9d,
            0x3b, 0x16, 0x15, 0x4b, 0x82, 0x46, 0x5e, 0xdd,
            0x62, 0x14, 0x4c, 0x0a, 0xc1, 0xfc, 0x5a, 0x18,
            0x50, 0x6a, 0x22, 0x44, 0xba, 0x44, 0x9a, 0xc4,
        ];

        let basepoint = [
            0xe6, 0xdb, 0x68, 0x67, 0x58, 0x30, 0x30, 0xdb,
            0x35, 0x94, 0xc1, 0xa4, 0x24, 0xb1, 0x5f, 0x7c,
            0x72, 0x66, 0x24, 0xec, 0x26, 0xb3, 0x35, 0x3b,
            0x10, 0xa9, 0x03, 0xa6, 0xd0, 0xab, 0x1c, 0x4c,
        ];

        let expected = [
            0xc3, 0xda, 0x55, 0x37, 0x9d, 0xe9, 0xc6, 0x90,
            0x8e, 0x94, 0xea, 0x4d, 0xf2, 0x8d, 0x08, 0x4f,
            0x32, 0xec, 0xcf, 0x03, 0x49, 0x1c, 0x71, 0xf7,
            0x54, 0xb4, 0x07, 0x55, 0x77, 0xa2, 0x85, 0x52,
        ];

        let result = diffie_hellman(&scalar, &basepoint);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_basepoint_mult() {
        let secret = [
            1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        ];

        let public = public_key(&secret);
        assert_ne!(public, [0u8; 32]);
    }

    #[test]
    fn test_dh_symmetry() {
        let alice_secret = [0x77u8; 32];
        let alice_public = public_key(&alice_secret);

        let bob_secret = [0x99u8; 32];
        let bob_public = public_key(&bob_secret);

        let alice_shared = diffie_hellman(&alice_secret, &bob_public);
        let bob_shared = diffie_hellman(&bob_secret, &alice_public);

        assert_eq!(alice_shared, bob_shared);
    }

    #[test]
    fn test_public_key_deterministic() {
        let secret = [0x42u8; 32];
        let public1 = public_key(&secret);
        let public2 = public_key(&secret);
        assert_eq!(public1, public2);
    }
}