y² = x³ + 7 (mod p)

type AffinePoint = {
  x: bigint;  // x-coordinate (mod p)
  y: bigint;  // y-coordinate (mod p)
  infinity: boolean;  // Point at infinity (identity)
};

// Generator point
const G: AffinePoint = {
  x: 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798n,
  y: 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8n,
  infinity: false,
};

// Point at infinity (identity element)
const O: AffinePoint = {
  x: 0n,
  y: 0n,
  infinity: true,
};

function isOnCurve(P: AffinePoint): boolean {
  if (P.infinity) return true;

  const p = SECP256K1_P;
  const y2 = (P.y * P.y) % p;
  const x3_plus_7 = (P.x * P.x * P.x + 7n) % p;

  return y2 === x3_plus_7;
}

function pointAdd(P: AffinePoint, Q: AffinePoint): AffinePoint {
  // Identity cases
  if (P.infinity) return Q;
  if (Q.infinity) return P;

  // Same x-coordinate
  if (P.x === Q.x) {
    if (P.y === Q.y) return pointDouble(P);  // P + P = 2P
    return { x: 0n, y: 0n, infinity: true };  // P + (-P) = O
  }

  const p = SECP256K1_P;

  // Slope: λ = (y₂ - y₁) / (x₂ - x₁)
  const dy = modSub(Q.y, P.y, p);
  const dx = modSub(Q.x, P.x, p);
  const dx_inv = modInv(dx, p);
  const lambda = modMul(dy, dx_inv, p);

  // x₃ = λ² - x₁ - x₂
  const lambda2 = modMul(lambda, lambda, p);
  const x3 = modSub(modSub(lambda2, P.x, p), Q.x, p);

  // y₃ = λ(x₁ - x₃) - y₁
  const x_diff = modSub(P.x, x3, p);
  const y3 = modSub(modMul(lambda, x_diff, p), P.y, p);

  return { x: x3, y: y3, infinity: false };
}

// G + G = 2G
const twoG = pointAdd(G, G);

// Verify result is on curve
assert(isOnCurve(twoG));

function pointDouble(P: AffinePoint): AffinePoint {
  if (P.infinity) return P;

  const p = SECP256K1_P;

  // Slope: λ = 3x² / 2y (derivative of curve equation)
  const x2 = modMul(P.x, P.x, p);
  const three_x2 = modMul(3n, x2, p);
  const two_y = modMul(2n, P.y, p);
  const two_y_inv = modInv(two_y, p);
  const lambda = modMul(three_x2, two_y_inv, p);

  // x₃ = λ² - 2x
  const lambda2 = modMul(lambda, lambda, p);
  const two_x = modMul(2n, P.x, p);
  const x3 = modSub(lambda2, two_x, p);

  // y₃ = λ(x - x₃) - y
  const x_diff = modSub(P.x, x3, p);
  const y3 = modSub(modMul(lambda, x_diff, p), P.y, p);

  return { x: x3, y: y3, infinity: false };
}

// Compute 2G, 4G, 8G, ...
let P = G;
for (let i = 1; i <= 8; i++) {
  console.log(`${1 << i}G:`, P);
  P = pointDouble(P);
}

-P = (x, -y mod p) = (x, p - y)

function pointNegate(P: AffinePoint): AffinePoint {
  if (P.infinity) return P;

  return {
    x: P.x,
    y: SECP256K1_P - P.y,
    infinity: false,
  };
}

// P + (-P) = O
const P = G;
const negP = pointNegate(P);
const sum = pointAdd(P, negP);
assert(sum.infinity === true);  // Identity

R = O  (point at infinity)
Q = P
while k > 0:
  if k is odd:
    R = R + Q
  Q = 2Q  (double)
  k = k >> 1  (shift right)
return R

