Investigate using complete formulas

[paulmillr]

Renes–Costello–Batina provides complete formulas for addition & doubling of short weistrass curves.

Here's the implementation. Note that it requires homogenous coordinates, not jacobian: x=x/z, y=y/z. Seems like a 10-20% slowdown. Benchmarks (Feb 7, 2022 @ M1 / macos 12):

getPublicKey(utils.randomPrivateKey()) x 5,182 ops/sec @ 192μs/op
sign x 4,035 ops/sec @ 247μs/op
verify x 868 ops/sec @ 1ms/op
recoverPublicKey x 462 ops/sec @ 2ms/op
getSharedSecret aka ecdh x 533 ops/sec @ 1ms/op
getSharedSecret (precomputed) x 5,525 ops/sec @ 180μs/op
Point.fromHex (decompression) x 12,176 ops/sec @ 82μs/op
schnorr.sign x 386 ops/sec @ 2ms/op
schnorr.verify x 469 ops/sec @ 2ms/op

Code:

  // Actually HomogenousPoint
  double(): JacobianPoint {
    const { x, y, z } = this;
    let X3 = _0n, Y3 = _0n, Z3 = _0n; // prettier-ignore
    let b3 = CURVE.b * 3n;
    let t0 = mod(y * y); // step 1
    Z3 = t0 + t0;
    Z3 = Z3 + Z3;
    Z3 = Z3 + Z3;
    let t1 = mod(y * z); // step 5
    let t2 = mod(z * z);
    t2 = mod(b3 * t2);
    X3 = mod(t2 * Z3);
    Y3 = t0 + t2;
    Z3 = mod(t1 * Z3); // step 10
    t1 = t2 + t2;
    t2 = t1 + t2;
    t0 = t0 - t2;
    Y3 = mod(t0 * Y3);
    Y3 = mod(X3 + Y3); // step 15
    t1 = mod(x * y);
    X3 = mod(t0 * t1);
    X3 = mod(X3 + X3);
    return new JacobianPoint(X3, Y3, Z3);
  }
  // Algorithm 7 of 2015/1060
  add(other: JacobianPoint): JacobianPoint {
    const { x: X1, y: Y1, z: Z1 } = this;
    const { x: X2, y: Y2, z: Z2 } = other;
    let X3 = _0n, Y3 = _0n, Z3 = _0n; // prettier-ignore
    let b3 = CURVE.b * 3n;
    //if (this.equals(other)) return this.double();
    let t0 = mod(X1 * X2); // step 1
    let t1 = mod(Y1 * Y2);
    let t2 = mod(Z1 * Z2);
    let t3 = X1 + Y1;
    let t4 = X2 + Y2; // step 5
    t3 = mod(t3 * t4);
    t4 = t0 + t1;
    t3 = t3 - t4;
    t4 = Y1 + Z1;
    X3 = Y2 + Z2; // step 10
    t4 = mod(t4 * X3);
    X3 = t1 + t2;
    t4 = t4 - X3;
    X3 = X1 + Z1;
    Y3 = X2 + Z2; // step 15
    X3 = mod(X3 * Y3);
    Y3 = t0 + t2;
    Y3 = X3 - Y3;
    X3 = t0 + t0;
    t0 = X3 + t0; // step 20
    t2 = mod(b3 * t2);
    Z3 = t1 + t2;
    t1 = t1 - t2;
    Y3 = mod(b3 * Y3);
    X3 = mod(t4 * Y3); // step 25
    t2 = mod(t3 * t1);
    X3 = mod(t2 - X3);
    Y3 = mod(Y3 * t0);
    t1 = mod(t1 * Z3);
    Y3 = mod(t1 + Y3); // step 30
    t0 = mod(t0 * t3);
    Z3 = mod(Z3 * t4);
    Z3 = mod(Z3 + t0);
    return new JacobianPoint(X3, Y3, Z3);
  }
  equals(other: JacobianPoint): boolean {
    if (!(other instanceof JacobianPoint)) throw new TypeError('JacobianPoint expected');
    const { x: X1, y: Y1, z: Z1 } = this;
    const { x: X2, y: Y2, z: Z2 } = other;
    const U1 = mod(X1 * Z2) === mod(X2 * Z1);
    const U2 = mod(Y1 * Z2) === mod(Y2 * Z1);
    return U1 && U2;
  }
  toAffine(invZ: bigint = invert(this.z)): Point {
    const { x, y, z } = this;
    const ax = mod(x * invZ);
    const ay = mod(y * invZ);
    const zz = mod(z * invZ);
    if (zz !== _1n) throw new Error('invZ was invalid');
    return new Point(ax, ay);
  }

[paulmillr]

General case for all short weistrass curves including non -secp256k1

    double(): JacobianPoint {
      const { x, y, z } = this;
      let X3 = _0n, Y3 = _0n, Z3 = _0n; // prettier-ignore
      let b3 = CURVE.b * 3n;
      let t0 = modP(y * y); // step 1
      Z3 = modP(t0 + t0);
      Z3 = modP(Z3 + Z3);
      Z3 = modP(Z3 + Z3);
      let t1 = modP(y * z); // step 5
      let t2 = modP(z * z);
      t2 = modP(b3 * t2);
      X3 = modP(t2 * Z3);
      Y3 = modP(t0 + t2);
      Z3 = modP(t1 * Z3); // step 10
      t1 = modP(t2 + t2);
      t2 = modP(t1 + t2);
      t0 = modP(t0 - t2);
      Y3 = modP(t0 * Y3);
      Y3 = modP(X3 + Y3); // step 15
      t1 = modP(x * y);
      X3 = modP(t0 * t1);
      X3 = modP(X3 + X3);
      return new JacobianPoint(X3, Y3, Z3);
    }
