X25519 Elliptic-curve Diffie-Hellman (ECDH) with Curve25519 (RFC7748)

x25519.c

#include <stdint.h>
#include <sys/random.h>

int x25519_keygen(uint8_t secret[32], uint8_t public[32]);
void x25519_secret_to_public(const uint8_t secret[32], uint8_t public[32]);
int x25519(const uint8_t secret[32], const uint8_t public[32], uint8_t out[32]);

/*
 * RFC7748 X25519 Elliptic Curve Diffie-Hellman (ECDH) with Curve25519.
 * Montgomery curve y² = x³ + 486662x² + x over GF(p) = GF(2^255 - 19).
 *
 * 256-bit arithmetic modulo 2^255 - 19 implemented with 8 x 32-bit limbs :3
 */

static uint32_t *import(const uint8_t in[32], uint32_t out[8])
{
    out[0] = in[ 0] + (in[ 1] << 8) + (in[ 2] << 16) + (in[ 3] << 24);
    out[1] = in[ 4] + (in[ 5] << 8) + (in[ 6] << 16) + (in[ 7] << 24);
    out[2] = in[ 8] + (in[ 9] << 8) + (in[10] << 16) + (in[11] << 24);
    out[3] = in[12] + (in[13] << 8) + (in[14] << 16) + (in[15] << 24);
    out[4] = in[16] + (in[17] << 8) + (in[18] << 16) + (in[19] << 24);
    out[5] = in[20] + (in[21] << 8) + (in[22] << 16) + (in[23] << 24);
    out[6] = in[24] + (in[25] << 8) + (in[26] << 16) + (in[27] << 24);
    out[7] = in[28] + (in[29] << 8) + (in[30] << 16) + (in[31] << 24);
    return out;
}

static uint8_t *export(const uint32_t i[8], uint8_t o[32])
{
    o[ 0] = i[0]; o[ 1] = i[0] >> 8; o[ 2] = i[0] >> 16; o[ 3] = i[0] >> 24;
    o[ 4] = i[1]; o[ 5] = i[1] >> 8; o[ 6] = i[1] >> 16; o[ 7] = i[1] >> 24;
    o[ 8] = i[2]; o[ 9] = i[2] >> 8; o[10] = i[2] >> 16; o[11] = i[2] >> 24;
    o[12] = i[3]; o[13] = i[3] >> 8; o[14] = i[3] >> 16; o[15] = i[3] >> 24;
    o[16] = i[4]; o[17] = i[4] >> 8; o[18] = i[4] >> 16; o[19] = i[4] >> 24;
    o[20] = i[5]; o[21] = i[5] >> 8; o[22] = i[5] >> 16; o[23] = i[5] >> 24;
    o[24] = i[6]; o[25] = i[6] >> 8; o[26] = i[6] >> 16; o[27] = i[6] >> 24;
    o[28] = i[7]; o[29] = i[7] >> 8; o[30] = i[7] >> 16; o[31] = i[7] >> 24;
    return o;
}

static uint32_t *clamp_scalar(uint32_t a[8])
{
    a[0] &= 0xfffffff8;
    a[7] &= 0x7fffffff;
    a[7] |= 0x40000000;
    return a;
}

static uint32_t *clamp_coordinate(uint32_t u[8])
{
    u[7] &= 0x7fffffff;
    return u;
}

/* x = a (mod p) */
static uint32_t *modp(const uint32_t a[8], uint32_t x[8])
{
    int i;
    uint64_t c = (uint64_t)a[0] + 19;
    for (i = 1; i < 8; i++) c = (c >> 32) + a[i];
    c = (c >> 31) * 19;
    for (i = 0; i < 7; i++) {
        c += a[i];
        x[i] = c;
        c >>= 32;
    }
    x[i] = ((uint32_t)c + a[i]) & 0x7fffffff;
    return x;
}

/* x = a + b (mod p) */
static uint32_t *add(const uint32_t a[8], const uint32_t b[8], uint32_t x[8])
{
    int i;
    uint64_t c = 0;
    for (i = 0; i < 8; i++) {
        c += a[i];
        c += b[i];
        x[i] = c;
        c >>= 32;
    }
    return x;
    /*
     * Note that we can skip modp() here because we do not chain more than one
     * addition or subtraction between multiplications. mul() handles unreduced
     * integers just fine and returns (mod p)-reduced product as the result.
     */
}

