https://paulmillr.com/posts/noble-secp256k1-fast-ecc/

ecc.ts

// PART 1
// secp256k1 curve parameters. Verify using https://www.secg.org/sec2-v2.pdf
const P = 2n ** 256n - 0x1000003d1n; // curve's field prime
const N = 2n ** 256n - 0x14551231950b75fc4402da1732fc9bebfn; // curve (group) order
const CURVE = { P: P, n: N, a: 0n, b: 7n };
const err = (m = ''): never => {
  throw new Error(m);
}; // error helper
const M = (a: bigint, b: bigint = P) => {
  // mod division
  const r = a % b;
  return r >= 0n ? r : b + r; // a % b = r
};
// prettier-ignore
const inv = (num: bigint, md: bigint): bigint => {      // modular inversion
  if (num === 0n || md <= 0n) err('no inverse n=' + num + ' mod=' + md); // no neg exponent for now
  let a = M(num, md), b = md, x = 0n, y = 1n, u = 1n, v = 0n;
  while (a !== 0n) {                                    // uses euclidean gcd algorithm
    const q = b / a, r = b % a;                         // not constant-time
    const m = x - u * q, n = y - v * q;
    b = a, a = r, x = u, y = v, u = m, v = n;
  }
  return b === 1n ? M(x, md) : err('no inverse');       // b is gcd at this point
};
// Point in 2d affine (x, y) coordinates
interface AffinePoint {
  x: bigint;
  y: bigint;
}
const affine = (x: bigint, y: bigint): AffinePoint => ({ x: x, y: y });
const equals = (a: AffinePoint, b: AffinePoint) => a.x === b.x && a.y === b.y;
const Point_ZERO = { x: 0n, y: 0n }; // Point at infinity aka identity point aka zero
// Adds point to itself. https://hyperelliptic.org/EFD/g1p/auto-shortw.html
const double = (a: AffinePoint) => {
  const [X1, Y1] = [a.x, a.y];
  // Calculates slope of the tangent line
  const lam = M(3n * X1 ** 2n * inv(2n * Y1, P)); // λ = (3x₁² + a) / (2y₁)
  const X3 = M(lam * lam - 2n * X1); // x₃ = λ² - 2x₁
  const Y3 = M(lam * (X1 - X3) - Y1); // y₃ = λ * (x₁ - x₃) - y₁
  return affine(X3, Y3);
};
// Adds point to other point. https://hyperelliptic.org/EFD/g1p/auto-shortw.html
const add = (a: AffinePoint, b: AffinePoint): AffinePoint => {
  const [X1, Y1, X2, Y2] = [a.x, a.y, b.x, b.y];
  if (X1 === 0n || Y1 === 0n) return b;
  if (X2 === 0n || Y2 === 0n) return a;
  if (X1 === X2 && Y1 === Y2) return double(a); // special case
  if (X1 === X2 && Y1 === M(-Y2)) return Point_ZERO; // special case
  const lam = M((Y2 - Y1) * inv(X2 - X1, P));
  const X3 = M(lam * lam - X1 - X2);
  const Y3 = M(lam * (X1 - X3) - Y1);
  return affine(X3, Y3);
};

/** Generator / base point */
// G x, y values taken from official secp256k1 document: https://www.secg.org/sec2-v2.pdf
const Gx = 55066263022277343669578718895168534326250603453777594175500187360389116729240n;
const Gy = 32670510020758816978083085130507043184471273380659243275938904335757337482424n;
const Point_BASE = affine(Gx, Gy);
const mul_unsafe = (q: AffinePoint, n: bigint) => {
  let curr = affine(q.x, q.y);
  let p = Point_ZERO;
  while (n > 0n) {
    if (n & 1n) p = add(p, curr);
    curr = double(curr);
    n >>= 1n;
  }
  return p;
};
function getPublicKey_unsafe(privateKey: bigint) {
  return mul_unsafe(Point_BASE, privateKey);
}

// PART 2