function scalarMul(k: bigint, P: AffinePoint): AffinePoint {
  if (k === 0n || P.infinity) {
    return { x: 0n, y: 0n, infinity: true };
  }

  let result = { x: 0n, y: 0n, infinity: true };  // O
  let addend = P;
  let scalar = k;

  while (scalar > 0n) {
    if (scalar & 1n) {
      result = pointAdd(result, addend);
    }
    addend = pointDouble(addend);
    scalar >>= 1n;
  }

  return result;
}

// Derive public key from private key
const privateKey = 0x123456789abcdef0123456789abcdef0123456789abcdef0123456789abcdefn;
const publicKey = scalarMul(privateKey, G);

console.log('Public key x:', publicKey.x.toString(16));
console.log('Public key y:', publicKey.y.toString(16));

// Precompute: P, 3P, 5P, 7P, ... (odd multiples)
function precompute(P: AffinePoint, windowSize: number): AffinePoint[] {
  const window = [];
  const twoP = pointDouble(P);

  window.push(P);  // 1P
  for (let i = 1; i < (1 << (windowSize - 1)); i++) {
    window.push(pointAdd(window[i - 1], twoP));
  }

  return window;
}

// Scalar multiplication with wNAF
function scalarMulwNAF(k: bigint, P: AffinePoint, windowSize: number = 4): AffinePoint {
  const precomputed = precompute(P, windowSize);
  const naf = toNAF(k, windowSize);

  let result = { x: 0n, y: 0n, infinity: true };

  for (let i = naf.length - 1; i >= 0; i--) {
    result = pointDouble(result);

    if (naf[i] > 0) {
      result = pointAdd(result, precomputed[(naf[i] - 1) / 2]);
    } else if (naf[i] < 0) {
      result = pointAdd(result, pointNegate(precomputed[(-naf[i] - 1) / 2]));
    }
  }

  return result;
}

// Precompute: G, 2G, 4G, 8G, ..., 2^255 G
const G_multiples: AffinePoint[] = [];
let current = G;
for (let i = 0; i < 256; i++) {
  G_multiples.push(current);
  current = pointDouble(current);
}

// Fast scalar multiplication using precomputed table
function fastMulG(k: bigint): AffinePoint {
  let result = { x: 0n, y: 0n, infinity: true };

  for (let i = 0; i < 256; i++) {
    if ((k >> BigInt(i)) & 1n) {
      result = pointAdd(result, G_multiples[i]);
    }
  }

  return result;
}

Affine (x, y) ↔ Projective (X, Y, Z) where x = X/Z, y = Y/Z

function compressPoint(P: AffinePoint): Uint8Array {
  const compressed = new Uint8Array(33);

  // Prefix: 0x02 (even y) or 0x03 (odd y)
  compressed[0] = (P.y & 1n) === 0n ? 0x02 : 0x03;

  // x-coordinate (32 bytes, big-endian)
  const xBytes = bigIntToBytes(P.x, 32);
  compressed.set(xBytes, 1);

  return compressed;
}

function decompressPoint(compressed: Uint8Array): AffinePoint {
  if (compressed.length !== 33) throw new Error('Invalid compressed point');

  const prefix = compressed[0];
  const x = bytesToBigInt(compressed.slice(1));

  // Solve for y: y² = x³ + 7 mod p
  const p = SECP256K1_P;
  const y2 = modAdd(modPow(x, 3n, p), 7n, p);

  // Compute y = y²^((p+1)/4) mod p (works because p ≡ 3 mod 4)
  const y = modPow(y2, (p + 1n) / 4n, p);

  // Verify y² = y2
  if (modMul(y, y, p) !== y2) throw new Error('Point not on curve');

  // Choose correct y based on prefix
  const yIsOdd = (y & 1n) === 1n;
  const prefixOdd = prefix === 0x03;
  const finalY = yIsOdd === prefixOdd ? y : p - y;

  return { x, y: finalY, infinity: false };
}

if (scalar_bit == 1) {
  result = pointAdd(result, current);  // Timing leak
}

// Conditional add without branching
mask = -(scalar_bit & 1);  // 0 or 0xFFFF...
temp = pointAdd(result, current);
result = conditionalMove(result, temp, mask);