/* x = a - b (mod p) */
static uint32_t *sub(const uint32_t a[8], const uint32_t b[8], uint32_t x[8])
{
    int i;
    int64_t c = -19;
    for (i = 0; i < 8; i++) {
        c += a[i];
        c -= b[i];
        x[i] = c;
        c >>= 32;
    }
    x[7] += 0x80000000;
    return x; /* Again, skip modp() */
}

/* x = a * b (mod p) */
static uint32_t *mul(const uint32_t a[8], const uint32_t b[8], uint32_t x[8])
{
    int i, j;
    uint32_t ab[8];
    uint64_t c, d, l, r;
    for (i = c = 0; i < 8; i++) {
        l = c & 0xffffffff;
        r = 0;
        c >>= 32;
        for (j = 0; j <= i; j++) {
            l += (uint64_t)a[j] * b[i - j];
            c += l >> 32;
            l &= 0xffffffff;
        }
        for (d = 0; j < 8; j++) {
            r += (uint64_t)a[j] * b[8 - (j - i)];
            d += r >> 32;
            r &= 0xffffffff;
        }
        l += 38 * r;
        c += 38 * d + (l >> 32);
        ab[i] = l;
    }
    c = 19 * (2*c + (ab[7] >> 31));
    ab[7] &= 0x7fffffff;
    for (i = 0; i < 8; i++) {
        c += ab[i];
        ab[i] = c;
        c >>= 32;
    }
    return modp(ab, x);
}

/* x = a² (mod p) */
static uint32_t *square(const uint32_t a[8], uint32_t x[8])
{
    return mul(a, a, x);
}

/* x = a^(p-2) (mod p) */
static uint32_t *inverse(const uint32_t a[8], uint32_t x[8])
{
    int i;
    uint32_t b[8];
    mul(a, square(a, b), x);
    mul(x, square(square(b, b), b), x);
    square(b, b);
    for (i = 5; i < 255; i++) mul(x, square(b, b), x);
    return x;
}

/* Swap a and b if c != 0 */
static void cswap(uint32_t a[8], uint32_t b[8], int c)
{
    int i;
    uint32_t mask = -(uint32_t)!!c;
    for (i = 0; i < 8; i++) {
        const uint32_t xor = mask & (a[i] ^ b[i]);
        a[i] ^= xor;
        b[i] ^= xor;
    }
}

/*
 * Constant-time Montgomery ladder for elliptic curve point multiplication with
 * a scalar k and point coordinate u. We expect k and u to be properly clamped.
 */