interface Signature {
  r: bigint;
  s: bigint;
  recovery?: number;
}
// Convert Uint8Array to bigint, big endian
// prettier-ignore
const bytesToNumBE = (bytes: Uint8Array): bigint => {
  const hex = Array.from(bytes).map((e) => e.toString(16).padStart(2, "0")).join("");
  return hex ? BigInt("0x" + hex) : 0n;
};
// Get random k from CSPRNG: eliminates 0 and keeps modulo bias small
function rand_k_from_1_to_N_small_bias(): bigint {
  // We can't just fetch 32 bytes from CSPRNG: need +16 more to keep modulo bias irrelevant.
  const kbytes = crypto.getRandomValues(new Uint8Array(48));
  // Follow FIPS 186 B.4.1 recomendation to remove 0:
  // instead of doing (x mod N), we do (x mod N-1)+1
  const num = M(bytesToNumBE(kbytes), N - 1n);
  return num + 1n;
}
function sign_slow_unsafe(msgh: Uint8Array, priv: bigint): Signature {
  const m = bytesToNumBE(msgh);
  const d = priv;
  let r = 0n;
  let s = 0n;
  let q;
  do {
    const k = rand_k_from_1_to_N_small_bias();
    const ik = inv(k, N);
    q = mul_unsafe(Point_BASE, k);
    r = M(q.x, N);
    s = M(ik * M(m + d * r, N), N);
  } while (r === 0n || s === 0n);
  if (s > N >> 1n) {
    s = M(-s, N);
  }
  return { r, s };
}
function verify_slow(sig: Signature, msgh: Uint8Array, pub: AffinePoint): boolean {
  const h = M(bytesToNumBE(msgh), N);
  const is = inv(sig.s, N);
  const u1 = M(h * is, N);
  const u2 = M(sig.r * is, N);
  const G_u1 = mul_unsafe(Point_BASE, u1);
  const P_u2 = mul_unsafe(pub, u2);
  const R = add(G_u1, P_u2);
  const v = M(R.x, N);
  return v === sig.r;
}

// PART 3

// Constant Time multiplication
const getPowers = (q: AffinePoint) => {
  let points: AffinePoint[] = [];
  for (let bit = 0, dbl = q; bit <= 256; bit++, dbl = double(dbl)) {
    points.push(dbl);
  }
  return points;
};
const mul_CT_slow = (q: AffinePoint, n: bigint) => {
  const pows = getPowers(q);
  let p = Point_ZERO;
  let f = Point_ZERO; // fake point
  for (let bit = 0; bit <= 256; bit++) {
    const at_bit_i = pows[bit];
    if (n & 1n) p = add(p, at_bit_i);
    else f = add(f, at_bit_i);
    n >>= 1n;
  }
  return p;
};

// PART 5

