Skip to content

Instantly share code, notes, and snippets.

@resilar
Last active March 19, 2024 10:27
Show Gist options
  • Save resilar/58b46d86675c0a124bb150de18330150 to your computer and use it in GitHub Desktop.
Save resilar/58b46d86675c0a124bb150de18330150 to your computer and use it in GitHub Desktop.
X25519 Elliptic-curve Diffie-Hellman (ECDH) with Curve25519 (RFC7748)
#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;
}
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment