ECDH using arbitrary-precision arithmetic

bn_ecc.c

/**
    Running test...key : 803D8AB2E5B6E6FCA715737C3A82F7CE3C783124F6D51CD0
    Session keys match...
    OK.
    Generating random keys...OK
    Private Key for Alice : 31FA1084
    Private Key for Bob   : 2D748885

    Generating public keys...
    Public Key for Alice : (10DEE015, 7A458BE8)
    Public Key for Bob   : (49DC9B2E, 660EDB3F)

    Generating session keys...
    Session Key for Alice : (41B46B35, 45594721)
    Session Key for Bob   : (41B46B35, 45594721)

    Shared secret: 41B46B35
*/
#if defined(_WIN32) || defined(_WIN64)
#define WINDOWS
#include <windows.h>
#include <wincrypt.h>
#pragma comment(lib, "advapi32")
#else
#define NIX
#include <sys/types.h>
#include <unistd.h>
#include <fcntl.h>
#include <errno.h>
#endif

#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <ctype.h>
#include <time.h>

#include "sha256.h"

typedef union w64_t {
    uint8_t  b[4];
    uint16_t h[2];
    uint32_t w;
} w64_t;

#if defined(NIX)

/* ---------- random helpers (POSIX) ------------------------------------ */
uint32_t random32(void) {
    int fd;
    ssize_t len, u = 0;
    const size_t outlen = 4;
    w64_t rnd = {0};

    fd = open("/dev/urandom", O_RDONLY);
    if (fd >= 0) {
        for (u = 0; u < (ssize_t)outlen;) {
            len = read(fd, rnd.b + u, outlen - u);
            if (len < 0) break;
            u += len;
        }
        close(fd);
    }
    return rnd.w;
}

/* renamed to avoid collision with libc random()               */
/*                                                             */
/* NOTE: the old name ‘random()’ is kept as an alias via macro */
/* so existing code (if any) still compiles.                   */
static void rand_bytes(void *out, size_t outlen) {            
    int fd;
    ssize_t len, u = 0;
    uint8_t *ptr = (uint8_t*)out;

    fd = open("/dev/urandom", O_RDONLY);
    if (fd >= 0) {
        for (u = 0; u < (ssize_t)outlen;) {
            len = read(fd, ptr + u, outlen - u);
            if (len < 0) break;
            u += len;
        }
        close(fd);
    }
}
#define random  rand_bytes    //  keep compatibility if older code used it

#else   /* -------------------- Windows version ------------------------- */

uint32_t random32(void) {
    HCRYPTPROV hp;
    BOOL ok = FALSE;
    const DWORD outlen = 4;
    w64_t rnd = {0};

    if (CryptAcquireContext(&hp, 0, 0,
            PROV_RSA_FULL, CRYPT_VERIFYCONTEXT | CRYPT_SILENT)) {
        ok = CryptGenRandom(hp, outlen, (PBYTE)&rnd);
        CryptReleaseContext(hp, 0);
    }
    (void)ok;                         /* suppress unused-value warning   */
    return rnd.w;
}

static void rand_bytes(void *out, size_t outlen) {            
    HCRYPTPROV hp;
    BOOL ok = FALSE;

    if (CryptAcquireContext(&hp, 0, 0,
            PROV_RSA_FULL, CRYPT_VERIFYCONTEXT | CRYPT_SILENT)) {
        ok = CryptGenRandom(hp, (DWORD)outlen, (PBYTE)out);
        CryptReleaseContext(hp, 0);
    }
    (void)ok;
}
#define random  rand_bytes

#endif  /* ---------------------------------------------------------------- */

#define BN_MAX_BITS  384
#define BN_MAX_BYTES (BN_MAX_BITS / 8)
#define BN_MAX_WORDS (BN_MAX_BYTES / sizeof(uint32_t))

typedef union _bn_t {
    uint8_t  b[BN_MAX_BYTES];
    uint32_t w[BN_MAX_BYTES / sizeof(uint32_t)];
    uint64_t q[BN_MAX_BYTES / sizeof(uint64_t)];
} bn;

/* ------ bignum forward decls (from earlier source) -------------------- */
void  bn_set_word(bn *x, uint32_t w);
const char* bn_bn2hex(bn *a);
void  bn_modexp(bn *x, bn *b, bn *e, bn *m);
void  bn_mulmod(bn *x, bn *b, bn *e, bn *m);
void  bn_addmod(bn *r, bn *a, bn *b, bn *m);
void  bn_submod(bn *r, bn *a, bn *b, bn *m);
void  bn_mod(bn *dividend, bn *divisor, bn *remainder);
int   bn_is_zero(bn*);
int   bn_cmp(bn *a, bn *b);
void  bn_copy(bn *d, bn *s);
void  bn_zero(bn *x);
int   bn_is_bit_set(bn *x, int bit_idx);
int   bn_hex2bn(bn *a, const char *hexStr);
int   bn_num_bits(bn *x);
const char* bn_bn2dec(bn *a);
int   bn_dec2bn(bn *a, const char *decStr);

void bn_set_hex(bn *x, const char *hexStr) {
    bn_zero(x);
    bn_hex2bn(x, hexStr);
}

/* ------------------ ECC helpers -------------------------------------- */
typedef struct _Point {
    bn x;
    bn y;
} Point;

/* Zero Point P */
void ecp_zero(Point *P) {
    for (int i = 0; i < BN_MAX_WORDS; i++) {
        P->x.w[i] = 0;
        P->y.w[i] = 0;
    }
}

/* Copy Q → P */
void ecp_copy(Point *P, Point *Q) {
    for (int i = 0; i < BN_MAX_WORDS; i++) {
        P->x.w[i] = Q->x.w[i];
        P->y.w[i] = Q->y.w[i];
    }
}

/* Test equality */
int ecp_equal(Point *P, Point *Q) {
    for (int i = 0; i < BN_MAX_WORDS; i++) {
        if (P->x.w[i] != Q->x.w[i] || P->y.w[i] != Q->y.w[i]) return 0;
    }
    return 1;
}

bn inv_mod_p;      /* global Fermat exponent for current field */

/* Point addition / doubling */
void ecp_add(Point *P, Point *Q, Point *R, bn *p, bn *a) {
    bn t1, t2, s;
    Point r;

    /* P at infinity? */
    if (bn_is_zero(&P->x) && bn_is_zero(&P->y)) { ecp_copy(R, Q); return; }

    /* Q at infinity? */
    if (bn_is_zero(&Q->x) && bn_is_zero(&Q->y)) { ecp_copy(R, P); return; }

    /* P ==  Q ?  -> doubling  */
    if (ecp_equal(P, Q)) {
        /* if y == 0 then tangent is vertical, result = infinity */
        if (bn_is_zero(&P->y)) {                         //  y==0 ⇒ O
            ecp_zero(R);
            return;
        }
        bn_mulmod(&t1, &P->x, &P->x, p);
        bn_addmod(&t2, &t1, &t1, p);
        bn_addmod(&t1, &t1, &t2, p);
        bn_addmod(&t1, &t1, a, p);
        bn_addmod(&t2, &P->y, &P->y, p);
    } else {
        /* xP == xQ but yP == -yQ  ->  result = infinity             */
        if (bn_cmp(&P->x, &Q->x) == 0) {                 //  P+(-P)=O
            ecp_zero(R);
            return;
        }
        /* general addition                                         */
        bn_submod(&t1, &P->y, &Q->y, p);
        bn_submod(&t2, &P->x, &Q->x, p);
    }

    /* slope = t1 / t2 mod p */
    bn_modexp(&t2, &t2, &inv_mod_p, p);      /* inverse of t2 */
    bn_mulmod(&s, &t1, &t2, p);

    /* x_r = (s * s) − xP − xQ */
    bn_mulmod(&t1, &s, &s, p);
    bn_addmod(&t2, &P->x, &Q->x, p);
    bn_submod(&r.x, &t1, &t2, p);

    /* y_r = s(xP − x_r) − yP */
    bn_submod(&t1, &P->x, &r.x, p);
    bn_mulmod(&t1, &s, &t1, p);
    bn_submod(&r.y, &t1, &P->y, p);

    ecp_copy(R, &r);
}

/* Point multiplication (double & add, LSB-first) */
void ecp_mul(Point *A, bn *k, Point *R, bn *p, bn *a) {
    Point N;
    ecp_copy(&N, A);
    ecp_zero(R);

    int bits = bn_num_bits(k);
    for (int i = 0; i < bits; i++) {
        if (bn_is_bit_set(k, i)) ecp_add(R, &N, R, p, a);
        ecp_add(&N, &N, &N, p, a);
    }
}