// Point in 3d projective (x, y, z) coordinates
class Point {
  px: bigint;
  py: bigint;
  pz: bigint;
  constructor(x: bigint, y: bigint, z: bigint) {
    this.px = x;
    this.py = y;
    this.pz = z;
  }
  /** Create 3d xyz point from 2d xy. Edge case: (0, 0) => (0, 1, 0), not (0, 0, 1) */
  static fromAffine(p: AffinePoint): Point {
    return p.x === 0n && p.y === 0n ? new Point(0n, 1n, 0n) : new Point(p.x, p.y, 1n);
  }
  toAffine(iz: bigint = inv(this.pz, P)): AffinePoint {
    const x = M(this.px * iz);
    const y = M(this.py * iz);
    return { x, y };
  }
  // prettier-ignore
  add(other: Point) {
    const { px: X1, py: Y1, pz: Z1 } = this;
    const { px: X2, py: Y2, pz: Z2 } = other;
    const { a, b } = CURVE;
    let X3 = 0n, Y3 = 0n, Z3 = 0n;
    const b3 = M(b * 3n); // step 1
    let t0 = M(X1 * X2), t1 = M(Y1 * Y2), t2 = M(Z1 * Z2), t3 = M(X1 + Y1);
    let t4 = M(X2 + Y2);                                // step 5
    t3 = M(t3 * t4); t4 = M(t0 + t1); t3 = M(t3 - t4); t4 = M(X1 + Z1);
    let t5 = M(X2 + Z2);                                // step 10
    t4 = M(t4 * t5); t5 = M(t0 + t2); t4 = M(t4 - t5); t5 = M(Y1 + Z1);
    X3 = M(Y2 + Z2);                                    // step 15
    t5 = M(t5 * X3); X3 = M(t1 + t2); t5 = M(t5 - X3); Z3 = M(a * t4);
    X3 = M(b3 * t2);                                    // step 20
    Z3 = M(X3 + Z3); X3 = M(t1 - Z3); Z3 = M(t1 + Z3); Y3 = M(X3 * Z3);
    t1 = M(t0 + t0);                                    // step 25
    t1 = M(t1 + t0); t2 = M(a * t2); t4 = M(b3 * t4); t1 = M(t1 + t2);
    t2 = M(t0 - t2);                                    // step 30
    t2 = M(a * t2); t4 = M(t4 + t2); t0 = M(t1 * t4); Y3 = M(Y3 + t0);
    t0 = M(t5 * t4);                                    // step 35
    X3 = M(t3 * X3); X3 = M(X3 - t0); t0 = M(t3 * t1); Z3 = M(t5 * Z3);
    Z3 = M(Z3 + t0);                                    // step 40
    return new Point(X3, Y3, Z3);
  }
  double() {
    return this.add(this);
  }
  negate() {
    return new Point(this.px, M(-this.py), this.pz);
  }
}
const Proj_ZERO = Point.fromAffine(Point_ZERO);
const Proj_BASE = Point.fromAffine(Point_BASE);
const getPowersProj = (qxy: AffinePoint): Point[] => {
  const q = Point.fromAffine(qxy);
  let points: Point[] = [];
  for (let bit = 0, dbl = q; bit <= 256; bit++, dbl = dbl.double()) {
    points.push(dbl);
  }
  return points;
};
const mul_CT = (q: AffinePoint, n: bigint) => {
  const pows = getPowersProj(q);
  let p = Proj_ZERO;
  let f = Proj_ZERO; // fake point
  for (let bit = 0; bit <= 256; bit++) {
    const at_bit_i = pows[bit];
    if (n & 1n) p = p.add(at_bit_i);
    else f = f.add(at_bit_i);
    n >>= 1n;
  }
  return p.toAffine();
};

// PART 6

const W = 8; // Precomputes-related code. W = window size
// prettier-ignore
const precompute = () => {                              // They give 12x faster getPublicKey(),
  const points: Point[] = [];                           // 10x sign(), 2x verify(). To achieve this,
  const windows = 256 / W + 1;                          // app needs to spend 40ms+ to calculate
  let p = Proj_BASE, b = p;                  // a lot of points related to base point G.
  for (let w = 0; w < windows; w++) {                   // Points are stored in array and used
    b = p;                                              // any time Gx multiplication is done.
    points.push(b);                                     // They consume 16-32 MiB of RAM.
    for (let i = 1; i < 2 ** (W - 1); i++) { b = b.add(p); points.push(b); }
    p = b.double();                                     // Precomputes don't speed-up getSharedKey,
  }                                                     // which multiplies user point by scalar,
  return points;                                        // when precomputes are using base point
}
let Gpows: Point[] | undefined = undefined; // precomputes for base point G
// prettier-ignore
const wNAF = (n: bigint): { p: Point; f: Point } => {   // w-ary non-adjacent form (wNAF) method.
                                                        // Compared to other point mult methods,
  const comp = Gpows || (Gpows = precompute());         // stores 2x less points using subtraction
  const neg = (cnd: boolean, p: Point) => { let n = p.negate(); return cnd ? n : p; } // negate
  let p = Proj_ZERO, f = Proj_BASE;// f must be G, or could become I in the end
  const windows = 1 + 256 / W;                          // W=8 17 windows
  const wsize = 2 ** (W - 1);                           // W=8 128 window size
  const mask = BigInt(2 ** W - 1);                      // W=8 will create mask 0b11111111
  const maxNum = 2 ** W;                                // W=8 256
  const shiftBy = BigInt(W);                            // W=8 8
  for (let w = 0; w < windows; w++) {
    const off = w * wsize;
    let wbits = Number(n & mask);                       // extract W bits.
    n >>= shiftBy;                                      // shift number by W bits.
    if (wbits > wsize) { wbits -= maxNum; n += 1n; }    // split if bits > max: +224 => 256-32
    const off1 = off, off2 = off + Math.abs(wbits) - 1; // offsets, evaluate both
    const cnd1 = w % 2 !== 0, cnd2 = wbits < 0;         // conditions, evaluate both
    if (wbits === 0) {
      f = f.add(neg(cnd1, comp[off1]));                 // bits are 0: add garbage to fake point
    } else {                                            //          ^ can't add off2, off2 = I
      p = p.add(neg(cnd2, comp[off2]));                 // bits are 1: add to result point
    }
  }
  return { p, f }                                       // return both real and fake points for JIT
}; // !! you can disable precomputes by commenting-out call of the wNAF() inside Point#mul()
const mul_G_wnaf = (n: bigint) => {
  const { p, f } = wNAF(n);
  f.toAffine(); // result ignored
  return p.toAffine();
};