// Cost: 8M + 3S + 2*b3 + 3*a + 15add.
        double2() {
            const { a, b } = CURVE;
            const b3 = modP(_3n * b);
            const { x: X1, y: Y1, z: Z1 } = this;
            // Should work for complete formulas
            //return this.add(this);
            let X3 = _0n, Y3 = _0n, Z3 = _0n; // prettier-ignore
            let t0 = modP(X1 * X1);
            let t1 = modP(Y1 * Y1);
            let t2 = modP(Z1 * Z1);
            let t3 = modP(X1 * Y1);
            Z3 = modP(X1 * Z1);
            t3 = modP(t3 + t3);
            Z3 = modP(Z3 + Z3);
            X3 = modP(a * Z3);
            Y3 = modP(b3 * t2);
            Y3 = modP(X3 + Y3);
            X3 = modP(t1 - Y3);
            Y3 = modP(t1 + Y3);
            Y3 = modP(X3 * Y3);
            X3 = modP(t3 * X3);
            Z3 = modP(b3 * Z3);
            t2 = modP(a * t2);
            t3 = modP(t0 - t2);
            t3 = modP(a * t3);
            t3 = modP(t3 + Z3);
            Z3 = modP(t0 + t0);
            t0 = modP(Z3 + t0);
            t0 = modP(t0 + t2);
            t0 = modP(t0 * t3);
            Y3 = modP(Y3 + t0);
            t2 = modP(Y1 * Z1);
            t2 = modP(t2 + t2);
            t0 = modP(t2 * t3);
            X3 = modP(X3 - t0);
            Z3 = modP(t2 * t1);
            Z3 = modP(Z3 + Z3);
            Z3 = modP(Z3 + Z3);
            return new ProjectivePoint(X3, Y3, Z3);
        }
    // Algorithm 1 of 2015/1060
    add(other: JacobianPoint): JacobianPoint {
      const { x: X1, y: Y1, z: Z1 } = this;
      const { x: X2, y: Y2, z: Z2 } = other;
      let X3 = _0n, Y3 = _0n, Z3 = _0n; // prettier-ignore
      const a = CURVE.a;
      const b3 = CURVE.b * 3n;
      //if (this.equals(other)) return this.double();
      let t0 = modP(X1 * X2); // step 1
      let t1 = modP(Y1 * Y2);
      let t2 = modP(Z1 * Z2);
      let t3 = modP(X1 + Y1);
      let t4 = modP(X2 + Y2); // step 5
      t3 = modP(t3 * t4);
      t4 = modP(t0 + t1);
      t3 = modP(t3 - t4);
      t4 = modP(X1 + Z1);
      let t5 = modP(X2 + Z2); // step 10
      t4 = modP(t4 * t5);
      t5 = modP(t0 + t2);
      t4 = modP(t4 - t5);
      t5 = modP(Y1 + Z1);
      X3 = modP(Y2 + Z2); // step 15
      t5 = modP(t5 * X3);
      X3 = modP(t1 + t2);
      t5 = modP(t5 - X3);
      Z3 = modP(a * t4);
      X3 = modP(b3 * t2); // step 20
      Z3 = modP(X3 + Z3);
      X3 = modP(t1 - Z3);
      Z3 = modP(t1 + Z3);
      Y3 = modP(X3 * Z3);
      t1 = modP(t0 + t0); // step 25
      t1 = modP(t1 + t0);
      t2 = modP(a * t2);
      t4 = modP(b3 * t4);
      t1 = modP(t1 + t2);
      t2 = modP(t0 - t2); // step 30
      t2 = modP(a * t2);
      t4 = modP(t4 + t2);
      t0 = modP(t1 * t4);
      Y3 = modP(Y3 + t0);
      t0 = modP(t5 * t4); // step 35
      X3 = modP(t3 * X3);
      X3 = modP(X3 - t0);
      t0 = modP(t3 * t1);
      Z3 = modP(t5 * Z3);
      Z3 = modP(Z3 + t0); // step 40
      return new JacobianPoint(X3, Y3, Z3);
    }

[paulmillr]

v2

double() {
  const { a, b } = CURVE;
  const b3 = mod(b * 3n);
  const { x: X1, y: Y1, z: Z1 } = this;
  let X3 = 0n, Y3 = 0n, Z3 = 0n; // prettier-ignore
  let t0 = mod(X1 * X1); // step 1
  let t1 = mod(Y1 * Y1);
  let t2 = mod(Z1 * Z1);
  let t3 = mod(X1 * Y1);
  t3 = t3 + t3; // step 5
  Z3 = mod(X1 * Z1);
  Z3 = Z3 + Z3;
  X3 = mod(a * Z3);
  Y3 = mod(b3 * t2);
  Y3 = X3 + Y3; // step 10
  X3 = t1 - Y3; Y3 = t1 + Y3; Y3 = mod(X3 * Y3); X3 = mod(t3 * X3);
  Z3 = mod(b3 * Z3); // step 15
  t2 = mod(a * t2); t3 = t0 - t2; t3 = mod(a * t3); t3 = t3 + Z3;
  Z3 = t0 + t0; // step 20
  t0 = Z3 + t0; t0 = t0 + t2; t0 = mod(t0 * t3); Y3 = mod(Y3 + t0);
  t2 = mod(Y1 * Z1); // step 25
  t2 = t2 + t2; t0 = mod(t2 * t3); X3 = mod(X3 - t0); Z3 = mod(t2 * t1);
  Z3 = Z3 + Z3; // step 30
  Z3 = mod(Z3 + Z3);
  return Point(X3, Y3, Z3);
}