/* ------------------- ECDSA self-test ---------------------------------- */
void test_ecdsa(void) {
    Point pk;                 /* public key */
    Point G;                  /* base point */
    bn tn, n, p, a = {0}, sk;

    /* secp256k1 params */
    bn_hex2bn(&p,
        "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f");
    bn_copy(&inv_mod_p, &p);
    inv_mod_p.b[0] -= 2;                     /* p−2 for Fermat inverse */

    bn_hex2bn(&G.x,
        "79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798");
    bn_hex2bn(&G.y,
        "483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8");

    bn_hex2bn(&n,
        "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141");
    bn_copy(&tn, &n);
    tn.b[0] -= 2;                            /* n−2 for Fermat inverse */

    /* private key */
    bn_hex2bn(&sk,
        "caafcef23ebbef8227a87df538da046bd560162f3dee6f661c649d29797ae29a");

    /* public key pk = k·G */
    ecp_mul(&G, &sk, &pk, &p, &a);
    printf("public key x : %s\n", bn_bn2hex(&pk.x));
    printf("public key y : %s\n", bn_bn2hex(&pk.y));

    /* ---- sign message ---- */
    const char msg[] = "Hello, World!";
    size_t msg_len = strlen(msg);

    sha256_ctx ctx;
    bn z = {0};

    sha256_init(&ctx);
    sha256_update(&ctx, msg, (uint32_t)msg_len);
    sha256_final(z.b, &ctx);                 //  pass pointer to first byte

    /* reverse hash to match vector (demo only) */
    int len = ((bn_num_bits(&z) + 7) & ~7) / 8;
    for (int i = 0; i < len / 2; i++) {
        uint8_t t = z.b[i];
        z.b[i] = z.b[len - i - 1];
        z.b[len - i - 1] = t;
    }

    bn r = {0}, s = {0};
    bn k, kinv, x, y, zs, tmp;
    Point t;
    bn_hex2bn(&k,
        "e26ee863c09bf9880f29f43d08db0c72d6751ece202b74c87b4047b473e862da");

    for (;;) {
        /*   rand_bytes(k.b, 32); */

        ecp_mul(&G, &k, &t, &p, &a);
        bn_mod(&t.x, &n, &r);
        if (bn_is_zero(&r)) continue;

        bn_modexp(&kinv, &k, &tn, &n);            /* kinv = k ^ -1 mod n */
        bn_mulmod(&tmp, &r, &sk, &n);
        bn_addmod(&zs, &z, &tmp, &n);
        bn_mulmod(&s, &kinv, &zs, &n);
        if (!bn_is_zero(&s)) break;
    }

    /* ---- verify ---- */
    bn w, u1, u2, m2;
    Point p1, p2, res;

    bn_modexp(&w, &s, &tn, &n);           /* w = s ^ -1 */
    bn_mulmod(&u1, &z, &w, &n);
    bn_mulmod(&u2, &r, &w, &n);
    ecp_mul(&G,  &u1, &p1, &p, &a);
    ecp_mul(&pk, &u2, &p2, &p, &a);
    ecp_add(&p1, &p2, &res, &p, &a);
    bn_mod(&res.x, &n, &m2);

    printf("Signature verified : %s\n",
           bn_cmp(&m2, &r) == 0 ? "OK" : "FAILED");
}

/* ------------------- ECDH self-test ----------------------------------- */
void test_ecdh(void) {
    printf("Running ECDH test...");

    bn sk, secret, p, a, b;
    Point pk = {0}, G, skey = {0};

    /* curve P-192 prime192v1 */
    bn_hex2bn(&p,
        "fffffffffffffffffffffffffffffffeffffffffffffffff");
    bn_copy(&inv_mod_p, &p);
    inv_mod_p.b[0] -= 2;

    bn_hex2bn(&a,
        "fffffffffffffffffffffffffffffffefffffffffffffffc");
    bn_hex2bn(&b,
        "64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1");

    /* private key for Alice */
    bn_hex2bn(&sk,
        "f17d3fea367b74d340851ca4270dcb24c271f445bed9d527");

    /* Bob's public key (unchecked!) */
    bn_hex2bn(&pk.x,
        "42ea6dd9969dd2a61fea1aac7f8e98edcc896c6e55857cc0");
    bn_hex2bn(&pk.y,
        "dfbe5d7c61fac88b11811bde328e8a0d12bf01a9d204b523");

    /* expected secret */
    bn_hex2bn(&secret,
        "803d8ab2e5b6e6fca715737c3a82f7ce3c783124f6d51cd0");

    /* session key  skey = sk · pk  */
    ecp_mul(&pk, &sk, &skey, &p, &a);

    printf("key : %s\n", bn_bn2hex(&skey.x));
    printf("Session keys %s...\n",
           bn_cmp(&skey.x, &secret) == 0 ? "match" : "don't match!");

    printf("OK.\n");
}

/* ---------------- main ----------------------------------------------- */
int main(int argc, char *argv[]) {
    (void)argc; (void)argv;

    test_ecdsa();
    test_ecdh();
    return 0;
}