// PART 7

const CURVE_beta = 0x7ae96a2b657c07106e64479eac3434e99cf0497512f58995c1396c28719501een;
const divNearest = (a: bigint, b: bigint) => (a + b / 2n) / b;
const splitScalar = (k: bigint) => {
  const n = N;
  const a1 = 0x3086d221a7d46bcde86c90e49284eb15n;
  const b1 = -1n * 0xe4437ed6010e88286f547fa90abfe4c3n;
  const a2 = 0x114ca50f7a8e2f3f657c1108d9d44cfd8n;
  const b2 = a1;
  const POW_2_128 = 0x100000000000000000000000000000000n; // (2n**128n).toString(16)

  const c1 = divNearest(b2 * k, n);
  const c2 = divNearest(-b1 * k, n);
  let k1 = M(k - c1 * a1 - c2 * a2, n);
  let k2 = M(-c1 * b1 - c2 * b2, n);
  const k1neg = k1 > POW_2_128;
  const k2neg = k2 > POW_2_128;
  if (k1neg) k1 = n - k1;
  if (k2neg) k2 = n - k2;
  if (k1 > POW_2_128 || k2 > POW_2_128) err('endomorphism failed, k=' + k);
  return { k1neg, k1, k2neg, k2 };
};
const mul_endo = (q: AffinePoint, n: bigint) => {
  let { k1neg, k1, k2neg, k2 } = splitScalar(n);
  let k1p = Proj_ZERO;
  let k2p = Proj_ZERO;
  let d: Point = Point.fromAffine(q);
  while (k1 > 0n || k2 > 0n) {
    if (k1 & 1n) k1p = k1p.add(d);
    if (k2 & 1n) k2p = k2p.add(d);
    d = d.double();
    k1 >>= 1n;
    k2 >>= 1n;
  }
  if (k1neg) k1p = k1p.negate();
  if (k2neg) k2p = k2p.negate();
  k2p = new Point(M(k2p.px * CURVE_beta), k2p.py, k2p.pz);
  return k1p.add(k2p).toAffine();
};

// PART END
const mul = (q: AffinePoint, n: bigint, safe = true) => {
  if (equals(q, Point_BASE)) return mul_G_wnaf(n);
  return safe ? mul_CT(q, n) : mul_endo(q, n);
};

function getPublicKey(privateKey: bigint) {
  return mul(Point_BASE, privateKey);
}

function sign(msgh: Uint8Array, priv: bigint): Signature {
  const m = bytesToNumBE(msgh);
  const d = priv;
  let r = 0n;
  let s = 0n;
  let q;
  do {
    const k = rand_k_from_1_to_N_small_bias();
    const ik = inv(k, N);
    q = mul(Point_BASE, k);
    r = M(q.x, N);
    s = M(ik * M(m + d * r, N), N);
  } while (r === 0n || s === 0n);
  let recovery = (q.x === r ? 0 : 2) | Number(q.y & 1n);
  if (s > N >> 1n) {
    recovery ^= 1;
    s = M(-s, N);
  }
  return { r, s, recovery };
}
function verify(sig: Signature, msgh: Uint8Array, pub: AffinePoint): boolean {
  const h = M(bytesToNumBE(msgh), N);
  const is = inv(sig.s, N);
  const u1 = M(h * is, N);
  const u2 = M(sig.r * is, N);
  const G_u1 = mul(Point_BASE, u1);
  const P_u2 = mul(pub, u2, false);
  const R = add(G_u1, P_u2);
  const v = M(R.x, N);
  return v === sig.r;
}
function getSharedSecret(privA: bigint, pubB: AffinePoint) {
  return mul(pubB, privA);
}