static uint32_t *monty(const uint32_t k[8], const uint32_t u[8], uint32_t x[8])
{
    int i, swap = 0;
    uint32_t x0[8] = { 1, 0, 0, 0, 0, 0, 0, 0 };
    uint32_t z0[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
    uint32_t x1[8];
    uint32_t z1[8] = { 1, 0, 0, 0, 0, 0, 0, 0 };
    const uint32_t a24[8] = { (486662 - 2) / 4, 0, 0, 0, 0, 0, 0, 0 };
    for (i = 0; i < 8; i++) x1[i] = u[i];
    for (i = 1; i <= 255; i++) {
        uint32_t A[8], AA[8], B[8], BB[8], C[8], CB[8], D[8], DA[8],E[8];
        const int bit = (k[(255-i) / 32] >> ((255-i) % 32)) & 1;
        swap ^= bit;
        cswap(x0, x1, swap);
        cswap(z0, z1, swap);
        swap = bit;

        /*
         * A = x_0 + z_0, AA = A^2
         * B = x_0 - z_0, BB = B^2
         * C = x_1 + z_1, CB = C * B
         * D = x_1 - z_1, DA = D * A
         * E = AA - BB
         */
        square(add(x0, z0, A), AA);
        square(sub(x0, z0, B), BB);
        mul(add(x1, z1, C), B, CB);
        mul(sub(x1, z1, D), A, DA);
        sub(AA, BB, E);

        /*
         * x_0 = AA * BB
         * x_1 = (DA + CB)^2
         */
        mul(AA, BB, x0);
        square(add(DA, CB, x1), x1);

        /*
         * z_0 = E * (AA + a24 * E) where a24 = (486662 - 2) / 4 = 121665
         * z_1 = (DA - CB)^2 * u
         */
        add(AA, mul(a24, E, z0), z0);
        mul(E, z0, z0);
        square(sub(DA, CB, z1), z1);
        mul(u, z1, z1);
    }

    cswap(x0, x1, swap);
    cswap(z0, z1, swap);
    return mul(x0, inverse(z0, x), x);
}

/*
 * Generate a public/private key pair for x25519 key exchange.
 *
 * The secret key is generated with getrandom(2) and clamped/masked according
 * to RFC7748. The public key is derived from the secret key by elliptic curve
 * point multiplication with the generator point G = { 9 } i.e. the Montgomery
 * curve point u=9.
 */
int x25519_keygen(uint8_t secret[32], uint8_t public[32])
{
    uint32_t k[8];
    const uint32_t G[8] = { 9, 0, 0, 0, 0, 0, 0, 0 };
    if (getrandom(k, sizeof(k), 0) != sizeof(k))
        return 0;
    clamp_scalar(k);
    export(k, secret);
    monty(k, G, k);
    export(k, public);
    return 1;
}

/*
 * Compute a public key from secret.
 */
void x25519_secret_to_public(const uint8_t secret[32], uint8_t public[32])
{
    uint32_t k[8];
    const uint32_t G[8] = { 9, 0, 0, 0, 0, 0, 0, 0 };
    import(secret, k);
    clamp_scalar(k);
    monty(k, G, k);
    export(k, public);
}

/*
 * X25519 Elliptic-curve Diffie-Hellman (ECDH).
 *
 * x25519(k, u) <=> U(kP) = x25519(k, U(P)) where P = (u, v) and U(P) = u.
 *
 *  K_a = clamp_scalar(urandom(32))             Alice's private key
 *  P_a = x25519(K_a, G) = K_a x G = K_a x 9    Alice's public key
 *  K_b = clamp_scalar(urandom(32))             Bob's private key
 *  P_b = x25519(K_b, G) = K_b x G = K_b x 9    Bob's public key
 *
 * Elliptic-curve Diffie-Hellman produces a shared key for Alice & Bob:
 *
 *  K = x25519(K_a, P_b) = x25519(K_b, P_a)     Alice's & Bob's shared key
 *
 * Proof of correctness:
 * (A) x25519(K_a, P_b) = x25519(K_a, x25519(K_b, 9)) = K_a x K_b x 9 = K
 * (B) x25519(K_b, P_a) = x25519(K_b, x25519(K_a, 9)) = K_b x K_a x 9 = K
 *  => Alice and Bob derive identical shared key K.
 *
 * The return value is non-zero on success. An elliptic curve point with small
 * order causes the function return to zero and out buffer receives only zeros.
 */
int x25519(const uint8_t secret[32], const uint8_t public[32], uint8_t out[32])
{
    uint32_t k[8], u[8];
    import(secret, k);
    import(public, u);
    clamp_scalar(k);
    clamp_coordinate(u);
    monty(k, u, k);
    export(k, out);
    return (k[0] | k[1] | k[2] | k[3] | k[4] | k[5] | k[6] | k[7]) != 0;
}

x25519_cli.c

#include <stdio.h>
#include <string.h>

#include "x25519.c"

void usage(const char *argv0)
{
    fprintf(stderr,
        "Key generation: %s -g\n"
        "Key derivation: %s [secret_key]\n"
        "Key exchange:   %s [secret_key] [public_key]\n",
        argv0, argv0, argv0
    );
}

int curve25519_atoi(char *a, uint8_t i[32])
{
    uint32_t ten[8] = { 10, 0, 0, 0, 0, 0, 0, 0 };
    uint32_t x[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
    uint32_t y[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
    while (*a) {
        if (*a < '0' || *a > '9')
            return 0;

        y[0] = *a - '0';
        mul(x, ten, x);
        add(x, y, x);
        a++;
    }
    export(x, i);
    return 1;
}

char *curve25519_itoa_to(const uint32_t x[8], char out[80])
{
    int i, j, c;
    unsigned char res[80] = {0};
    unsigned char acc[80] = {1};
    const unsigned int w = sizeof(*x) * 8;

    for (i = 0; i < 256; i++) {
        if (x[i / w] & ((uint32_t)1 << (i % w))) {
            for (j = c = 0; j < 80; j++) {
                res[j] += acc[j] + c;
                if ((c = res[j] >= 10))
                    res[j] -= 10;
            }
        }
        for (j = c = 0; j < 80; j++) {
            acc[j] += acc[j] + c;
            if ((c = acc[j] >= 10))
                acc[j] -= 10;
        }
    }

    for (i = j = 0; i < 80; i++) if (res[i] > 0) j = i;
    for (i = 0; i <= j; i++) out[i] = res[j-i] + '0';
    out[i] = '\0';
    return out;
}

char *curve25519_itoa(const uint8_t x[32])
{
    static char buf[80];
    return curve25519_itoa_to((uint32_t *)x, buf);
}

char *curve25519_hexify_to(const uint32_t x[8], char out[80])
{
    int i;
    char *p = out;
    for (i = 0; i < 32; i++) {
        const uint8_t byte = x[i / 4] >> 8*(i % 4);
        const uint8_t lo = byte & 0x0f;
        const uint8_t hi = (byte & 0xf0) >> 4;
        *p++ = hi < 10 ? '0' + hi : 'a' + (hi - 10);
        *p++ = lo < 10 ? '0' + lo : 'a' + (lo - 10);
    }
    for  (i *= 2; i < 80; *p++ = '\0', i++);
    return out;
}

char *curve25519_hexify(const uint32_t x[8])
{
    static char buf[80];
    return curve25519_hexify_to(x, buf);
}

uint32_t *curve25519_unhexify(const char *in, uint32_t x[8])
{
    int i; 
    if (*in == '0' && (in[1] == 'x' || in[1] == 'X')) {
        in += 2;
    }
    for (i = 0; i < 8; i++) x[i] = 0;
    for (i = 0; i < 256 && *in; i += 4) {
        const char c = *in++;
        const int j = i / 32;
        const int s = (i ^ 4) % 32;
        if ('0' <= c && c <= '9') {
            x[j] |= (c - '0') << s;
        } else if ('a' <= c && c <= 'f') {
            x[j] |= (c - 'a' + 10) << s;
        } else if ('A' <= c && c <= 'F') {
            x[j] |= (c - 'A' + 10) << s;
        }
    }
    return x;
}

int main(int argc, char *argv[])
{
    uint8_t pk[32], sk[32], k[32];

    if (argc <= 1 || (argc == 2 && !strcmp(argv[1], "-h"))) {
        usage((argc >= 1 && argv[0]) ? argv[0] : "./x25519");
    } else if (argc == 2 && !strcmp(argv[1], "-g")) {
        if (!x25519_keygen(sk, pk)) {
            fprintf(stderr, "x25519_keygen: getrandom() error\n");
            return 2;
        }
        fprintf(stderr, "secret key: ");
        fflush(stderr);
        printf("%s\n", curve25519_itoa(sk));
        fprintf(stderr, "public key: ");
        fflush(stderr);
        printf("%s\n", curve25519_itoa(pk));
    } else if (argc == 2 && curve25519_atoi(argv[1], sk)) {
        uint8_t G[32] = { 9 };
        x25519(sk, G, pk);
        fprintf(stderr, "secret key: ");
        fflush(stderr);
        printf("%s\n", curve25519_itoa(sk));
        fprintf(stderr, "public key: ");
        fflush(stderr);
        printf("%s\n", curve25519_itoa(pk));
    } else if (argc == 3
            && curve25519_atoi(argv[1], sk)
            && curve25519_atoi(argv[2], pk))
    {
        x25519(sk, pk, k);
        printf("%s\n", curve25519_itoa(k));
    } else {
        usage((argc >= 1 && argv[0]) ? argv[0] : "./x25519");
        return 1;
    }

    return 0;
}