// PART MISC

// PART TEST

const measure = (label: string, fn: any) => {
  console.time(label);
  fn();
  console.timeEnd(label);
};
const test = () => {
  let displayAllLogs = false;

  const P = affine(
    112711660439710606056748659173929673102114977341539408544630613555209775888121n,
    25583027980570883691656905877401976406448868254816295069919888960541586679410n
  );
  const Q = affine(
    21505829891763648114329055987619236494102133314575206970830385799158076338148n,
    98003708678762621233683240503080860129026887322874138805529884920309963580118n
  );
  console.log('test_basic');
  console.log('P + Q', add(P, Q));
  console.log('P + P', double(P));
  // P+Q {
  //   x: 21262057306151627953595685090280431278183829487175876377991189246716355947009n,
  //   y: 41749993296225487051377864631615517161996906063147759678534462689479575333124n
  // }
  // P+P {
  //   x: 115780575977492633039504758427830329241728645270042306223540962614150928364886n,
  //   y: 78735063515800386211891312544505775871260717697865196436804966483607426560663n
  // }

  const G5 = {
    x: 21505829891763648114329055987619236494102133314575206970830385799158076338148n,
    y: 98003708678762621233683240503080860129026887322874138805529884920309963580118n,
  };

  const assert = (a, b) => {
    if (!equals(a, b)) throw new Error(`expected a == b, got a=${a}, b=${b}`);
  };

  function header(label, point) {
    if (displayAllLogs) console.log(label, point);
    assert(point, G5);
  }

  header('mul_DA', mul_unsafe(Point_BASE, 5n));
  measure('small DA', () => mul_unsafe(Point_BASE, 2n));
  measure('large DA', () => mul_unsafe(Point_BASE, 2n ** 255n - 19n));
  // small: 0.078ms
  // large: 14.687ms
  // BASE * 5, Point {
  //   x: 21505829891763648114329055987619236494102133314575206970830385799158076338148n,
  //   y: 98003708678762621233683240503080860129026887322874138805529884920309963580118n
  // }

  header('mul_CT', mul_CT_slow(Point_BASE, 5n));
  measure('small CT', () => mul_CT_slow(Point_BASE, 2n));
  measure('large CT', () => mul_CT_slow(Point_BASE, 2n ** 255n - 19n));
  // small CT: 13.197ms
  // large CT: 13.182ms

  header('mul_CT_proj', mul_CT(Point_BASE, 5n));
  measure('CT proj', () => mul_CT(Point_BASE, 2n ** 255n - 19n));
  // console.log("BASE * 5", mul_CT_proj(Point_BASE, 5n));
  // small CT: 13.197ms
  // large CT: 13.182ms
  header('mul_G_wnaf', mul_G_wnaf(5n));
  measure('mul_G_wnaf', () => mul_G_wnaf(5n));

  header('mul_endo', mul_endo(Point_BASE, 5n));
  measure('mul_endo', () => mul_endo(Point_BASE, 5n));

  header('mul_combined', mul(Point_BASE, 5n));
  measure('mul_combined', () => mul(Point_BASE, 5n));

  const priv = 2n ** 253n - 1230n;
  const msgh = new Uint8Array(32).fill(0xca);
  const pub = getPublicKey(priv);
  const sig = sign(msgh, priv);
  const is_valid = verify_slow(sig, msgh, pub);
  // console.log({
  //   priv,
  //   msgh,
  //   pub,
  //   sig,
  //   is_valid,
  // });
  console.log('is valid sig', is_valid);

  console.log();
  measure('getPublicKey_unsafe 1', () => getPublicKey_unsafe(priv));
  measure('getPublicKey_unsafe 2', () => getPublicKey_unsafe(5n));
  measure('sign_slow_unsafe', () => sign_slow_unsafe(msgh, priv));
  measure('verify_slow', () => verify_slow(sig, msgh, pub));

  console.log();
  measure('getPublicKey', () => getPublicKey(priv));
  measure('sign', () => sign(msgh, priv));
  measure('verify', () => verify(sig, msgh, pub));
};
test();