Release 0.5.2.

Fix lint
tests: Use describe()
2023-01-13 16:23:52 +01:00 · 2023-01-13 16:02:07 +01:00 · 2023-01-13 16:00:13 +01:00 · 2023-01-13 15:58:04 +01:00 · 2023-01-13 01:29:54 +01:00 · 2023-01-13 01:26:00 +01:00
28 changed files with 3993 additions and 3875 deletions
--- a/README.md
+++ b/README.md
@@ -7,8 +7,8 @@ Minimal, auditable JS implementation of elliptic curve cryptography.
 - [hash to curve](https://datatracker.ietf.org/doc/draft-irtf-cfrg-hash-to-curve/)
  for encoding or hashing an arbitrary string to a point on an elliptic curve
 - Auditable, [fast](#speed)
- 🔻 Tree-shaking-friendly: there is no entry point, which ensures small size of your app
 - 🔍 Unique tests ensure correctness. Wycheproof vectors included
+- 🔻 Tree-shaking-friendly: there is no entry point, which ensures small size of your app

 There are two parts of the package:

@@ -84,11 +84,12 @@ To define a custom curve, check out API below.
 ## API

 - [Overview](#overview)
- [abstract/edwards: Twisted Edwards curve](#abstract/edwards-twisted-edwards-curve)
- [abstract/montgomery: Montgomery curve](#abstract/montgomery-montgomery-curve)
- [abstract/weierstrass: Short Weierstrass curve](#abstract/weierstrass-short-weierstrass-curve)
- [abstract/modular](#abstract/modular)
- [abstract/utils](#abstract/utils)
+- [abstract/edwards: Twisted Edwards curve](#abstractedwards-twisted-edwards-curve)
+- [abstract/montgomery: Montgomery curve](#abstractmontgomery-montgomery-curve)
+- [abstract/weierstrass: Short Weierstrass curve](#abstractweierstrass-short-weierstrass-curve)
+- [abstract/hash-to-curve: Hashing strings to curve points](#abstracthash-to-curve-hashing-strings-to-curve-points)
+- [abstract/modular](#abstractmodular)
+- [abstract/utils](#abstractutils)

 ### Overview

@@ -200,8 +201,6 @@ export type CurveFn = {
  ExtendedPoint: ExtendedPointConstructor;
  Signature: SignatureConstructor;
  utils: {
-    mod: (a: bigint, b?: bigint) => bigint;
-    invert: (number: bigint, modulo?: bigint) => bigint;
    randomPrivateKey: () => Uint8Array;
    getExtendedPublicKey: (key: PrivKey) => {
      head: Uint8Array;
@@ -305,6 +304,7 @@ export type CurveFn = {
  getPublicKey: (privateKey: PrivKey, isCompressed?: boolean) => Uint8Array;
  getSharedSecret: (privateA: PrivKey, publicB: PubKey, isCompressed?: boolean) => Uint8Array;
  sign: (msgHash: Hex, privKey: PrivKey, opts?: SignOpts) => SignatureType;
+  signUnhashed: (msg: Uint8Array, privKey: PrivKey, opts?: SignOpts) => SignatureType;
  verify: (
    signature: Hex | SignatureType,
    msgHash: Hex,
@@ -315,8 +315,6 @@ export type CurveFn = {
  ProjectivePoint: ProjectivePointConstructor;
  Signature: SignatureConstructor;
  utils: {
-    mod: (a: bigint) => bigint;
-    invert: (number: bigint) => bigint;
    isValidPrivateKey(privateKey: PrivKey): boolean;
    hashToPrivateKey: (hash: Hex) => Uint8Array;
    randomPrivateKey: () => Uint8Array;
@@ -324,20 +322,70 @@ export type CurveFn = {
 };
 ```

+### abstract/hash-to-curve: Hashing strings to curve points
+
+The module allows to hash arbitrary strings to elliptic curve points.
+
+- `expand_message_xmd` [(spec)](https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-5.4.1) produces a uniformly random byte string using a cryptographic hash function H that outputs b bits..
+
+    ```ts
+    function expand_message_xmd(
+      msg: Uint8Array, DST: Uint8Array, lenInBytes: number, H: CHash
+    ): Uint8Array;
+    function expand_message_xof(
+      msg: Uint8Array, DST: Uint8Array, lenInBytes: number, k: number, H: CHash
+    ): Uint8Array;
+    ```
+
+- `hash_to_field(msg, count, options)` [(spec)](https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-5.3)
+hashes arbitrary-length byte strings to a list of one or more elements of a finite field F.
+    * `msg` a byte string containing the message to hash
+    * `count` the number of elements of F to output
+    * `options` `{DST: string, p: bigint, m: number, k: number, expand: 'xmd' | 'xof', hash: H}`
+    * Returns `[u_0, ..., u_(count - 1)]`, a list of field elements.
+
+    ```ts
+    function hash_to_field(msg: Uint8Array, count: number, options: htfOpts): bigint[][];
+    type htfOpts = {
+      // DST: a domain separation tag
+      // defined in section 2.2.5
+      DST: string;
+      // p: the characteristic of F
+      //    where F is a finite field of characteristic p and order q = p^m
+      p: bigint;
+      // m: the extension degree of F, m >= 1
+      //     where F is a finite field of characteristic p and order q = p^m
+      m: number;
+      // k: the target security level for the suite in bits
+      // defined in section 5.1
+      k: number;
+      // option to use a message that has already been processed by
+      // expand_message_xmd
+      expand?: 'xmd' | 'xof';
+      // Hash functions for: expand_message_xmd is appropriate for use with a
+      // wide range of hash functions, including SHA-2, SHA-3, BLAKE2, and others.
+      // BBS+ uses blake2: https://github.com/hyperledger/aries-framework-go/issues/2247
+      // TODO: verify that hash is shake if expand==='xof' via types
+      hash: CHash;
+    };
+    ```
+
 ### abstract/modular

 Modular arithmetics utilities.

 ```typescript
-import { mod, invert, div, invertBatch, sqrt, Fp } from '@noble/curves/abstract/modular';
+import { Fp, mod, invert, div, invertBatch, sqrt } from '@noble/curves/abstract/modular';
+const fp = Fp(2n ** 255n - 19n); // Finite field over 2^255-19
+fp.mul(591n, 932n);
+fp.pow(481n, 11024858120n);
+
+// Generic non-FP utils are also available
 mod(21n, 10n); // 21 mod 10 == 1n; fixed version of 21 % 10
 invert(17n, 10n); // invert(17) mod 10; modular multiplicative inverse
 div(5n, 17n, 10n); // 5/17 mod 10 == 5 * invert(17) mod 10; division
 invertBatch([1n, 2n, 4n], 21n); // => [1n, 11n, 16n] in one inversion
 sqrt(21n, 73n); // √21 mod 73; square root
-const fp = Fp(2n ** 255n - 19n); // Finite field over 2^255-19
-fp.mul(591n, 932n);
-fp.pow(481n, 11024858120n);
 ```

 ### abstract/utils
--- a/benchmark/index.js
+++ b/benchmark/index.js
@@ -96,6 +96,10 @@ export const CURVES = {
        old_secp.recoverPublicKey(msg, new old_secp.Signature(sig.r, sig.s), sig.recovery),
      secp256k1: ({ sig, msg }) => sig.recoverPublicKey(msg),
    },
+    hashToCurve: {
+      samples: 500,
+      noble: () => secp256k1.Point.hashToCurve('abcd'),
+    },
  },
  ed25519: {
    data: () => {
@@ -124,6 +128,10 @@ export const CURVES = {
      old: ({ sig, msg, pub }) => noble_ed25519.sync.verify(sig, msg, pub),
      noble: ({ sig, msg, pub }) => ed25519.verify(sig, msg, pub),
    },
+    hashToCurve: {
+      samples: 500,
+      noble: () => ed25519.Point.hashToCurve('abcd'),
+    },
  },
  ed448: {
    data: () => {
@@ -145,6 +153,10 @@ export const CURVES = {
      samples: 500,
      noble: ({ sig, msg, pub }) => ed448.verify(sig, msg, pub),
    },
+    hashToCurve: {
+      samples: 500,
+      noble: () => ed448.Point.hashToCurve('abcd'),
+    },
  },
  nist: {
    data: () => {
@@ -168,6 +180,12 @@ export const CURVES = {
      P384: ({ p384: { sig, msg, pub } }) => P384.verify(sig, msg, pub),
      P521: ({ p521: { sig, msg, pub } }) => P521.verify(sig, msg, pub),
    },
+    hashToCurve: {
+      samples: 500,
+      P256: () => P256.Point.hashToCurve('abcd'),
+      P384: () => P384.Point.hashToCurve('abcd'),
+      P521: () => P521.Point.hashToCurve('abcd'),
+    },
  },
  stark: {
    data: () => {
--- a/benchmark/package.json
+++ b/benchmark/package.json
@@ -12,13 +12,15 @@
    "author": "",
    "license": "MIT",
    "devDependencies": {
-        "micro-bmark": "0.2.0"
+        "micro-bmark": "0.2.1"
    },
    "dependencies": {
        "@noble/bls12-381": "^1.4.0",
        "@noble/ed25519": "^1.7.1",
        "@noble/hashes": "^1.1.5",
        "@noble/secp256k1": "^1.7.0",
-        "@starkware-industries/starkware-crypto-utils": "^0.0.2"
+        "@starkware-industries/starkware-crypto-utils": "^0.0.2",
+        "calculate-correlation": "^1.2.3",
+        "elliptic": "^6.5.4"
    }
 }
--- a/package.json
+++ b/package.json
@@ -1,6 +1,6 @@
 {
  "name": "@noble/curves",
-  "version": "0.5.0",
+  "version": "0.5.2",
  "description": "Minimal, auditable JS implementation of elliptic curve cryptography",
  "files": [
    "lib"
@@ -31,7 +31,7 @@
    "@types/node": "18.11.3",
    "fast-check": "3.0.0",
    "micro-bmark": "0.2.0",
-    "micro-should": "0.2.0",
+    "micro-should": "0.3.0",
    "prettier": "2.6.2",
    "rollup": "2.75.5",
    "typescript": "4.7.3"
--- a/src/abstract/bls.ts
+++ b/src/abstract/bls.ts
@@ -1,13 +1,26 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
-// Barreto-Lynn-Scott Curves. A family of pairing friendly curves, with embedding degree = 12 or 24
-// NOTE: only 12 supported for now
-// Constructed from pair of weierstrass curves, based pairing logic
+/**
+ * BLS (Barreto-Lynn-Scott) family of pairing-friendly curves.
+ * Implements BLS (Boneh-Lynn-Shacham) signatures.
+ * Consists of two curves: G1 and G2:
+ * - G1 is a subgroup of (x, y) E(Fq) over y² = x³ + 4.
+ * - G2 is a subgroup of ((x₁, x₂+i), (y₁, y₂+i)) E(Fq²) over y² = x³ + 4(1 + i) where i is √-1
+ * - Gt, created by bilinear (ate) pairing e(G1, G2), consists of p-th roots of unity in
+ *   Fq^k where k is embedding degree. Only degree 12 is currently supported, 24 is not.
+ * Pairing is used to aggregate and verify signatures.
+ * We are using Fp for private keys (shorter) and Fp₂ for signatures (longer).
+ * Some projects may prefer to swap this relation, it is not supported for now.
+ */
 import * as mod from './modular.js';
-import { ensureBytes, numberToBytesBE, bytesToNumberBE, bitLen, bitGet } from './utils.js';
-import * as utils from './utils.js';
-// Types
-import { hexToBytes, bytesToHex, Hex, PrivKey } from './utils.js';
-import { htfOpts, stringToBytes, hash_to_field, expand_message_xmd } from './hash-to-curve.js';
+import * as ut from './utils.js';
+// Types require separate import
+import { Hex, PrivKey } from './utils.js';
+import {
+  htfOpts,
+  stringToBytes,
+  hash_to_field as hashToField,
+  expand_message_xmd as expandMessageXMD,
+} from './hash-to-curve.js';
 import { CurvePointsType, PointType, CurvePointsRes, weierstrassPoints } from './weierstrass.js';

 type Fp = bigint; // Can be different field?
@@ -39,7 +52,7 @@ export type CurveType<Fp, Fp2, Fp6, Fp12> = {
    finalExponentiate(num: Fp12): Fp12;
  };
  htfDefaults: htfOpts;
-  hash: utils.CHash; // Because we need outputLen for DRBG
+  hash: ut.CHash; // Because we need outputLen for DRBG
  randomBytes: (bytesLength?: number) => Uint8Array;
 };

@@ -80,12 +93,9 @@ export type CurveFn<Fp, Fp2, Fp6, Fp12> = {
    publicKeys: (Hex | PointType<Fp>)[]
  ) => boolean;
  utils: {
-    bytesToHex: typeof utils.bytesToHex;
-    hexToBytes: typeof utils.hexToBytes;
    stringToBytes: typeof stringToBytes;
-    hashToField: typeof hash_to_field;
-    expandMessageXMD: typeof expand_message_xmd;
-    mod: typeof mod.mod;
+    hashToField: typeof hashToField;
+    expandMessageXMD: typeof expandMessageXMD;
    getDSTLabel: () => string;
    setDSTLabel(newLabel: string): void;
  };
@@ -95,12 +105,9 @@ export function bls<Fp2, Fp6, Fp12>(
  CURVE: CurveType<Fp, Fp2, Fp6, Fp12>
 ): CurveFn<Fp, Fp2, Fp6, Fp12> {
  // Fields looks pretty specific for curve, so for now we need to pass them with options
-  const Fp = CURVE.Fp;
-  const Fr = CURVE.Fr;
-  const Fp2 = CURVE.Fp2;
-  const Fp6 = CURVE.Fp6;
-  const Fp12 = CURVE.Fp12;
-  const BLS_X_LEN = bitLen(CURVE.x);
+  const { Fp, Fr, Fp2, Fp6, Fp12 } = CURVE;
+  const BLS_X_LEN = ut.bitLen(CURVE.x);
+  const groupLen = 32; // TODO: calculate; hardcoded for now

  // Pre-compute coefficients for sparse multiplication
  // Point addition and point double calculations is reused for coefficients
@@ -125,7 +132,7 @@ export function bls<Fp2, Fp6, Fp12>(
      Rx = Fp2.div(Fp2.mul(Fp2.mul(Fp2.sub(t0, t3), Rx), Ry), 2n); // ((T0 - T3) * Rx * Ry) / 2
      Ry = Fp2.sub(Fp2.square(Fp2.div(Fp2.add(t0, t3), 2n)), Fp2.mul(Fp2.square(t2), 3n)); // ((T0 + T3) / 2)² - 3 * T2²
      Rz = Fp2.mul(t0, t4); // T0 * T4
-      if (bitGet(CURVE.x, i)) {
+      if (ut.bitGet(CURVE.x, i)) {
        // Addition
        let t0 = Fp2.sub(Ry, Fp2.mul(Qy, Rz)); // Ry - Qy * Rz
        let t1 = Fp2.sub(Rx, Fp2.mul(Qx, Rz)); // Rx - Qx * Rz
@@ -147,13 +154,14 @@ export function bls<Fp2, Fp6, Fp12>(
  }

  function millerLoop(ell: [Fp2, Fp2, Fp2][], g1: [Fp, Fp]): Fp12 {
+    const { x } = CURVE;
    const Px = g1[0];
    const Py = g1[1];
    let f12 = Fp12.ONE;
    for (let j = 0, i = BLS_X_LEN - 2; i >= 0; i--, j++) {
      const E = ell[j];
      f12 = Fp12.multiplyBy014(f12, E[0], Fp2.mul(E[1], Px), Fp2.mul(E[2], Py));
-      if (bitGet(CURVE.x, i)) {
+      if (ut.bitGet(x, i)) {
        j += 1;
        const F = ell[j];
        f12 = Fp12.multiplyBy014(f12, F[0], Fp2.mul(F[1], Px), Fp2.mul(F[2], Py));
@@ -163,81 +171,31 @@ export function bls<Fp2, Fp6, Fp12>(
    return Fp12.conjugate(f12);
  }

-  // bls12-381 is a construction of two curves:
-  // 1. Fp: (x, y)
-  // 2. Fp₂: ((x₁, x₂+i), (y₁, y₂+i)) - (complex numbers)
-  //
-  // Bilinear Pairing (ate pairing) is used to combine both elements into a paired one:
-  //   Fp₁₂ = e(Fp, Fp2)
-  //   where Fp₁₂ = 12-degree polynomial
-  // Pairing is used to verify signatures.
-  //
-  // We are using Fp for private keys (shorter) and Fp2 for signatures (longer).
-  // Some projects may prefer to swap this relation, it is not supported for now.
-
-  const htfDefaults = { ...CURVE.htfDefaults };
-
-  function isWithinCurveOrder(num: bigint): boolean {
-    return 0 < num && num < CURVE.r;
-  }
-
  const utils = {
-    hexToBytes: hexToBytes,
-    bytesToHex: bytesToHex,
-    mod: mod.mod,
-    stringToBytes,
+    hexToBytes: ut.hexToBytes,
+    bytesToHex: ut.bytesToHex,
+    stringToBytes: stringToBytes,
    // TODO: do we need to export it here?
-    hashToField: (msg: Uint8Array, count: number, options: Partial<typeof htfDefaults> = {}) =>
-      hash_to_field(msg, count, { ...CURVE.htfDefaults, ...options }),
+    hashToField: (
+      msg: Uint8Array,
+      count: number,
+      options: Partial<typeof CURVE.htfDefaults> = {}
+    ) => hashToField(msg, count, { ...CURVE.htfDefaults, ...options }),
    expandMessageXMD: (msg: Uint8Array, DST: Uint8Array, lenInBytes: number, H = CURVE.hash) =>
-      expand_message_xmd(msg, DST, lenInBytes, H),
-
-    /**
-     * Can take 40 or more bytes of uniform input e.g. from CSPRNG or KDF
-     * and convert them into private key, with the modulo bias being negligible.
-     * As per FIPS 186 B.1.1.
-     * https://research.kudelskisecurity.com/2020/07/28/the-definitive-guide-to-modulo-bias-and-how-to-avoid-it/
-     * @param hash hash output from sha512, or a similar function
-     * @returns valid private key
-     */
-    hashToPrivateKey: (hash: Hex): Uint8Array => {
-      hash = ensureBytes(hash);
-      if (hash.length < 40 || hash.length > 1024)
-        throw new Error('Expected 40-1024 bytes of private key as per FIPS 186');
-      //     hashToPrivateScalar(hash, CURVE.r)
-      // NOTE: doesn't add +/-1
-      const num = mod.mod(bytesToNumberBE(hash), CURVE.r);
-      // This should never happen
-      if (num === 0n || num === 1n) throw new Error('Invalid private key');
-      return numberToBytesBE(num, 32);
-    },
-
-    randomBytes: (bytesLength: number = 32): Uint8Array => CURVE.randomBytes(bytesLength),
-    // NIST SP 800-56A rev 3, section 5.6.1.2.2
-    // https://research.kudelskisecurity.com/2020/07/28/the-definitive-guide-to-modulo-bias-and-how-to-avoid-it/
-    randomPrivateKey: (): Uint8Array => utils.hashToPrivateKey(utils.randomBytes(40)),
-    getDSTLabel: () => htfDefaults.DST,
+      expandMessageXMD(msg, DST, lenInBytes, H),
+    hashToPrivateKey: (hash: Hex): Uint8Array => Fr.toBytes(ut.hashToPrivateScalar(hash, CURVE.r)),
+    randomBytes: (bytesLength: number = groupLen): Uint8Array => CURVE.randomBytes(bytesLength),
+    randomPrivateKey: (): Uint8Array => utils.hashToPrivateKey(utils.randomBytes(groupLen + 8)),
+    getDSTLabel: () => CURVE.htfDefaults.DST,
    setDSTLabel(newLabel: string) {
      // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-3.1
      if (typeof newLabel !== 'string' || newLabel.length > 2048 || newLabel.length === 0) {
        throw new TypeError('Invalid DST');
      }
-      htfDefaults.DST = newLabel;
+      CURVE.htfDefaults.DST = newLabel;
    },
  };

-  function normalizePrivKey(key: PrivKey): bigint {
-    let int: bigint;
-    if (key instanceof Uint8Array && key.length === 32) int = bytesToNumberBE(key);
-    else if (typeof key === 'string' && key.length === 64) int = BigInt(`0x${key}`);
-    else if (typeof key === 'number' && key > 0 && Number.isSafeInteger(key)) int = BigInt(key);
-    else if (typeof key === 'bigint' && key > 0n) int = key;
-    else throw new TypeError('Expected valid private key');
-    int = mod.mod(int, CURVE.r);
-    if (!isWithinCurveOrder(int)) throw new Error('Private key must be 0 < key < CURVE.r');
-    return int;
-  }
-
  // Point on G1 curve: (x, y)
  const G1 = weierstrassPoints({
    n: Fr.ORDER,
@@ -309,7 +267,7 @@ export function bls<Fp2, Fp6, Fp12>(
  function sign(message: G2Hex, privateKey: PrivKey): Uint8Array | G2 {
    const msgPoint = normP2Hash(message);
    msgPoint.assertValidity();
-    const sigPoint = msgPoint.multiply(normalizePrivKey(privateKey));
+    const sigPoint = msgPoint.multiply(G1.normalizePrivateKey(privateKey));
    if (message instanceof G2.Point) return sigPoint;
    return Signature.encode(sigPoint);
  }
--- a/src/abstract/edwards.ts
+++ b/src/abstract/edwards.ts
@@ -2,28 +2,18 @@
 // Twisted Edwards curve. The formula is: ax² + y² = 1 + dx²y²

 // Differences from @noble/ed25519 1.7:
-// 1. Different field element lengths in ed448:
+// 1. Variable field element lengths between EDDSA/ECDH:
 //   EDDSA (RFC8032) is 456 bits / 57 bytes, ECDH (RFC7748) is 448 bits / 56 bytes
 // 2. Different addition formula (doubling is same)
 // 3. uvRatio differs between curves (half-expected, not only pow fn changes)
-// 4. Point decompression code is different too (unexpected), now using generalized formula
+// 4. Point decompression code is different (unexpected), now using generalized formula
 // 5. Domain function was no-op for ed25519, but adds some data even with empty context for ed448

 import * as mod from './modular.js';
-import {
-  bytesToHex,
-  concatBytes,
-  ensureBytes,
-  numberToBytesLE,
-  bytesToNumberLE,
-  hashToPrivateScalar,
-  BasicCurve,
-  validateOpts as utilOpts,
-  Hex,
-  PrivKey,
-} from './utils.js'; // TODO: import * as u from './utils.js'?
+import * as ut from './utils.js';
+import { ensureBytes, Hex, PrivKey } from './utils.js';
 import { Group, GroupConstructor, wNAF } from './group.js';
-import { hash_to_field, htfOpts, validateHTFOpts } from './hash-to-curve.js';
+import { hash_to_field as hashToField, htfOpts, validateHTFOpts } from './hash-to-curve.js';

 // Be friendly to bad ECMAScript parsers by not using bigint literals like 123n
 const _0n = BigInt(0);
@@ -31,49 +21,41 @@ const _1n = BigInt(1);
 const _2n = BigInt(2);
 const _8n = BigInt(8);

-export type CHash = {
-  (message: Uint8Array | string): Uint8Array;
-  blockLen: number;
-  outputLen: number;
-  create(): any;
-};
-
-export type CurveType = BasicCurve<bigint> & {
+// Edwards curves must declare params a & d.
+export type CurveType = ut.BasicCurve<bigint> & {
  // Params: a, d
  a: bigint;
  d: bigint;
  // Hashes
-  hash: CHash; // Because we need outputLen for DRBG
+  // The interface, because we need outputLen for DRBG
+  hash: ut.CHash;
+  // CSPRNG
  randomBytes: (bytesLength?: number) => Uint8Array;
+  // Probably clears bits in a byte array to produce a valid field element
  adjustScalarBytes?: (bytes: Uint8Array) => Uint8Array;
+  // Used during hashing
  domain?: (data: Uint8Array, ctx: Uint8Array, phflag: boolean) => Uint8Array;
+  // Ratio √(u/v)
  uvRatio?: (u: bigint, v: bigint) => { isValid: boolean; value: bigint };
-  preHash?: CHash;
-  clearCofactor?: (c: ExtendedPointConstructor, point: ExtendedPointType) => ExtendedPointType;
-  // Hash to field opts
+  // RFC 8032 pre-hashing of messages to sign() / verify()
+  preHash?: ut.CHash;
+  // Hash to field options
  htfDefaults?: htfOpts;
  mapToCurve?: (scalar: bigint[]) => { x: bigint; y: bigint };
 };

-// Should be separate from overrides, since overrides can use information about curve (for example nBits)
 function validateOpts(curve: CurveType) {
-  const opts = utilOpts(curve);
-  if (typeof opts.hash !== 'function' || !Number.isSafeInteger(opts.hash.outputLen))
+  const opts = ut.validateOpts(curve);
+  if (typeof opts.hash !== 'function' || !ut.isPositiveInt(opts.hash.outputLen))
    throw new Error('Invalid hash function');
  for (const i of ['a', 'd'] as const) {
-    if (typeof opts[i] !== 'bigint')
-      throw new Error(`Invalid curve param ${i}=${opts[i]} (${typeof opts[i]})`);
+    const val = opts[i];
+    if (typeof val !== 'bigint') throw new Error(`Invalid curve param ${i}=${val} (${typeof val})`);
  }
  for (const fn of ['randomBytes'] as const) {
    if (typeof opts[fn] !== 'function') throw new Error(`Invalid ${fn} function`);
  }
-  for (const fn of [
-    'adjustScalarBytes',
-    'domain',
-    'uvRatio',
-    'mapToCurve',
-    'clearCofactor',
-  ] as const) {
+  for (const fn of ['adjustScalarBytes', 'domain', 'uvRatio', 'mapToCurve'] as const) {
    if (opts[fn] === undefined) continue; // Optional
    if (typeof opts[fn] !== 'function') throw new Error(`Invalid ${fn} function`);
  }
@@ -96,7 +78,7 @@ export type SignatureConstructor = {
  fromHex(hex: Hex): SignatureType;
 };

-// Instance
+// Instance of Extended Point with coordinates in X, Y, Z, T
 export interface ExtendedPointType extends Group<ExtendedPointType> {
  readonly x: bigint;
  readonly y: bigint;
@@ -109,7 +91,7 @@ export interface ExtendedPointType extends Group<ExtendedPointType> {
  toAffine(invZ?: bigint): PointType;
  clearCofactor(): ExtendedPointType;
 }
-// Static methods
+// Static methods of Extended Point with coordinates in X, Y, Z, T
 export interface ExtendedPointConstructor extends GroupConstructor<ExtendedPointType> {
  new (x: bigint, y: bigint, z: bigint, t: bigint): ExtendedPointType;
  fromAffine(p: PointType): ExtendedPointType;
@@ -117,7 +99,7 @@ export interface ExtendedPointConstructor extends GroupConstructor<ExtendedPoint
  normalizeZ(points: ExtendedPointType[]): ExtendedPointType[];
 }

-// Instance
+// Instance of Affine Point with coordinates in X, Y
 export interface PointType extends Group<PointType> {
  readonly x: bigint;
  readonly y: bigint;
@@ -127,7 +109,7 @@ export interface PointType extends Group<PointType> {
  isTorsionFree(): boolean;
  clearCofactor(): PointType;
 }
-// Static methods
+// Static methods of Affine Point with coordinates in X, Y
 export interface PointConstructor extends GroupConstructor<PointType> {
  new (x: bigint, y: bigint): PointType;
  fromHex(hex: Hex): PointType;
@@ -148,8 +130,6 @@ export type CurveFn = {
  ExtendedPoint: ExtendedPointConstructor;
  Signature: SignatureConstructor;
  utils: {
-    mod: (a: bigint) => bigint;
-    invert: (number: bigint) => bigint;
    randomPrivateKey: () => Uint8Array;
    getExtendedPublicKey: (key: PrivKey) => {
      head: Uint8Array;
@@ -164,36 +144,31 @@ export type CurveFn = {
 // NOTE: it is not generic twisted curve for now, but ed25519/ed448 generic implementation
 export function twistedEdwards(curveDef: CurveType): CurveFn {
  const CURVE = validateOpts(curveDef) as ReturnType<typeof validateOpts>;
-  const Fp = CURVE.Fp as mod.Field<bigint>;
+  const Fp = CURVE.Fp;
  const CURVE_ORDER = CURVE.n;
-  const fieldLen = Fp.BYTES; // 32 (length of one field element)
-  if (fieldLen > 2048) throw new Error('Field lengths over 2048 are not supported');
-  const groupLen = CURVE.nByteLength;
-  // (2n ** 256n).toString(16);
-  const maxGroupElement = _2n ** BigInt(groupLen * 8); // previous POW_2_256
+  const maxGroupElement = _2n ** BigInt(CURVE.nByteLength * 8);

  // Function overrides
  const { randomBytes } = CURVE;
  const modP = Fp.create;

  // sqrt(u/v)
-  function _uvRatio(u: bigint, v: bigint) {
+  const uvRatio =
+    CURVE.uvRatio ||
+    ((u: bigint, v: bigint) => {
      try {
-      const value = Fp.sqrt(u * Fp.invert(v));
-      return { isValid: true, value };
+        return { isValid: true, value: Fp.sqrt(u * Fp.invert(v)) };
      } catch (e) {
        return { isValid: false, value: _0n };
      }
-  }
-  const uvRatio = CURVE.uvRatio || _uvRatio;
-
-  const _adjustScalarBytes = (bytes: Uint8Array) => bytes; // NOOP
-  const adjustScalarBytes = CURVE.adjustScalarBytes || _adjustScalarBytes;
-  function _domain(data: Uint8Array, ctx: Uint8Array, phflag: boolean) {
+    });
+  const adjustScalarBytes = CURVE.adjustScalarBytes || ((bytes: Uint8Array) => bytes); // NOOP
+  const domain =
+    CURVE.domain ||
+    ((data: Uint8Array, ctx: Uint8Array, phflag: boolean) => {
      if (ctx.length || phflag) throw new Error('Contexts/pre-hash are not supported');
      return data;
-  }
-  const domain = CURVE.domain || _domain; // NOOP
+    }); // NOOP

  /**
   * Extended Point works in extended coordinates: (x, y, z, t) ∋ (x=x/z, y=y/z, t=xy).
@@ -336,25 +311,27 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    // Non-constant-time multiplication. Uses double-and-add algorithm.
    // It's faster, but should only be used when you don't care about
    // an exposed private key e.g. sig verification.
-    // Allows scalar bigger than curve order, but less than 2^256
    multiplyUnsafe(scalar: number | bigint): ExtendedPoint {
      let n = normalizeScalar(scalar, CURVE_ORDER, false);
-      const G = ExtendedPoint.BASE;
      const P0 = ExtendedPoint.ZERO;
      if (n === _0n) return P0;
      if (this.equals(P0) || n === _1n) return this;
-      if (this.equals(G)) return this.wNAF(n);
+      if (this.equals(ExtendedPoint.BASE)) return this.wNAF(n);
      return wnaf.unsafeLadder(this, n);
    }

+    // Checks if point is of small order.
+    // If you add something to small order point, you will have "dirty"
+    // point with torsion component.
    // Multiplies point by cofactor and checks if the result is 0.
    isSmallOrder(): boolean {
      return this.multiplyUnsafe(CURVE.h).equals(ExtendedPoint.ZERO);
    }

-    // Multiplies point by a very big scalar n and checks if the result is 0.
+    // Multiplies point by curve order (very big scalar CURVE.n) and checks if the result is 0.
+    // Returns `false` is the point is dirty.
    isTorsionFree(): boolean {
-      return this.multiplyUnsafe(CURVE_ORDER).equals(ExtendedPoint.ZERO);
+      return wnaf.unsafeLadder(this, CURVE_ORDER).equals(ExtendedPoint.ZERO);
    }

    // Converts Extended point to default (x, y) coordinates.
@@ -371,14 +348,12 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
      return new Point(ax, ay);
    }
    clearCofactor(): ExtendedPoint {
-      if (CURVE.h === _1n) return this; // Fast-path
-      // clear_cofactor(P) := h_eff * P
-      // hEff = h for ed25519/ed448. Maybe worth moving to params?
-      if (CURVE.clearCofactor) return CURVE.clearCofactor(ExtendedPoint, this) as ExtendedPoint;
-      return this.multiplyUnsafe(CURVE.h);
+      const { h: cofactor } = CURVE;
+      if (cofactor === _1n) return this;
+      return this.multiplyUnsafe(cofactor);
    }
  }
-  const wnaf = wNAF(ExtendedPoint, groupLen * 8);
+  const wnaf = wNAF(ExtendedPoint, CURVE.nByteLength * 8);

  function assertExtPoint(other: unknown) {
    if (!(other instanceof ExtendedPoint)) throw new TypeError('ExtendedPoint expected');
@@ -413,19 +388,20 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    // Uses algo from RFC8032 5.1.3.
    static fromHex(hex: Hex, strict = true) {
      const { d, a } = CURVE;
-      hex = ensureBytes(hex, fieldLen);
+      const len = Fp.BYTES;
+      hex = ensureBytes(hex, len);
      // 1.  First, interpret the string as an integer in little-endian
      // representation. Bit 255 of this number is the least significant
      // bit of the x-coordinate and denote this value x_0.  The
      // y-coordinate is recovered simply by clearing this bit.  If the
      // resulting value is >= p, decoding fails.
      const normed = hex.slice();
-      const lastByte = hex[fieldLen - 1];
-      normed[fieldLen - 1] = lastByte & ~0x80;
-      const y = bytesToNumberLE(normed);
+      const lastByte = hex[len - 1];
+      normed[len - 1] = lastByte & ~0x80;
+      const y = ut.bytesToNumberLE(normed);

      if (strict && y >= Fp.ORDER) throw new Error('Expected 0 < hex < P');
-      if (!strict && y >= maxGroupElement) throw new Error('Expected 0 < hex < 2**256');
+      if (!strict && y >= maxGroupElement) throw new Error('Expected 0 < hex < CURVE.n');

      // 2.  To recover the x-coordinate, the curve equation implies
      // Ed25519: x² = (y² - 1) / (d y² + 1) (mod p).
@@ -459,16 +435,18 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    // When compressing point, it's enough to only store its y coordinate
    // and use the last byte to encode sign of x.
    toRawBytes(): Uint8Array {
-      const bytes = numberToBytesLE(this.y, fieldLen);
-      bytes[fieldLen - 1] |= this.x & _1n ? 0x80 : 0;
+      const bytes = ut.numberToBytesLE(this.y, Fp.BYTES);
+      bytes[Fp.BYTES - 1] |= this.x & _1n ? 0x80 : 0;
      return bytes;
    }

    // Same as toRawBytes, but returns string.
    toHex(): string {
-      return bytesToHex(this.toRawBytes());
+      return ut.bytesToHex(this.toRawBytes());
    }

+    // Determines if point is in prime-order subgroup.
+    // Returns `false` is the point is dirty.
    isTorsionFree(): boolean {
      return ExtendedPoint.fromAffine(this).isTorsionFree();
    }
@@ -509,20 +487,20 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    // Encodes byte string to elliptic curve
    // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-3
    static hashToCurve(msg: Hex, options?: Partial<htfOpts>) {
-      if (!CURVE.mapToCurve) throw new Error('No mapToCurve defined for curve');
-      msg = ensureBytes(msg);
-      const u = hash_to_field(msg, 2, { ...CURVE.htfDefaults, ...options } as htfOpts);
-      const { x: x0, y: y0 } = CURVE.mapToCurve(u[0]);
-      const { x: x1, y: y1 } = CURVE.mapToCurve(u[1]);
+      const { mapToCurve, htfDefaults } = CURVE;
+      if (!mapToCurve) throw new Error('No mapToCurve defined for curve');
+      const u = hashToField(ensureBytes(msg), 2, { ...htfDefaults, ...options } as htfOpts);
+      const { x: x0, y: y0 } = mapToCurve(u[0]);
+      const { x: x1, y: y1 } = mapToCurve(u[1]);
      const p = new Point(x0, y0).add(new Point(x1, y1)).clearCofactor();
      return p;
    }
    // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-16#section-3
    static encodeToCurve(msg: Hex, options?: Partial<htfOpts>) {
-      if (!CURVE.mapToCurve) throw new Error('No mapToCurve defined for curve');
-      msg = ensureBytes(msg);
-      const u = hash_to_field(msg, 1, { ...CURVE.htfDefaults, ...options } as htfOpts);
-      const { x, y } = CURVE.mapToCurve(u[0]);
+      const { mapToCurve, htfDefaults } = CURVE;
+      if (!mapToCurve) throw new Error('No mapToCurve defined for curve');
+      const u = hashToField(ensureBytes(msg), 1, { ...htfDefaults, ...options } as htfOpts);
+      const { x, y } = mapToCurve(u[0]);
      return new Point(x, y).clearCofactor();
    }
  }
@@ -536,9 +514,10 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    }

    static fromHex(hex: Hex) {
-      const bytes = ensureBytes(hex, 2 * fieldLen);
-      const r = Point.fromHex(bytes.slice(0, fieldLen), false);
-      const s = bytesToNumberLE(bytes.slice(fieldLen, 2 * fieldLen));
+      const len = Fp.BYTES;
+      const bytes = ensureBytes(hex, 2 * len);
+      const r = Point.fromHex(bytes.slice(0, len), false);
+      const s = ut.bytesToNumberLE(bytes.slice(len, 2 * len));
      return new Signature(r, s);
    }

@@ -551,17 +530,17 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    }

    toRawBytes() {
-      return concatBytes(this.r.toRawBytes(), numberToBytesLE(this.s, fieldLen));
+      return ut.concatBytes(this.r.toRawBytes(), ut.numberToBytesLE(this.s, Fp.BYTES));
    }

    toHex() {
-      return bytesToHex(this.toRawBytes());
+      return ut.bytesToHex(this.toRawBytes());
    }
  }

  // Little-endian SHA512 with modulo n
-  function modlLE(hash: Uint8Array): bigint {
-    return mod.mod(bytesToNumberLE(hash), CURVE_ORDER);
+  function modnLE(hash: Uint8Array): bigint {
+    return mod.mod(ut.bytesToNumberLE(hash), CURVE_ORDER);
  }

  /**
@@ -572,7 +551,7 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
   */
  function normalizeScalar(num: number | bigint, max: bigint, strict = true): bigint {
    if (!max) throw new TypeError('Specify max value');
-    if (typeof num === 'number' && Number.isSafeInteger(num)) num = BigInt(num);
+    if (ut.isPositiveInt(num)) num = BigInt(num);
    if (typeof num === 'bigint' && num < max) {
      if (strict) {
        if (_0n < num) return num;
@@ -580,37 +559,32 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
        if (_0n <= num) return num;
      }
    }
-    throw new TypeError('Expected valid scalar: 0 < scalar < max');
-  }
-
-  function checkPrivateKey(key: PrivKey) {
-    // Normalize bigint / number / string to Uint8Array
-    key =
-      typeof key === 'bigint' || typeof key === 'number'
-        ? numberToBytesLE(normalizeScalar(key, maxGroupElement), groupLen)
-        : ensureBytes(key);
-    if (key.length !== groupLen) throw new Error(`Expected ${groupLen} bytes, got ${key.length}`);
-    return key;
-  }
-
-  // Takes 64 bytes
-  function getKeyFromHash(hashed: Uint8Array) {
-    // First 32 bytes of 64b uniformingly random input are taken,
-    // clears 3 bits of it to produce a random field element.
-    const head = adjustScalarBytes(hashed.slice(0, groupLen));
-    // Second 32 bytes is called key prefix (5.1.6)
-    const prefix = hashed.slice(groupLen, 2 * groupLen);
-    // The actual private scalar
-    const scalar = modlLE(head);
-    // Point on Edwards curve aka public key
-    const point = Point.BASE.multiply(scalar);
-    const pointBytes = point.toRawBytes();
-    return { head, prefix, scalar, point, pointBytes };
+    throw new TypeError(`Expected valid scalar: 0 < scalar < ${max}`);
  }

  /** Convenience method that creates public key and other stuff. RFC8032 5.1.5 */
  function getExtendedPublicKey(key: PrivKey) {
-    return getKeyFromHash(CURVE.hash(checkPrivateKey(key)));
+    const groupLen = CURVE.nByteLength;
+    // Normalize bigint / number / string to Uint8Array
+    const keyb =
+      typeof key === 'bigint' || typeof key === 'number'
+        ? ut.numberToBytesLE(normalizeScalar(key, maxGroupElement), groupLen)
+        : key;
+    // Hash private key with curve's hash function to produce uniformingly random input
+    // We check byte lengths e.g.: ensureBytes(64, hash(ensureBytes(32, key)))
+    const hashed = ensureBytes(CURVE.hash(ensureBytes(keyb, groupLen)), 2 * groupLen);
+
+    // First half's bits are cleared to produce a random field element.
+    const head = adjustScalarBytes(hashed.slice(0, groupLen));
+    // Second half is called key prefix (5.1.6)
+    const prefix = hashed.slice(groupLen, 2 * groupLen);
+    // The actual private scalar
+    const scalar = modnLE(head);
+    // Point on Edwards curve aka public key
+    const point = Point.BASE.multiply(scalar);
+    // Uint8Array representation
+    const pointBytes = point.toRawBytes();
+    return { head, prefix, scalar, point, pointBytes };
  }

  /**
@@ -625,7 +599,7 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
  const EMPTY = new Uint8Array();
  function hashDomainToScalar(message: Uint8Array, context: Hex = EMPTY) {
    context = ensureBytes(context);
-    return modlLE(CURVE.hash(domain(message, context, !!CURVE.preHash)));
+    return modnLE(CURVE.hash(domain(message, context, !!CURVE.preHash)));
  }

  /** Signs message with privateKey. RFC8032 5.1.6 */
@@ -633,9 +607,9 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    message = ensureBytes(message);
    if (CURVE.preHash) message = CURVE.preHash(message);
    const { prefix, scalar, pointBytes } = getExtendedPublicKey(privateKey);
-    const r = hashDomainToScalar(concatBytes(prefix, message), context);
+    const r = hashDomainToScalar(ut.concatBytes(prefix, message), context);
    const R = Point.BASE.multiply(r); // R = rG
-    const k = hashDomainToScalar(concatBytes(R.toRawBytes(), pointBytes, message), context); // k = hash(R+P+msg)
+    const k = hashDomainToScalar(ut.concatBytes(R.toRawBytes(), pointBytes, message), context); // k = hash(R+P+msg)
    const s = mod.mod(r + k * scalar, CURVE_ORDER); // s = r + kp
    return new Signature(R, s).toRawBytes();
  }
@@ -672,13 +646,13 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {
    const { r, s } = sig;
    const SB = ExtendedPoint.BASE.multiplyUnsafe(s);
    const k = hashDomainToScalar(
-      concatBytes(r.toRawBytes(), publicKey.toRawBytes(), message),
+      ut.concatBytes(r.toRawBytes(), publicKey.toRawBytes(), message),
      context
    );
    const kA = ExtendedPoint.fromAffine(publicKey).multiplyUnsafe(k);
    const RkA = ExtendedPoint.fromAffine(r).add(kA);
    // [8][S]B = [8]R + [8][k]A'
-    return RkA.subtract(SB).multiplyUnsafe(CURVE.h).equals(ExtendedPoint.ZERO);
+    return RkA.subtract(SB).clearCofactor().equals(ExtendedPoint.ZERO);
  }

  // Enable precomputes. Slows down first publicKey computation by 20ms.
@@ -686,19 +660,16 @@ export function twistedEdwards(curveDef: CurveType): CurveFn {

  const utils = {
    getExtendedPublicKey,
-    mod: modP,
-    invert: Fp.invert,
-
    /**
     * Not needed for ed25519 private keys. Needed if you use scalars directly (rare).
     */
-    hashToPrivateScalar: (hash: Hex): bigint => hashToPrivateScalar(hash, CURVE_ORDER, true),
+    hashToPrivateScalar: (hash: Hex): bigint => ut.hashToPrivateScalar(hash, CURVE_ORDER, true),

    /**
     * ed25519 private keys are uniform 32-bit strings. We do not need to check for
     * modulo bias like we do in secp256k1 randomPrivateKey()
     */
-    randomPrivateKey: (): Uint8Array => randomBytes(fieldLen),
+    randomPrivateKey: (): Uint8Array => randomBytes(Fp.BYTES),

    /**
     * We're doing scalar multiplication (used in getPublicKey etc) with precomputed BASE_POINT
--- a/src/abstract/hash-to-curve.ts
+++ b/src/abstract/hash-to-curve.ts
@@ -17,10 +17,11 @@ export type htfOpts = {
  k: number;
  // option to use a message that has already been processed by
  // expand_message_xmd
-  expand: boolean;
+  expand?: 'xmd' | 'xof';
  // Hash functions for: expand_message_xmd is appropriate for use with a
  // wide range of hash functions, including SHA-2, SHA-3, BLAKE2, and others.
  // BBS+ uses blake2: https://github.com/hyperledger/aries-framework-go/issues/2247
+  // TODO: verify that hash is shake if expand==='xof' via types
  hash: CHash;
 };

@@ -29,7 +30,8 @@ export function validateHTFOpts(opts: htfOpts) {
  if (typeof opts.p !== 'bigint') throw new Error('Invalid htf/p');
  if (typeof opts.m !== 'number') throw new Error('Invalid htf/m');
  if (typeof opts.k !== 'number') throw new Error('Invalid htf/k');
-  if (typeof opts.expand !== 'boolean') throw new Error('Invalid htf/expand');
+  if (opts.expand !== 'xmd' && opts.expand !== 'xof' && opts.expand !== undefined)
+    throw new Error('Invalid htf/expand');
  if (typeof opts.hash !== 'function' || !Number.isSafeInteger(opts.hash.outputLen))
    throw new Error('Invalid htf/hash function');
 }
@@ -101,13 +103,40 @@ export function expand_message_xmd(
  return pseudo_random_bytes.slice(0, lenInBytes);
 }

-// hashes arbitrary-length byte strings to a list of one or more elements of a finite field F
-// https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-5.3
-// Inputs:
-// msg - a byte string containing the message to hash.
-// count - the number of elements of F to output.
-// Outputs:
-// [u_0, ..., u_(count - 1)], a list of field elements.
+export function expand_message_xof(
+  msg: Uint8Array,
+  DST: Uint8Array,
+  lenInBytes: number,
+  k: number,
+  H: CHash
+): Uint8Array {
+  // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-16#section-5.3.3
+  // DST = H('H2C-OVERSIZE-DST-' || a_very_long_DST, Math.ceil((lenInBytes * k) / 8));
+  if (DST.length > 255) {
+    const dkLen = Math.ceil((2 * k) / 8);
+    DST = H.create({ dkLen }).update(stringToBytes('H2C-OVERSIZE-DST-')).update(DST).digest();
+  }
+  if (lenInBytes > 65535 || DST.length > 255)
+    throw new Error('expand_message_xof: invalid lenInBytes');
+  return (
+    H.create({ dkLen: lenInBytes })
+      .update(msg)
+      .update(i2osp(lenInBytes, 2))
+      // 2. DST_prime = DST || I2OSP(len(DST), 1)
+      .update(DST)
+      .update(i2osp(DST.length, 1))
+      .digest()
+  );
+}
+
+/**
+ * Hashes arbitrary-length byte strings to a list of one or more elements of a finite field F
+ * https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-5.3
+ * @param msg a byte string containing the message to hash
+ * @param count the number of elements of F to output
+ * @param options `{DST: string, p: bigint, m: number, k: number, expand: 'xmd' | 'xof', hash: H}`
+ * @returns [u_0, ..., u_(count - 1)], a list of field elements.
+ */
 export function hash_to_field(msg: Uint8Array, count: number, options: htfOpts): bigint[][] {
  // if options is provided but incomplete, fill any missing fields with the
  // value in hftDefaults (ie hash to G2).
@@ -116,8 +145,10 @@ export function hash_to_field(msg: Uint8Array, count: number, options: htfOpts):
  const len_in_bytes = count * options.m * L;
  const DST = stringToBytes(options.DST);
  let pseudo_random_bytes = msg;
-  if (options.expand) {
+  if (options.expand === 'xmd') {
    pseudo_random_bytes = expand_message_xmd(msg, DST, len_in_bytes, options.hash);
+  } else if (options.expand === 'xof') {
+    pseudo_random_bytes = expand_message_xof(msg, DST, len_in_bytes, options.k, options.hash);
  }
  const u = new Array(count);
  for (let i = 0; i < count; i++) {
--- a/src/abstract/modular.ts
+++ b/src/abstract/modular.ts
@@ -1,10 +1,11 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
+// TODO: remove circular imports
 import * as utils from './utils.js';
 // Utilities for modular arithmetics and finite fields
 // prettier-ignore
 const _0n = BigInt(0), _1n = BigInt(1), _2n = BigInt(2), _3n = BigInt(3);
 // prettier-ignore
-const _4n = BigInt(4), _5n = BigInt(5), _7n = BigInt(7), _8n = BigInt(8);
+const _4n = BigInt(4), _5n = BigInt(5), _8n = BigInt(8);
 // prettier-ignore
 const _9n = BigInt(9), _16n = BigInt(16);

@@ -54,6 +55,7 @@ export function invert(number: bigint, modulo: bigint): bigint {
  // prettier-ignore
  let x = _0n, y = _1n, u = _1n, v = _0n;
  while (a !== _0n) {
+    // JIT applies optimization if those two lines follow each other
    const q = b / a;
    const r = b % a;
    const m = x - u * q;
@@ -66,26 +68,68 @@ export function invert(number: bigint, modulo: bigint): bigint {
  return mod(x, modulo);
 }

-/**
- * Calculates Legendre symbol (a | p), which denotes the value of a^((p-1)/2) (mod p).
- * * (a | p) ≡ 1    if a is a square (mod p)
- * * (a | p) ≡ -1   if a is not a square (mod p)
- * * (a | p) ≡ 0    if a ≡ 0 (mod p)
- */
-export function legendre(num: bigint, fieldPrime: bigint): bigint {
-  return pow(num, (fieldPrime - _1n) / _2n, fieldPrime);
+// Tonelli-Shanks algorithm
+// Paper 1: https://eprint.iacr.org/2012/685.pdf (page 12)
+// Paper 2: Square Roots from 1; 24, 51, 10 to Dan Shanks
+export function tonelliShanks(P: bigint) {
+  // Legendre constant: used to calculate Legendre symbol (a | p),
+  // which denotes the value of a^((p-1)/2) (mod p).
+  // (a | p) ≡ 1    if a is a square (mod p)
+  // (a | p) ≡ -1   if a is not a square (mod p)
+  // (a | p) ≡ 0    if a ≡ 0 (mod p)
+  const legendreC = (P - _1n) / _2n;
+
+  let Q: bigint, S: number, Z: bigint;
+  // Step 1: By factoring out powers of 2 from p - 1,
+  // find q and s such that p - 1 = q*(2^s) with q odd
+  for (Q = P - _1n, S = 0; Q % _2n === _0n; Q /= _2n, S++);
+
+  // Step 2: Select a non-square z such that (z | p) ≡ -1 and set c ≡ zq
+  for (Z = _2n; Z < P && pow(Z, legendreC, P) !== P - _1n; Z++);
+
+  // Fast-path
+  if (S === 1) {
+    const p1div4 = (P + _1n) / _4n;
+    return function tonelliFast<T>(Fp: Field<T>, n: T) {
+      const root = Fp.pow(n, p1div4);
+      if (!Fp.equals(Fp.square(root), n)) throw new Error('Cannot find square root');
+      return root;
+    };
  }

-/**
- * Calculates square root of a number in a finite field.
- * √a mod P
- */
-// TODO: rewrite as generic Fp function && remove bls versions
-export function sqrt(number: bigint, modulo: bigint): bigint {
-  // prettier-ignore
-  const n = number;
-  const P = modulo;
-  const p1div4 = (P + _1n) / _4n;
+  // Slow-path
+  const Q1div2 = (Q + _1n) / _2n;
+  return function tonelliSlow<T>(Fp: Field<T>, n: T): T {
+    // Step 0: Check that n is indeed a square: (n | p) should not be ≡ -1
+    if (Fp.pow(n, legendreC) === Fp.negate(Fp.ONE)) throw new Error('Cannot find square root');
+    let r = S;
+    // TODO: will fail at Fp2/etc
+    let g = Fp.pow(Fp.mul(Fp.ONE, Z), Q); // will update both x and b
+    let x = Fp.pow(n, Q1div2); // first guess at the square root
+    let b = Fp.pow(n, Q); // first guess at the fudge factor
+
+    while (!Fp.equals(b, Fp.ONE)) {
+      if (Fp.equals(b, Fp.ZERO)) return Fp.ZERO; // https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm (4. If t = 0, return r = 0)
+      // Find m such b^(2^m)==1
+      let m = 1;
+      for (let t2 = Fp.square(b); m < r; m++) {
+        if (Fp.equals(t2, Fp.ONE)) break;
+        t2 = Fp.square(t2); // t2 *= t2
+      }
+      // NOTE: r-m-1 can be bigger than 32, need to convert to bigint before shift, otherwise there will be overflow
+      const ge = Fp.pow(g, _1n << BigInt(r - m - 1)); // ge = 2^(r-m-1)
+      g = Fp.square(ge); // g = ge * ge
+      x = Fp.mul(x, ge); // x *= ge
+      b = Fp.mul(b, g); // b *= g
+      r = m;
+    }
+    return x;
+  };
+}
+
+export function FpSqrt(P: bigint) {
+  // NOTE: different algorithms can give different roots, it is up to user to decide which one they want.
+  // For example there is FpSqrtOdd/FpSqrtEven to choice root based on oddness (used for hash-to-curve).

  // P ≡ 3 (mod 4)
  // √n = n^((P+1)/4)
@@ -94,48 +138,54 @@ export function sqrt(number: bigint, modulo: bigint): bigint {
    // const ORDER =
    //   0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaabn;
    // const NUM = 72057594037927816n;
-    // TODO: fix sqrtMod in secp256k1
-    const root = pow(n, p1div4, P);
-    if (mod(root * root, modulo) !== number) throw new Error('Cannot find square root');
+    const p1div4 = (P + _1n) / _4n;
+    return function sqrt3mod4<T>(Fp: Field<T>, n: T) {
+      const root = Fp.pow(n, p1div4);
+      // Throw if root**2 != n
+      if (!Fp.equals(Fp.square(root), n)) throw new Error('Cannot find square root');
      return root;
+    };
  }

-  // P ≡ 5 (mod 8)
+  // Atkin algorithm for q ≡ 5 (mod 8), https://eprint.iacr.org/2012/685.pdf (page 10)
  if (P % _8n === _5n) {
-    const n2 = mod(n * _2n, P);
-    const v = pow(n2, (P - _5n) / _8n, P);
-    const nv = mod(n * v, P);
-    const i = mod(_2n * nv * v, P);
-    const r = mod(nv * (i - _1n), P);
-    return r;
+    const c1 = (P - _5n) / _8n;
+    return function sqrt5mod8<T>(Fp: Field<T>, n: T) {
+      const n2 = Fp.mul(n, _2n);
+      const v = Fp.pow(n2, c1);
+      const nv = Fp.mul(n, v);
+      const i = Fp.mul(Fp.mul(nv, _2n), v);
+      const root = Fp.mul(nv, Fp.sub(i, Fp.ONE));
+      if (!Fp.equals(Fp.square(root), n)) throw new Error('Cannot find square root');
+      return root;
+    };
+  }
+
+  // P ≡ 9 (mod 16)
+  if (P % _16n === _9n) {
+    // NOTE: tonelli is too slow for bls-Fp2 calculations even on start
+    // Means we cannot use sqrt for constants at all!
+    //
+    // const c1 = Fp.sqrt(Fp.negate(Fp.ONE)); //  1. c1 = sqrt(-1) in F, i.e., (c1^2) == -1 in F
+    // const c2 = Fp.sqrt(c1);                //  2. c2 = sqrt(c1) in F, i.e., (c2^2) == c1 in F
+    // const c3 = Fp.sqrt(Fp.negate(c1));     //  3. c3 = sqrt(-c1) in F, i.e., (c3^2) == -c1 in F
+    // const c4 = (P + _7n) / _16n;           //  4. c4 = (q + 7) / 16        # Integer arithmetic
+    // sqrt = (x) => {
+    //   let tv1 = Fp.pow(x, c4);             //  1. tv1 = x^c4
+    //   let tv2 = Fp.mul(c1, tv1);           //  2. tv2 = c1 * tv1
+    //   const tv3 = Fp.mul(c2, tv1);         //  3. tv3 = c2 * tv1
+    //   let tv4 = Fp.mul(c3, tv1);           //  4. tv4 = c3 * tv1
+    //   const e1 = Fp.equals(Fp.square(tv2), x); //  5.  e1 = (tv2^2) == x
+    //   const e2 = Fp.equals(Fp.square(tv3), x); //  6.  e2 = (tv3^2) == x
+    //   tv1 = Fp.cmov(tv1, tv2, e1); //  7. tv1 = CMOV(tv1, tv2, e1)  # Select tv2 if (tv2^2) == x
+    //   tv2 = Fp.cmov(tv4, tv3, e2); //  8. tv2 = CMOV(tv4, tv3, e2)  # Select tv3 if (tv3^2) == x
+    //   const e3 = Fp.equals(Fp.square(tv2), x); //  9.  e3 = (tv2^2) == x
+    //   return Fp.cmov(tv1, tv2, e3); //  10.  z = CMOV(tv1, tv2, e3)  # Select the sqrt from tv1 and tv2
+    // }
  }

  // Other cases: Tonelli-Shanks algorithm
-  if (legendre(n, P) !== _1n) throw new Error('Cannot find square root');
-  let q: bigint, s: number, z: bigint;
-  for (q = P - _1n, s = 0; q % _2n === _0n; q /= _2n, s++);
-  if (s === 1) return pow(n, p1div4, P);
-  for (z = _2n; z < P && legendre(z, P) !== P - _1n; z++);
-
-  let c = pow(z, q, P);
-  let r = pow(n, (q + _1n) / _2n, P);
-  let t = pow(n, q, P);
-
-  let t2 = _0n;
-  while (mod(t - _1n, P) !== _0n) {
-    t2 = mod(t * t, P);
-    let i;
-    for (i = 1; i < s; i++) {
-      if (mod(t2 - _1n, P) === _0n) break;
-      t2 = mod(t2 * t2, P);
-    }
-    let b = pow(c, BigInt(1 << (s - i - 1)), P);
-    r = mod(r * b, P);
-    c = mod(b * b, P);
-    t = mod(t * c, P);
-    s = i;
-  }
-  return r;
+  return tonelliShanks(P);
 }

 // Little-endian check for first LE bit (last BE bit);
@@ -176,6 +226,7 @@ export interface Field<T> {

  // Optional
  // Should be same as sgn0 function in https://datatracker.ietf.org/doc/draft-irtf-cfrg-hash-to-curve/
+  // NOTE: sgn0 is 'negative in LE', which is same as odd. And negative in LE is kinda strange definition anyway.
  isOdd?(num: T): boolean; // Odd instead of even since we have it for Fp2
  legendre?(num: T): T;
  pow(lhs: T, power: bigint): T;
@@ -246,21 +297,31 @@ export function FpDiv<T>(f: Field<T>, lhs: T, rhs: T | bigint): T {
  return f.mul(lhs, typeof rhs === 'bigint' ? invert(rhs, f.ORDER) : f.invert(rhs));
 }

+// This function returns True whenever the value x is a square in the field F.
+export function FpIsSquare<T>(f: Field<T>) {
+  const legendreConst = (f.ORDER - _1n) / _2n; // Integer arithmetic
+  return (x: T): boolean => {
+    const p = f.pow(x, legendreConst);
+    return f.equals(p, f.ZERO) || f.equals(p, f.ONE);
+  };
+}
+
 // NOTE: very fragile, always bench. Major performance points:
 // - NonNormalized ops
 // - Object.freeze
 // - same shape of object (don't add/remove keys)
+type FpField = Field<bigint> & Required<Pick<Field<bigint>, 'isOdd'>>;
 export function Fp(
  ORDER: bigint,
  bitLen?: number,
  isLE = false,
  redef: Partial<Field<bigint>> = {}
-): Readonly<Field<bigint>> {
+): Readonly<FpField> {
  if (ORDER <= _0n) throw new Error(`Expected Fp ORDER > 0, got ${ORDER}`);
  const { nBitLength: BITS, nByteLength: BYTES } = utils.nLength(ORDER, bitLen);
  if (BYTES > 2048) throw new Error('Field lengths over 2048 bytes are not supported');
-  const sqrtP = (num: bigint) => sqrt(num, ORDER);
-  const f: Field<bigint> = Object.freeze({
+  const sqrtP = FpSqrt(ORDER);
+  const f: Readonly<FpField> = Object.freeze({
    ORDER,
    BITS,
    BYTES,
@@ -292,7 +353,7 @@ export function Fp(
    mulN: (lhs, rhs) => lhs * rhs,

    invert: (num) => invert(num, ORDER),
-    sqrt: redef.sqrt || sqrtP,
+    sqrt: redef.sqrt || ((n) => sqrtP(f, n)),
    invertBatch: (lst) => FpInvertBatch(f, lst),
    // TODO: do we really need constant cmov?
    // We don't have const-time bigints anyway, so probably will be not very useful
@@ -305,87 +366,18 @@ export function Fp(
        throw new Error(`Fp.fromBytes: expected ${BYTES}, got ${bytes.length}`);
      return isLE ? utils.bytesToNumberLE(bytes) : utils.bytesToNumberBE(bytes);
    },
-  } as Field<bigint>);
+  } as FpField);
  return Object.freeze(f);
 }

-// TODO: re-use in bls/generic sqrt for field/etc?
-// Something like sqrtUnsafe which always returns value, but sqrt throws exception if non-square
-// From draft-irtf-cfrg-hash-to-curve-16
-export function FpSqrt<T>(Fp: Field<T>) {
-  // NOTE: it requires another sqrt for constant precomputes, but no need for roots of unity,
-  // probably we can simply bls code using it
-  const q = Fp.ORDER;
-  const squareConst = (q - _1n) / _2n;
-  // is_square(x) := { True,  if x^((q - 1) / 2) is 0 or 1 in F;
-  //                 { False, otherwise.
-  let isSquare: (x: T) => boolean = (x) => {
-    const p = Fp.pow(x, squareConst);
-    return Fp.equals(p, Fp.ZERO) || Fp.equals(p, Fp.ONE);
-  };
-  // Constant-time Tonelli-Shanks algorithm
-  let l = _0n;
-  for (let o = q - _1n; o % _2n === _0n; o /= _2n) l += _1n;
-  const c1 = l; // 1. c1, the largest integer such that 2^c1 divides q - 1.
-  const c2 = (q - _1n) / _2n ** c1; // 2. c2 = (q - 1) / (2^c1)        # Integer arithmetic
-  const c3 = (c2 - _1n) / _2n; // 3. c3 = (c2 - 1) / 2            # Integer arithmetic
-  //  4. c4, a non-square value in F
-  //  5. c5 = c4^c2 in F
-  let c4 = Fp.ONE;
-  while (isSquare(c4)) c4 = Fp.add(c4, Fp.ONE);
-  const c5 = Fp.pow(c4, c2);
+export function FpSqrtOdd<T>(Fp: Field<T>, elm: T) {
+  if (!Fp.isOdd) throw new Error(`Field doesn't have isOdd`);
+  const root = Fp.sqrt(elm);
+  return Fp.isOdd(root) ? root : Fp.negate(root);
+}

-  let sqrt: (x: T) => T = (x) => {
-    let z = Fp.pow(x, c3); // 1.  z = x^c3
-    let t = Fp.square(z); // 2.  t = z * z
-    t = Fp.mul(t, x); // 3.  t = t * x
-    z = Fp.mul(z, x); // 4.  z = z * x
-    let b = t; // 5.  b = t
-    let c = c5; // 6.  c = c5
-    // 7.  for i in (c1, c1 - 1, ..., 2):
-    for (let i = c1; i > 1; i--) {
-      // 8.    for j in (1, 2, ..., i - 2):
-      // 9.      b = b * b
-      for (let j = _1n; j < i - _1n; i++) b = Fp.square(b);
-      const e = Fp.equals(b, Fp.ONE); // 10.   e = b == 1
-      const zt = Fp.mul(z, c); // 11.   zt = z * c
-      z = Fp.cmov(zt, z, e); // 12.   z = CMOV(zt, z, e)
-      c = Fp.square(c); // 13.   c = c * c
-      let tt = Fp.mul(t, c); // 14.   tt = t * c
-      t = Fp.cmov(tt, t, e); // 15.   t = CMOV(tt, t, e)
-      b = t; // 16.   b = t
-    }
-    return z; // 17. return z
-  };
-  if (q % _4n === _3n) {
-    const c1 = (q + _1n) / _4n; // 1. c1 = (q + 1) / 4     # Integer arithmetic
-    sqrt = (x) => Fp.pow(x, c1);
-  } else if (q % _8n === _5n) {
-    const c1 = Fp.sqrt(Fp.negate(Fp.ONE)); //  1. c1 = sqrt(-1) in F, i.e., (c1^2) == -1 in F
-    const c2 = (q + _3n) / _8n; //  2. c2 = (q + 3) / 8     # Integer arithmetic
-    sqrt = (x) => {
-      let tv1 = Fp.pow(x, c2); //  1. tv1 = x^c2
-      let tv2 = Fp.mul(tv1, c1); //  2. tv2 = tv1 * c1
-      let e = Fp.equals(Fp.square(tv1), x); //  3.   e = (tv1^2) == x
-      return Fp.cmov(tv2, tv1, e); //  4.   z = CMOV(tv2, tv1, e)
-    };
-  } else if (Fp.ORDER % _16n === _9n) {
-    const c1 = Fp.sqrt(Fp.negate(Fp.ONE)); //  1. c1 = sqrt(-1) in F, i.e., (c1^2) == -1 in F
-    const c2 = Fp.sqrt(c1); //  2. c2 = sqrt(c1) in F, i.e., (c2^2) == c1 in F
-    const c3 = Fp.sqrt(Fp.negate(c1)); //  3. c3 = sqrt(-c1) in F, i.e., (c3^2) == -c1 in F
-    const c4 = (Fp.ORDER + _7n) / _16n; //  4. c4 = (q + 7) / 16         # Integer arithmetic
-    sqrt = (x) => {
-      let tv1 = Fp.pow(x, c4); //  1. tv1 = x^c4
-      let tv2 = Fp.mul(c1, tv1); //  2. tv2 = c1 * tv1
-      const tv3 = Fp.mul(c2, tv1); //  3. tv3 = c2 * tv1
-      let tv4 = Fp.mul(c3, tv1); //  4. tv4 = c3 * tv1
-      const e1 = Fp.equals(Fp.square(tv2), x); //  5.  e1 = (tv2^2) == x
-      const e2 = Fp.equals(Fp.square(tv3), x); //  6.  e2 = (tv3^2) == x
-      tv1 = Fp.cmov(tv1, tv2, e1); //  7. tv1 = CMOV(tv1, tv2, e1)  # Select tv2 if (tv2^2) == x
-      tv2 = Fp.cmov(tv4, tv3, e2); //  8. tv2 = CMOV(tv4, tv3, e2)  # Select tv3 if (tv3^2) == x
-      const e3 = Fp.equals(Fp.square(tv2), x); //  9.  e3 = (tv2^2) == x
-      return Fp.cmov(tv1, tv2, e3); //  10.  z = CMOV(tv1, tv2, e3)  # Select the sqrt from tv1 and tv2
-    };
-  }
-  return { sqrt, isSquare };
+export function FpSqrtEven<T>(Fp: Field<T>, elm: T) {
+  if (!Fp.isOdd) throw new Error(`Field doesn't have isOdd`);
+  const root = Fp.sqrt(elm);
+  return Fp.isOdd(root) ? Fp.negate(root) : root;
 }
--- a/src/abstract/montgomery.ts
+++ b/src/abstract/montgomery.ts
@@ -1,11 +1,6 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
 import * as mod from './modular.js';
-import {
-  ensureBytes,
-  numberToBytesLE,
-  bytesToNumberLE,
-  // nLength,
-} from './utils.js';
+import { ensureBytes, numberToBytesLE, bytesToNumberLE, isPositiveInt } from './utils.js';

 const _0n = BigInt(0);
 const _1n = BigInt(1);
@@ -38,7 +33,7 @@ function validateOpts(curve: CurveType) {
  }
  for (const i of ['montgomeryBits', 'nByteLength'] as const) {
    if (curve[i] === undefined) continue; // Optional
-    if (!Number.isSafeInteger(curve[i]))
+    if (!isPositiveInt(curve[i]))
      throw new Error(`Invalid curve param ${i}=${curve[i]} (${typeof curve[i]})`);
  }
  for (const fn of ['adjustScalarBytes', 'domain', 'powPminus2'] as const) {
--- a/src/abstract/utils.ts
+++ b/src/abstract/utils.ts
@@ -12,7 +12,7 @@ export type CHash = {
  (message: Uint8Array | string): Uint8Array;
  blockLen: number;
  outputLen: number;
-  create(): any;
+  create(opts?: { dkLen?: number }): any; // For shake
 };

 // NOTE: these are generic, even if curve is on some polynominal field (bls), it will still have P/n/h
@@ -40,19 +40,24 @@ export type BasicCurve<T> = {
  allowInfinityPoint?: boolean;
 };

+// Bans floats and integers above 2^53-1
+export function isPositiveInt(num: any): num is number {
+  return typeof num === 'number' && Number.isSafeInteger(num) && num > 0;
+}
+
 export function validateOpts<FP, T>(curve: BasicCurve<FP> & T) {
  mod.validateField(curve.Fp);
  for (const i of ['n', 'h'] as const) {
-    if (typeof curve[i] !== 'bigint')
-      throw new Error(`Invalid curve param ${i}=${curve[i]} (${typeof curve[i]})`);
+    const val = curve[i];
+    if (typeof val !== 'bigint') throw new Error(`Invalid curve param ${i}=${val} (${typeof val})`);
  }
  if (!curve.Fp.isValid(curve.Gx)) throw new Error('Invalid generator X coordinate Fp element');
  if (!curve.Fp.isValid(curve.Gy)) throw new Error('Invalid generator Y coordinate Fp element');

  for (const i of ['nBitLength', 'nByteLength'] as const) {
-    if (curve[i] === undefined) continue; // Optional
-    if (!Number.isSafeInteger(curve[i]))
-      throw new Error(`Invalid curve param ${i}=${curve[i]} (${typeof curve[i]})`);
+    const val = curve[i];
+    if (val === undefined) continue; // Optional
+    if (!isPositiveInt(val)) throw new Error(`Invalid curve param ${i}=${val} (${typeof val})`);
  }
  // Set defaults
  return Object.freeze({ ...nLength(curve.n, curve.nBitLength), ...curve } as const);
@@ -144,21 +149,22 @@ export function nLength(n: bigint, nBitLength?: number) {
 }

 /**
+ * FIPS 186 B.4.1-compliant "constant-time" private key generation utility.
 * Can take (n+8) or more bytes of uniform input e.g. from CSPRNG or KDF
 * and convert them into private scalar, with the modulo bias being neglible.
- * As per FIPS 186 B.4.1.
+ * Needs at least 40 bytes of input for 32-byte private key.
 * https://research.kudelskisecurity.com/2020/07/28/the-definitive-guide-to-modulo-bias-and-how-to-avoid-it/
- * @param hash hash output from sha512, or a similar function
+ * @param hash hash output from SHA3 or a similar function
 * @returns valid private scalar
 */
-export function hashToPrivateScalar(hash: Hex, CURVE_ORDER: bigint, isLE = false): bigint {
+export function hashToPrivateScalar(hash: Hex, groupOrder: bigint, isLE = false): bigint {
  hash = ensureBytes(hash);
-  const orderLen = nLength(CURVE_ORDER).nByteLength;
-  const minLen = orderLen + 8;
-  if (orderLen < 16 || hash.length < minLen || hash.length > 1024)
-    throw new Error('Expected valid bytes of private key as per FIPS 186');
+  const hashLen = hash.length;
+  const minLen = nLength(groupOrder).nByteLength + 8;
+  if (minLen < 24 || hashLen < minLen || hashLen > 1024)
+    throw new Error(`hashToPrivateScalar: expected ${minLen}-1024 bytes of input, got ${hashLen}`);
  const num = isLE ? bytesToNumberLE(hash) : bytesToNumberBE(hash);
-  return mod.mod(num, CURVE_ORDER - _1n) + _1n;
+  return mod.mod(num, groupOrder - _1n) + _1n;
 }

 export function equalBytes(b1: Uint8Array, b2: Uint8Array) {
--- a/src/abstract/weierstrass.ts
+++ b/src/abstract/weierstrass.ts
@@ -1,28 +1,16 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
 // Short Weierstrass curve. The formula is: y² = x³ + ax + b

-// TODO: sync vs async naming
-// TODO: default randomBytes
 // Differences from @noble/secp256k1 1.7:
 // 1. Different double() formula (but same addition)
 // 2. Different sqrt() function
 // 3. truncateHash() truncateOnly mode
 // 4. DRBG supports outputLen bigger than outputLen of hmac
+// 5. Support for different hash functions

 import * as mod from './modular.js';
-import {
-  bytesToHex,
-  bytesToNumberBE,
-  concatBytes,
-  ensureBytes,
-  hexToBytes,
-  hexToNumber,
-  numberToHexUnpadded,
-  hashToPrivateScalar,
-  Hex,
-  PrivKey,
-} from './utils.js';
-import * as utils from './utils.js';
+import * as ut from './utils.js';
+import { bytesToHex, Hex, PrivKey } from './utils.js';
 import { hash_to_field, htfOpts, validateHTFOpts } from './hash-to-curve.js';
 import { Group, GroupConstructor, wNAF } from './group.js';

@@ -31,39 +19,45 @@ type EndomorphismOpts = {
  beta: bigint;
  splitScalar: (k: bigint) => { k1neg: boolean; k1: bigint; k2neg: boolean; k2: bigint };
 };
-export type BasicCurve<T> = utils.BasicCurve<T> & {
+export type BasicCurve<T> = ut.BasicCurve<T> & {
  // Params: a, b
  a: T;
  b: T;
-  // TODO: move into options?

+  // Optional params
+  // Executed before privkey validation. Useful for P521 with var-length priv key
  normalizePrivateKey?: (key: PrivKey) => PrivKey;
+  // Whether to execute modular division on a private key, useful for bls curves with cofactor > 1
+  wrapPrivateKey?: boolean;
  // Endomorphism options for Koblitz curves
  endo?: EndomorphismOpts;
-  // Torsions, can be optimized via endomorphisms
+  // When a cofactor != 1, there can be an effective methods to:
+  // 1. Determine whether a point is torsion-free
  isTorsionFree?: (c: ProjectiveConstructor<T>, point: ProjectivePointType<T>) => boolean;
+  // 2. Clear torsion component
  clearCofactor?: (
    c: ProjectiveConstructor<T>,
    point: ProjectivePointType<T>
  ) => ProjectivePointType<T>;
-  // Hash to field opts
+  // Hash to field options
  htfDefaults?: htfOpts;
  mapToCurve?: (scalar: bigint[]) => { x: T; y: T };
 };
-// DER encoding utilities
+
+// ASN.1 DER encoding utilities
 class DERError extends Error {
  constructor(message: string) {
    super(message);
  }
 }

-function sliceDER(s: string): string {
+const DER = {
+  slice(s: string): string {
    // Proof: any([(i>=0x80) == (int(hex(i).replace('0x', '').zfill(2)[0], 16)>=8)  for i in range(0, 256)])
    // Padding done by numberToHex
    return Number.parseInt(s[0], 16) >= 8 ? '00' + s : s;
-}
-
-function parseDERInt(data: Uint8Array) {
+  },
+  parseInt(data: Uint8Array): { data: bigint; left: Uint8Array } {
    if (data.length < 2 || data[0] !== 0x02) {
      throw new DERError(`Invalid signature integer tag: ${bytesToHex(data)}`);
    }
@@ -76,31 +70,26 @@ function parseDERInt(data: Uint8Array) {
    if (res[0] === 0x00 && res[1] <= 0x7f) {
      throw new DERError('Invalid signature integer: trailing length');
    }
-  return { data: bytesToNumberBE(res), left: data.subarray(len + 2) };
-}
-
-function parseDERSignature(data: Uint8Array) {
+    return { data: ut.bytesToNumberBE(res), left: data.subarray(len + 2) };
+  },
+  parseSig(data: Uint8Array): { r: bigint; s: bigint } {
    if (data.length < 2 || data[0] != 0x30) {
      throw new DERError(`Invalid signature tag: ${bytesToHex(data)}`);
    }
    if (data[1] !== data.length - 2) {
      throw new DERError('Invalid signature: incorrect length');
    }
-  const { data: r, left: sBytes } = parseDERInt(data.subarray(2));
-  const { data: s, left: rBytesLeft } = parseDERInt(sBytes);
+    const { data: r, left: sBytes } = DER.parseInt(data.subarray(2));
+    const { data: s, left: rBytesLeft } = DER.parseInt(sBytes);
    if (rBytesLeft.length) {
      throw new DERError(`Invalid signature: left bytes after parsing: ${bytesToHex(rBytesLeft)}`);
    }
    return { r, s };
-}
-
-// Be friendly to bad ECMAScript parsers by not using bigint literals like 123n
-const _0n = BigInt(0);
-const _1n = BigInt(1);
-const _3n = BigInt(3);
+  },
+};

 type Entropy = Hex | true;
-type SignOpts = { lowS?: boolean; extraEntropy?: Entropy };
+export type SignOpts = { lowS?: boolean; extraEntropy?: Entropy };

 /**
 * ### Design rationale for types
@@ -124,7 +113,7 @@ type SignOpts = { lowS?: boolean; extraEntropy?: Entropy };
 * TODO: https://www.typescriptlang.org/docs/handbook/release-notes/typescript-2-7.html#unique-symbol
 */

-// Instance
+// Instance for 3d XYZ points
 export interface ProjectivePointType<T> extends Group<ProjectivePointType<T>> {
  readonly x: T;
  readonly y: T;
@@ -133,14 +122,14 @@ export interface ProjectivePointType<T> extends Group<ProjectivePointType<T>> {
  multiplyUnsafe(scalar: bigint): ProjectivePointType<T>;
  toAffine(invZ?: T): PointType<T>;
 }
-// Static methods
+// Static methods for 3d XYZ points
 export interface ProjectiveConstructor<T> extends GroupConstructor<ProjectivePointType<T>> {
  new (x: T, y: T, z: T): ProjectivePointType<T>;
  fromAffine(p: PointType<T>): ProjectivePointType<T>;
  toAffineBatch(points: ProjectivePointType<T>[]): PointType<T>[];
  normalizeZ(points: ProjectivePointType<T>[]): ProjectivePointType<T>[];
 }
-// Instance
+// Instance for 2d XY points
 export interface PointType<T> extends Group<PointType<T>> {
  readonly x: T;
  readonly y: T;
@@ -151,7 +140,7 @@ export interface PointType<T> extends Group<PointType<T>> {
  assertValidity(): void;
  multiplyAndAddUnsafe(Q: PointType<T>, a: bigint, b: bigint): PointType<T> | undefined;
 }
-// Static methods
+// Static methods for 2d XY points
 export interface PointConstructor<T> extends GroupConstructor<PointType<T>> {
  new (x: T, y: T): PointType<T>;
  fromHex(hex: Hex): PointType<T>;
@@ -167,7 +156,7 @@ export type CurvePointsType<T> = BasicCurve<T> & {
 };

 function validatePointOpts<T>(curve: CurvePointsType<T>) {
-  const opts = utils.validateOpts(curve);
+  const opts = ut.validateOpts(curve);
  const Fp = opts.Fp;
  for (const i of ['a', 'b'] as const) {
    if (!Fp.isValid(curve[i]))
@@ -206,22 +195,14 @@ export type CurvePointsRes<T> = {
  isWithinCurveOrder: (num: bigint) => boolean;
 };

+// Be friendly to bad ECMAScript parsers by not using bigint literals like 123n
+const _0n = BigInt(0);
+const _1n = BigInt(1);
+const _3n = BigInt(3);
+
 export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
  const CURVE = validatePointOpts(opts);
-  const Fp = CURVE.Fp;
-  // Lengths
-  // All curves has same field / group length as for now, but it can be different for other curves
-  const { nByteLength, nBitLength } = CURVE;
-  const groupLen = nByteLength;
-
-  // Not using ** operator with bigints for old engines.
-  // 2n ** (8n * 32n) == 2n << (8n * 32n - 1n)
-  //const FIELD_MASK = _2n << (_8n * BigInt(fieldLen) - _1n);
-  // function numToFieldStr(num: bigint): string {
-  //   if (typeof num !== 'bigint') throw new Error('Expected bigint');
-  //   if (!(_0n <= num && num < FIELD_MASK)) throw new Error(`Expected number < 2^${fieldLen * 8}`);
-  //   return num.toString(16).padStart(2 * fieldLen, '0');
-  // }
+  const { Fp } = CURVE; // All curves has same field / group length as for now, but they can differ

  /**
   * y² = x³ + ax + b: Short weierstrass curve formula
@@ -234,35 +215,48 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
    return Fp.add(Fp.add(x3, Fp.mul(x, a)), b); // x3 + a * x + b
  }

+  // Valid group elements reside in range 1..n-1
  function isWithinCurveOrder(num: bigint): boolean {
    return _0n < num && num < CURVE.n;
  }

+  /**
+   * Validates if a private key is valid and converts it to bigint form.
+   * Supports two options, that are passed when CURVE is initialized:
+   * - `normalizePrivateKey()` executed before all checks
+   * - `wrapPrivateKey` when true, executed after most checks, but before `0 < key < n`
+   */
  function normalizePrivateKey(key: PrivKey): bigint {
-    if (typeof CURVE.normalizePrivateKey === 'function') {
-      key = CURVE.normalizePrivateKey(key);
-    }
+    const { normalizePrivateKey: custom, nByteLength: groupLen, wrapPrivateKey, n: order } = CURVE;
+    if (typeof custom === 'function') key = custom(key);
    let num: bigint;
    if (typeof key === 'bigint') {
+      // Curve order check is done below
      num = key;
-    } else if (typeof key === 'number' && Number.isSafeInteger(key) && key > 0) {
+    } else if (ut.isPositiveInt(key)) {
      num = BigInt(key);
    } else if (typeof key === 'string') {
      if (key.length !== 2 * groupLen) throw new Error(`Expected ${groupLen} bytes of private key`);
-      num = hexToNumber(key);
+      // Validates individual octets
+      num = ut.hexToNumber(key);
    } else if (key instanceof Uint8Array) {
      if (key.length !== groupLen) throw new Error(`Expected ${groupLen} bytes of private key`);
-      num = bytesToNumberBE(key);
+      num = ut.bytesToNumberBE(key);
    } else {
      throw new TypeError('Expected valid private key');
    }
-    if (CURVE.wrapPrivateKey) num = mod.mod(num, CURVE.n);
+    // Useful for curves with cofactor != 1
+    if (wrapPrivateKey) num = mod.mod(num, order);
    if (!isWithinCurveOrder(num)) throw new Error('Expected private key: 0 < key < n');
    return num;
  }

+  /**
+   * Validates if a scalar ("private number") is valid.
+   * Scalars are valid only if they are less than curve order.
+   */
  function normalizeScalar(num: number | bigint): bigint {
-    if (typeof num === 'number' && Number.isSafeInteger(num) && num > 0) return BigInt(num);
+    if (ut.isPositiveInt(num)) return BigInt(num);
    if (typeof num === 'bigint' && isWithinCurveOrder(num)) return num;
    throw new TypeError('Expected valid private scalar: 0 < scalar < curve.n');
  }
@@ -289,7 +283,7 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {

    /**
     * Takes a bunch of Projective Points but executes only one
-     * invert on all of them. invert is very slow operation,
+     * inversion on all of them. Inversion is very slow operation,
     * so this improves performance massively.
     */
    static toAffineBatch(points: ProjectivePoint[]): Point[] {
@@ -297,6 +291,9 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
      return points.map((p, i) => p.toAffine(toInv[i]));
    }

+    /**
+     * Optimization: converts a list of projective points to a list of identical points with Z=1.
+     */
    static normalizeZ(points: ProjectivePoint[]): ProjectivePoint[] {
      return ProjectivePoint.toAffineBatch(points).map(ProjectivePoint.fromAffine);
    }
@@ -320,10 +317,6 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
      return new ProjectivePoint(this.x, Fp.negate(this.y), this.z);
    }

-    doubleAdd(): ProjectivePoint {
-      return this.add(this);
-    }
-
    // Renes-Costello-Batina exception-free doubling formula.
    // There is 30% faster Jacobian formula, but it is not complete.
    // https://eprint.iacr.org/2015/1060, algorithm 3
@@ -525,24 +518,25 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
      return new Point(ax, ay);
    }
    isTorsionFree(): boolean {
-      if (CURVE.h === _1n) return true; // No subgroups, always torsion fee
-      if (CURVE.isTorsionFree) return CURVE.isTorsionFree(ProjectivePoint, this);
-      // is multiplyUnsafe(CURVE.n) is always ok, same as for edwards?
-      throw new Error('Unsupported!');
+      const { h: cofactor, isTorsionFree } = CURVE;
+      if (cofactor === _1n) return true; // No subgroups, always torsion-free
+      if (isTorsionFree) return isTorsionFree(ProjectivePoint, this);
+      throw new Error('isTorsionFree() has not been declared for the elliptic curve');
    }
-    // Clear cofactor of G1
-    // https://eprint.iacr.org/2019/403
    clearCofactor(): ProjectivePoint {
-      if (CURVE.h === _1n) return this; // Fast-path
-      if (CURVE.clearCofactor) return CURVE.clearCofactor(ProjectivePoint, this) as ProjectivePoint;
+      const { h: cofactor, clearCofactor } = CURVE;
+      if (cofactor === _1n) return this; // Fast-path
+      if (clearCofactor) return clearCofactor(ProjectivePoint, this) as ProjectivePoint;
      return this.multiplyUnsafe(CURVE.h);
    }
  }
-  const wnaf = wNAF(ProjectivePoint, CURVE.endo ? nBitLength / 2 : nBitLength);
+  const _bits = CURVE.nBitLength;
+  const wnaf = wNAF(ProjectivePoint, CURVE.endo ? Math.ceil(_bits / 2) : _bits);

  function assertPrjPoint(other: unknown) {
    if (!(other instanceof ProjectivePoint)) throw new TypeError('ProjectivePoint expected');
  }
+
  // Stores precomputed values for points.
  const pointPrecomputes = new WeakMap<Point, ProjectivePoint[]>();

@@ -551,11 +545,11 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
   */
  class Point implements PointType<T> {
    /**
-     * Base point aka generator. public_key = Point.BASE * private_key
+     * Base point aka generator. Any public_key = Point.BASE * private_key
     */
    static BASE: Point = new Point(CURVE.Gx, CURVE.Gy);
    /**
-     * Identity point aka point at infinity. point = point + zero_point
+     * Identity point aka point at infinity. p - p = zero_p; p + zero_p = p
     */
    static ZERO: Point = new Point(Fp.ZERO, Fp.ZERO);

@@ -583,7 +577,7 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
     * @param hex short/long ECDSA hex
     */
    static fromHex(hex: Hex): Point {
-      const { x, y } = CURVE.fromBytes(ensureBytes(hex));
+      const { x, y } = CURVE.fromBytes(ut.ensureBytes(hex));
      const point = new Point(x, y);
      point.assertValidity();
      return point;
@@ -607,16 +601,16 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
      // Zero is valid point too!
      if (this.equals(Point.ZERO)) {
        if (CURVE.allowInfinityPoint) return;
-        throw new Error('Point is infinity');
+        throw new Error('Point at infinity');
      }
      // Some 3rd-party test vectors require different wording between here & `fromCompressedHex`
      const msg = 'Point is not on elliptic curve';
      const { x, y } = this;
+      // Check if x, y are valid field elements
      if (!Fp.isValid(x) || !Fp.isValid(y)) throw new Error(msg);
-      const left = Fp.square(y);
-      const right = weierstrassEquation(x);
+      const left = Fp.square(y); // y²
+      const right = weierstrassEquation(x); // x³ + ax + b
      if (!Fp.equals(left, right)) throw new Error(msg);
-      // TODO: flag to disable this?
      if (!this.isTorsionFree()) throw new Error('Point must be of prime-order subgroup');
    }

@@ -630,34 +624,37 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
      return new Point(this.x, Fp.negate(this.y));
    }

+    protected toProj() {
+      return ProjectivePoint.fromAffine(this);
+    }
+
    // Adds point to itself
    double() {
-      return ProjectivePoint.fromAffine(this).double().toAffine();
+      return this.toProj().double().toAffine();
    }

-    // Adds point to other point
    add(other: Point) {
-      return ProjectivePoint.fromAffine(this).add(ProjectivePoint.fromAffine(other)).toAffine();
+      return this.toProj().add(ProjectivePoint.fromAffine(other)).toAffine();
    }

-    // Subtracts other point from the point
    subtract(other: Point) {
      return this.add(other.negate());
    }

    multiply(scalar: number | bigint) {
-      return ProjectivePoint.fromAffine(this).multiply(scalar, this).toAffine();
+      return this.toProj().multiply(scalar, this).toAffine();
    }

    multiplyUnsafe(scalar: bigint) {
-      return ProjectivePoint.fromAffine(this).multiplyUnsafe(scalar).toAffine();
+      return this.toProj().multiplyUnsafe(scalar).toAffine();
    }
+
    clearCofactor() {
-      return ProjectivePoint.fromAffine(this).clearCofactor().toAffine();
+      return this.toProj().clearCofactor().toAffine();
    }

    isTorsionFree(): boolean {
-      return ProjectivePoint.fromAffine(this).isTorsionFree();
+      return this.toProj().isTorsionFree();
    }

    /**
@@ -667,7 +664,7 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
     * @returns non-zero affine point
     */
    multiplyAndAddUnsafe(Q: Point, a: bigint, b: bigint): Point | undefined {
-      const P = ProjectivePoint.fromAffine(this);
+      const P = this.toProj();
      const aP =
        a === _0n || a === _1n || this !== Point.BASE ? P.multiplyUnsafe(a) : P.multiply(a);
      const bQ = ProjectivePoint.fromAffine(Q).multiplyUnsafe(b);
@@ -678,23 +675,26 @@ export function weierstrassPoints<T>(opts: CurvePointsType<T>) {
    // Encodes byte string to elliptic curve
    // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-3
    static hashToCurve(msg: Hex, options?: Partial<htfOpts>) {
-      if (!CURVE.mapToCurve) throw new Error('No mapToCurve defined for curve');
-      msg = ensureBytes(msg);
+      const { mapToCurve } = CURVE;
+      if (!mapToCurve) throw new Error('CURVE.mapToCurve() has not been defined');
+      msg = ut.ensureBytes(msg);
      const u = hash_to_field(msg, 2, { ...CURVE.htfDefaults, ...options } as htfOpts);
-      const { x: x0, y: y0 } = CURVE.mapToCurve(u[0]);
-      const { x: x1, y: y1 } = CURVE.mapToCurve(u[1]);
-      const p = new Point(x0, y0).add(new Point(x1, y1)).clearCofactor();
-      return p;
+      const { x: x0, y: y0 } = mapToCurve(u[0]);
+      const { x: x1, y: y1 } = mapToCurve(u[1]);
+      return new Point(x0, y0).add(new Point(x1, y1)).clearCofactor();
    }
+
    // https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-16#section-3
    static encodeToCurve(msg: Hex, options?: Partial<htfOpts>) {
-      if (!CURVE.mapToCurve) throw new Error('No mapToCurve defined for curve');
-      msg = ensureBytes(msg);
+      const { mapToCurve } = CURVE;
+      if (!mapToCurve) throw new Error('CURVE.mapToCurve() has not been defined');
+      msg = ut.ensureBytes(msg);
      const u = hash_to_field(msg, 1, { ...CURVE.htfDefaults, ...options } as htfOpts);
-      const { x, y } = CURVE.mapToCurve(u[0]);
+      const { x, y } = mapToCurve(u[0]);
      return new Point(x, y).clearCofactor();
    }
  }
+
  return {
    Point: Point as PointConstructor<T>,
    ProjectivePoint: ProjectivePoint as ProjectiveConstructor<T>,
@@ -733,15 +733,15 @@ export type CurveType = BasicCurve<bigint> & {
  // Default options
  lowS?: boolean;
  // Hashes
-  hash: utils.CHash; // Because we need outputLen for DRBG
+  hash: ut.CHash; // Because we need outputLen for DRBG
  hmac: HmacFnSync;
  randomBytes: (bytesLength?: number) => Uint8Array;
  truncateHash?: (hash: Uint8Array, truncateOnly?: boolean) => bigint;
 };

 function validateOpts(curve: CurveType) {
-  const opts = utils.validateOpts(curve);
-  if (typeof opts.hash !== 'function' || !Number.isSafeInteger(opts.hash.outputLen))
+  const opts = ut.validateOpts(curve);
+  if (typeof opts.hash !== 'function' || !ut.isPositiveInt(opts.hash.outputLen))
    throw new Error('Invalid hash function');
  if (typeof opts.hmac !== 'function') throw new Error('Invalid hmac function');
  if (typeof opts.randomBytes !== 'function') throw new Error('Invalid randomBytes function');
@@ -754,6 +754,7 @@ export type CurveFn = {
  getPublicKey: (privateKey: PrivKey, isCompressed?: boolean) => Uint8Array;
  getSharedSecret: (privateA: PrivKey, publicB: PubKey, isCompressed?: boolean) => Uint8Array;
  sign: (msgHash: Hex, privKey: PrivKey, opts?: SignOpts) => SignatureType;
+  signUnhashed: (msg: Uint8Array, privKey: PrivKey, opts?: SignOpts) => SignatureType;
  verify: (
    signature: Hex | SignatureType,
    msgHash: Hex,
@@ -766,8 +767,6 @@ export type CurveFn = {
  ProjectivePoint: ProjectiveConstructor<bigint>;
  Signature: SignatureConstructor;
  utils: {
-    mod: (a: bigint, b?: bigint) => bigint;
-    invert: (number: bigint, modulo?: bigint) => bigint;
    _bigintToBytes: (num: bigint) => Uint8Array;
    _bigintToString: (num: bigint) => string;
    _normalizePrivateKey: (key: PrivKey) => bigint;
@@ -824,19 +823,17 @@ class HmacDrbg {
      out.push(sl);
      len += this.v.length;
    }
-    return concatBytes(...out);
+    return ut.concatBytes(...out);
  }
-  // There is no need in clean() method
-  // It's useless, there are no guarantees with JS GC
-  // whether bigints are removed even if you clean Uint8Arrays.
+  // There are no guarantees with JS GC whether bigints are removed even if you clean Uint8Arrays.
 }

 export function weierstrass(curveDef: CurveType): CurveFn {
  const CURVE = validateOpts(curveDef) as ReturnType<typeof validateOpts>;
  const CURVE_ORDER = CURVE.n;
  const Fp = CURVE.Fp;
-  const compressedLen = Fp.BYTES + 1; // 33
-  const uncompressedLen = 2 * Fp.BYTES + 1; // 65
+  const compressedLen = Fp.BYTES + 1; // e.g. 33 for 32
+  const uncompressedLen = 2 * Fp.BYTES + 1; // e.g. 65 for 32

  function isValidFieldElement(num: bigint): boolean {
    // 0 is disallowed by arbitrary reasons. Probably because infinity point?
@@ -847,10 +844,12 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    weierstrassPoints({
      ...CURVE,
      toBytes(c, point, isCompressed: boolean): Uint8Array {
+        const x = Fp.toBytes(point.x);
+        const cat = ut.concatBytes;
        if (isCompressed) {
-          return concatBytes(new Uint8Array([point.hasEvenY() ? 0x02 : 0x03]), Fp.toBytes(point.x));
+          return cat(Uint8Array.from([point.hasEvenY() ? 0x02 : 0x03]), x);
        } else {
-          return concatBytes(new Uint8Array([0x04]), Fp.toBytes(point.x), Fp.toBytes(point.y));
+          return cat(Uint8Array.from([0x04]), x, Fp.toBytes(point.y));
        }
      },
      fromBytes(bytes: Uint8Array) {
@@ -858,7 +857,7 @@ export function weierstrass(curveDef: CurveType): CurveFn {
        const header = bytes[0];
        // this.assertValidity() is done inside of fromHex
        if (len === compressedLen && (header === 0x02 || header === 0x03)) {
-          const x = bytesToNumberBE(bytes.subarray(1));
+          const x = ut.bytesToNumberBE(bytes.subarray(1));
          if (!isValidFieldElement(x)) throw new Error('Point is not on curve');
          const y2 = weierstrassEquation(x); // y² = x³ + ax + b
          let y = Fp.sqrt(y2); // y = y² ^ (p+1)/4
@@ -910,15 +909,18 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    return isBiggerThanHalfOrder(s) ? mod.mod(-s, CURVE_ORDER) : s;
  }

+  function bits2int_2(bytes: Uint8Array): bigint {
+    const delta = bytes.length * 8 - CURVE.nBitLength;
+    const num = ut.bytesToNumberBE(bytes);
+    return delta > 0 ? num >> BigInt(delta) : num;
+  }
+
  // Ensures ECDSA message hashes are 32 bytes and < curve order
  function _truncateHash(hash: Uint8Array, truncateOnly = false): bigint {
-    const { n, nBitLength } = CURVE;
-    const byteLength = hash.length;
-    const delta = byteLength * 8 - nBitLength; // size of curve.n (252 bits)
-    let h = bytesToNumberBE(hash);
-    if (delta > 0) h = h >> BigInt(delta);
-    if (!truncateOnly && h >= n) h -= n;
-    return h;
+    const h = bits2int_2(hash);
+    if (truncateOnly) return h;
+    const { n } = CURVE;
+    return h >= n ? h - n : h;
  }
  const truncateHash = CURVE.truncateHash || _truncateHash;

@@ -930,15 +932,17 @@ export function weierstrass(curveDef: CurveType): CurveFn {
      this.assertValidity();
    }

-    // pair (32 bytes of r, 32 bytes of s)
+    // pair (bytes of r, bytes of s)
    static fromCompact(hex: Hex) {
      const arr = hex instanceof Uint8Array;
      const name = 'Signature.fromCompact';
      if (typeof hex !== 'string' && !arr)
        throw new TypeError(`${name}: Expected string or Uint8Array`);
      const str = arr ? bytesToHex(hex) : hex;
-      if (str.length !== 128) throw new Error(`${name}: Expected 64-byte hex`);
-      return new Signature(hexToNumber(str.slice(0, 64)), hexToNumber(str.slice(64, 128)));
+      const gl = CURVE.nByteLength * 2; // group length in hex, not ui8a
+      if (str.length !== 2 * gl) throw new Error(`${name}: Expected ${gl / 2}-byte hex`);
+      const slice = (from: number, to: number) => ut.hexToNumber(str.slice(from, to));
+      return new Signature(slice(0, gl), slice(gl, 2 * gl));
    }

    // DER encoded ECDSA signature
@@ -947,7 +951,7 @@ export function weierstrass(curveDef: CurveType): CurveFn {
      const arr = hex instanceof Uint8Array;
      if (typeof hex !== 'string' && !arr)
        throw new TypeError(`Signature.fromDER: Expected string or Uint8Array`);
-      const { r, s } = parseDERSignature(arr ? hex : hexToBytes(hex));
+      const { r, s } = DER.parseSig(arr ? hex : ut.hexToBytes(hex));
      return new Signature(r, s);
    }

@@ -964,6 +968,7 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    /**
     * Recovers public key from signature with recovery bit. Throws on invalid hash.
     * https://en.wikipedia.org/wiki/Elliptic_Curve_Digital_Signature_Algorithm#Public_key_recovery
+     * It's also possible to recover key without bit: try all 4 bit values and check for sig match.
     *
     * ```
     * recover(r, s, h) where
@@ -978,16 +983,18 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    recoverPublicKey(msgHash: Hex): Point {
      const { r, s, recovery } = this;
      if (recovery == null) throw new Error('Cannot recover: recovery bit is not present');
-      if (recovery !== 0 && recovery !== 1) throw new Error('Cannot recover: invalid recovery bit');
-      const h = truncateHash(ensureBytes(msgHash));
+      if (![0, 1, 2, 3].includes(recovery)) throw new Error('Cannot recover: invalid recovery bit');
+      const h = truncateHash(ut.ensureBytes(msgHash));
      const { n } = CURVE;
-      const rinv = mod.invert(r, n);
+      const radj = recovery === 2 || recovery === 3 ? r + n : r;
+      if (radj >= Fp.ORDER) throw new Error('Cannot recover: bit 2/3 is invalid with current r');
+      const rinv = mod.invert(radj, n);
      // Q = u1⋅G + u2⋅R
      const u1 = mod.mod(-h * rinv, n);
      const u2 = mod.mod(s * rinv, n);
      const prefix = recovery & 1 ? '03' : '02';
-      const R = Point.fromHex(prefix + numToFieldStr(r));
-      const Q = Point.BASE.multiplyAndAddUnsafe(R, u1, u2);
+      const R = Point.fromHex(prefix + numToFieldStr(radj));
+      const Q = Point.BASE.multiplyAndAddUnsafe(R, u1, u2); // unsafe is fine: no priv data leaked
      if (!Q) throw new Error('Cannot recover: point at infinify');
      Q.assertValidity();
      return Q;
@@ -1009,22 +1016,24 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    }

    // DER-encoded
-    toDERRawBytes(isCompressed = false) {
-      return hexToBytes(this.toDERHex(isCompressed));
+    toDERRawBytes() {
+      return ut.hexToBytes(this.toDERHex());
    }
-    toDERHex(isCompressed = false) {
-      const sHex = sliceDER(numberToHexUnpadded(this.s));
-      if (isCompressed) return sHex;
-      const rHex = sliceDER(numberToHexUnpadded(this.r));
-      const rLen = numberToHexUnpadded(rHex.length / 2);
-      const sLen = numberToHexUnpadded(sHex.length / 2);
-      const length = numberToHexUnpadded(rHex.length / 2 + sHex.length / 2 + 4);
+    toDERHex() {
+      const { numberToHexUnpadded: toHex } = ut;
+      const sHex = DER.slice(toHex(this.s));
+      const rHex = DER.slice(toHex(this.r));
+      const sHexL = sHex.length / 2;
+      const rHexL = rHex.length / 2;
+      const sLen = toHex(sHexL);
+      const rLen = toHex(rHexL);
+      const length = toHex(rHexL + sHexL + 4);
      return `30${length}02${rLen}${rHex}02${sLen}${sHex}`;
    }

-    // 32 bytes of r, then 32 bytes of s
+    // padded bytes of r, then padded bytes of s
    toCompactRawBytes() {
-      return hexToBytes(this.toCompactHex());
+      return ut.hexToBytes(this.toCompactHex());
    }
    toCompactHex() {
      return numToFieldStr(this.r) + numToFieldStr(this.s);
@@ -1032,8 +1041,6 @@ export function weierstrass(curveDef: CurveType): CurveFn {
  }

  const utils = {
-    mod: (n: bigint, modulo = Fp.ORDER) => mod.mod(n, modulo),
-    invert: Fp.invert,
    isValidPrivateKey(privateKey: PrivKey) {
      try {
        normalizePrivateKey(privateKey);
@@ -1053,7 +1060,8 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    /**
     * Converts some bytes to a valid private key. Needs at least (nBitLength+64) bytes.
     */
-    hashToPrivateKey: (hash: Hex): Uint8Array => numToField(hashToPrivateScalar(hash, CURVE_ORDER)),
+    hashToPrivateKey: (hash: Hex): Uint8Array =>
+      numToField(ut.hashToPrivateScalar(hash, CURVE_ORDER)),

    /**
     * Produces cryptographically secure private key from random of size (nBitLength+64)
@@ -1078,10 +1086,10 @@ export function weierstrass(curveDef: CurveType): CurveFn {
  };

  /**
-   * Computes public key for a private key.
+   * Computes public key for a private key. Checks for validity of the private key.
   * @param privateKey private key
-   * @param isCompressed whether to return compact, or full key
-   * @returns Public key, full by default; short when isCompressed=true
+   * @param isCompressed whether to return compact (default), or full key
+   * @returns Public key, full when isCompressed=false; short when isCompressed=true
   */
  function getPublicKey(privateKey: PrivKey, isCompressed = false): Uint8Array {
    return Point.fromPrivateKey(privateKey).toRawBytes(isCompressed);
@@ -1101,12 +1109,12 @@ export function weierstrass(curveDef: CurveType): CurveFn {
  }

  /**
-   * ECDH (Elliptic Curve Diffie Hellman) implementation.
-   * 1. Checks for validity of private key
-   * 2. Checks for the public key of being on-curve
+   * ECDH (Elliptic Curve Diffie Hellman).
+   * Computes shared public key from private key and public key.
+   * Checks: 1) private key validity 2) shared key is on-curve
   * @param privateA private key
   * @param publicB different public key
-   * @param isCompressed whether to return compact (33-byte), or full (65-byte) key
+   * @param isCompressed whether to return compact (default), or full key
   * @returns shared public key
   */
  function getSharedSecret(privateA: PrivKey, publicB: PubKey, isCompressed = false): Uint8Array {
@@ -1118,9 +1126,21 @@ export function weierstrass(curveDef: CurveType): CurveFn {
  }

  // RFC6979 methods
-  function bits2int(bytes: Uint8Array) {
-    const slice = bytes.length > Fp.BYTES ? bytes.slice(0, Fp.BYTES) : bytes;
-    return bytesToNumberBE(slice);
+  function bits2int(bytes: Uint8Array): bigint {
+    const { nByteLength } = CURVE;
+    if (!(bytes instanceof Uint8Array)) throw new Error('Expected Uint8Array');
+    const slice = bytes.length > nByteLength ? bytes.slice(0, nByteLength) : bytes;
+    // const slice = bytes; nByteLength; nBitLength;
+    let num = ut.bytesToNumberBE(slice);
+    // const { nBitLength } = CURVE;
+    // const delta = (bytes.length * 8) - nBitLength;
+    // if (delta > 0) {
+    //   // console.log('bits=', bytes.length*8, 'CURVE n=', nBitLength, 'delta=', delta);
+    //   // console.log(bytes.length, nBitLength, delta);
+    //   // console.log(bytes, new Error().stack);
+    //   num >>= BigInt(delta);
+    // }
+    return num;
  }
  function bits2octets(bytes: Uint8Array): Uint8Array {
    const z1 = bits2int(bytes);
@@ -1128,28 +1148,28 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    return int2octets(z2 < _0n ? z1 : z2);
  }
  function int2octets(num: bigint): Uint8Array {
-    return numToField(num); // prohibits >32 bytes
+    return numToField(num); // prohibits >nByteLength bytes
  }
  // Steps A, D of RFC6979 3.2
  // Creates RFC6979 seed; converts msg/privKey to numbers.
  function initSigArgs(msgHash: Hex, privateKey: PrivKey, extraEntropy?: Entropy) {
    if (msgHash == null) throw new Error(`sign: expected valid message hash, not "${msgHash}"`);
    // Step A is ignored, since we already provide hash instead of msg
-    const h1 = numToField(truncateHash(ensureBytes(msgHash)));
+    const h1 = numToField(truncateHash(ut.ensureBytes(msgHash)));
    const d = normalizePrivateKey(privateKey);
    // K = HMAC_K(V || 0x00 || int2octets(x) || bits2octets(h1) || k')
    const seedArgs = [int2octets(d), bits2octets(h1)];
    // RFC6979 3.6: additional k' could be provided
    if (extraEntropy != null) {
      if (extraEntropy === true) extraEntropy = CURVE.randomBytes(Fp.BYTES);
-      const e = ensureBytes(extraEntropy);
+      const e = ut.ensureBytes(extraEntropy);
      if (e.length !== Fp.BYTES) throw new Error(`sign: Expected ${Fp.BYTES} bytes of extra data`);
      seedArgs.push(e);
    }
    // seed is constructed from private key and message
    // Step D
    // V, 0x00 are done in HmacDRBG constructor.
-    const seed = concatBytes(...seedArgs);
+    const seed = ut.concatBytes(...seedArgs);
    const m = bits2int(h1);
    return { seed, m, d };
  }
@@ -1172,9 +1192,10 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    // r = x mod n
    const r = mod.mod(q.x, n);
    if (r === _0n) return;
-    // s = (1/k * (m + dr) mod n
+    // s = (m + dr)/k mod n where x/k == x*inv(k)
    const s = mod.mod(kinv * mod.mod(m + mod.mod(d * r, n), n), n);
    if (s === _0n) return;
+    // recovery bit is usually 0 or 1; rarely it's 2 or 3, when q.x > n
    let recovery = (q.x === r ? 0 : 2) | Number(q.y & _1n);
    let normS = s;
    if (lowS && isBiggerThanHalfOrder(s)) {
@@ -1184,11 +1205,19 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    return new Signature(r, normS, recovery);
  }

+  const defaultSigOpts: SignOpts = { lowS: CURVE.lowS };
+
  /**
   * Signs message hash (not message: you need to hash it by yourself).
+   * ```
+   * sign(m, d, k) where
+   *   (x, y) = G × k
+   *   r = x mod n
+   *   s = (m + dr)/k mod n
+   * ```
   * @param opts `lowS, extraEntropy`
   */
-  function sign(msgHash: Hex, privKey: PrivKey, opts: SignOpts = { lowS: CURVE.lowS }): Signature {
+  function sign(msgHash: Hex, privKey: PrivKey, opts = defaultSigOpts): Signature {
    // Steps A, D of RFC6979 3.2.
    const { seed, m, d } = initSigArgs(msgHash, privKey, opts.extraEntropy);
    // Steps B, C, D, E, F, G
@@ -1199,6 +1228,14 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    while (!(sig = kmdToSig(drbg.generateSync(), m, d, opts.lowS))) drbg.reseedSync();
    return sig;
  }
+
+  /**
+   * Signs a message (not message hash).
+   */
+  function signUnhashed(msg: Uint8Array, privKey: PrivKey, opts = defaultSigOpts): Signature {
+    return sign(CURVE.hash(ut.ensureBytes(msg)), privKey, opts);
+  }
+
  // Enable precomputes. Slows down first publicKey computation by 20ms.
  Point.BASE._setWindowSize(8);

@@ -1234,7 +1271,7 @@ export function weierstrass(curveDef: CurveType): CurveFn {
          signature = Signature.fromCompact(signature as Hex);
        }
      }
-      msgHash = ensureBytes(msgHash);
+      msgHash = ut.ensureBytes(msgHash);
    } catch (error) {
      return false;
    }
@@ -1265,6 +1302,7 @@ export function weierstrass(curveDef: CurveType): CurveFn {
    getPublicKey,
    getSharedSecret,
    sign,
+    signUnhashed,
    verify,
    Point,
    ProjectivePoint,
--- a/src/bls12-381.ts
+++ b/src/bls12-381.ts
@@ -1,4 +1,16 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
+
+// The pairing-friendly Barreto-Lynn-Scott elliptic curve construction allows to:
+// - Construct zk-SNARKs at the 128-bit security
+// - Use threshold signatures, which allows a user to sign lots of messages with one signature and verify them swiftly in a batch, using Boneh-Lynn-Shacham signature scheme.
+// Differences from @noble/bls12-381 1.4:
+// - PointG1 -> G1.Point
+// - PointG2 -> G2.Point
+// - PointG2.fromSignature -> Signature.decode
+// - PointG2.toSignature ->  Signature.encode
+// - Fixed Fp2 ORDER
+// - Points now have only two coordinates
+
 import { sha256 } from '@noble/hashes/sha256';
 import { randomBytes } from '@noble/hashes/utils';
 import { bls, CurveFn } from './abstract/bls.js';
@@ -23,14 +35,6 @@ import {
 } from './abstract/weierstrass.js';
 import { isogenyMap } from './abstract/hash-to-curve.js';

-// Differences from bls12-381:
-// - PointG1 -> G1.Point
-// - PointG2 -> G2.Point
-// - PointG2.fromSignature -> Signature.decode
-// - PointG2.toSignature ->  Signature.encode
-// - Fixed Fp2 ORDER
-// Points now have only two coordinates
-
 // CURVE FIELDS
 // Finite field over p.
 const Fp =
@@ -131,6 +135,7 @@ const Fp2: mod.Field<Fp2> & Fp2Utils = {
    return { c0: Fp.mul(factor, Fp.create(a)), c1: Fp.mul(factor, Fp.create(-b)) };
  },
  sqrt: (num) => {
+    if (Fp2.equals(num, Fp2.ZERO)) return Fp2.ZERO; // Algo doesn't handles this case
    // TODO: Optimize this line. It's extremely slow.
    // Speeding this up would boost aggregateSignatures.
    // https://eprint.iacr.org/2012/685.pdf applicable?
@@ -926,12 +931,12 @@ const htfDefaults = {
  k: 128,
  // option to use a message that has already been processed by
  // expand_message_xmd
-  expand: true,
+  expand: 'xmd',
  // Hash functions for: expand_message_xmd is appropriate for use with a
  // wide range of hash functions, including SHA-2, SHA-3, BLAKE2, and others.
  // BBS+ uses blake2: https://github.com/hyperledger/aries-framework-go/issues/2247
  hash: sha256,
-};
+} as const;

 // Encoding utils
 // Point on G1 curve: (x, y)
--- a/src/ed25519.ts
+++ b/src/ed25519.ts
@@ -3,7 +3,7 @@ import { sha512 } from '@noble/hashes/sha512';
 import { concatBytes, randomBytes, utf8ToBytes } from '@noble/hashes/utils';
 import { twistedEdwards, ExtendedPointType } from './abstract/edwards.js';
 import { montgomery } from './abstract/montgomery.js';
-import { mod, pow2, isNegativeLE, Fp as Field } from './abstract/modular.js';
+import { mod, pow2, isNegativeLE, Fp as Field, FpSqrtEven } from './abstract/modular.js';
 import {
  ensureBytes,
  equalBytes,
@@ -91,7 +91,80 @@ export const ED25519_TORSION_SUBGROUP = [
  'c7176a703d4dd84fba3c0b760d10670f2a2053fa2c39ccc64ec7fd7792ac03fa',
 ];

-const Fp = Field(ED25519_P);
+const Fp = Field(ED25519_P, undefined, true);
+
+// Hash To Curve Elligator2 Map (NOTE: different from ristretto255 elligator)
+// NOTE: very important part is usage of FpSqrtEven for ELL2_C1_EDWARDS, since
+// SageMath returns different root first and everything falls apart
+
+const ELL2_C1 = (Fp.ORDER + BigInt(3)) / BigInt(8); // 1. c1 = (q + 3) / 8       # Integer arithmetic
+
+const ELL2_C2 = Fp.pow(_2n, ELL2_C1); // 2. c2 = 2^c1
+const ELL2_C3 = Fp.sqrt(Fp.negate(Fp.ONE)); // 3. c3 = sqrt(-1)
+const ELL2_C4 = (Fp.ORDER - BigInt(5)) / BigInt(8); // 4. c4 = (q - 5) / 8       # Integer arithmetic
+const ELL2_J = BigInt(486662);
+
+// prettier-ignore
+function map_to_curve_elligator2_curve25519(u: bigint) {
+  let tv1 = Fp.square(u);       //  1.  tv1 = u^2
+  tv1 = Fp.mul(tv1, _2n);       //  2.  tv1 = 2 * tv1
+  let xd = Fp.add(tv1, Fp.ONE); //  3.   xd = tv1 + 1         # Nonzero: -1 is square (mod p), tv1 is not
+  let x1n = Fp.negate(ELL2_J);  //  4.  x1n = -J              # x1 = x1n / xd = -J / (1 + 2 * u^2)
+  let tv2 = Fp.square(xd);      //  5.  tv2 = xd^2
+  let gxd = Fp.mul(tv2, xd);    //  6.  gxd = tv2 * xd        # gxd = xd^3
+  let gx1 = Fp.mul(tv1, ELL2_J); //  7.  gx1 = J * tv1         # x1n + J * xd
+  gx1 = Fp.mul(gx1, x1n);       //  8.  gx1 = gx1 * x1n       # x1n^2 + J * x1n * xd
+  gx1 = Fp.add(gx1, tv2);       //  9.  gx1 = gx1 + tv2       # x1n^2 + J * x1n * xd + xd^2
+  gx1 = Fp.mul(gx1, x1n);       //  10. gx1 = gx1 * x1n       # x1n^3 + J * x1n^2 * xd + x1n * xd^2
+  let tv3 = Fp.square(gxd);     //  11. tv3 = gxd^2
+  tv2 = Fp.square(tv3);         //  12. tv2 = tv3^2           # gxd^4
+  tv3 = Fp.mul(tv3, gxd);       //  13. tv3 = tv3 * gxd       # gxd^3
+  tv3 = Fp.mul(tv3, gx1);       //  14. tv3 = tv3 * gx1       # gx1 * gxd^3
+  tv2 = Fp.mul(tv2, tv3);       //  15. tv2 = tv2 * tv3       # gx1 * gxd^7
+  let y11 = Fp.pow(tv2, ELL2_C4); //  16. y11 = tv2^c4        # (gx1 * gxd^7)^((p - 5) / 8)
+  y11 = Fp.mul(y11, tv3);       //  17. y11 = y11 * tv3       # gx1*gxd^3*(gx1*gxd^7)^((p-5)/8)
+  let y12 = Fp.mul(y11, ELL2_C3); //  18. y12 = y11 * c3
+  tv2 = Fp.square(y11);         //  19. tv2 = y11^2
+  tv2 = Fp.mul(tv2, gxd);       //  20. tv2 = tv2 * gxd
+  let e1 = Fp.equals(tv2, gx1); //  21.  e1 = tv2 == gx1
+  let y1 = Fp.cmov(y12, y11, e1); //  22.  y1 = CMOV(y12, y11, e1)  # If g(x1) is square, this is its sqrt
+  let x2n = Fp.mul(x1n, tv1);   //  23. x2n = x1n * tv1       # x2 = x2n / xd = 2 * u^2 * x1n / xd
+  let y21 = Fp.mul(y11, u);     //  24. y21 = y11 * u
+  y21 = Fp.mul(y21, ELL2_C2);   //  25. y21 = y21 * c2
+  let y22 = Fp.mul(y21, ELL2_C3); //  26. y22 = y21 * c3
+  let gx2 = Fp.mul(gx1, tv1);   //  27. gx2 = gx1 * tv1       # g(x2) = gx2 / gxd = 2 * u^2 * g(x1)
+  tv2 = Fp.square(y21);         //  28. tv2 = y21^2
+  tv2 = Fp.mul(tv2, gxd);       //  29. tv2 = tv2 * gxd
+  let e2 = Fp.equals(tv2, gx2); //  30.  e2 = tv2 == gx2
+  let y2 = Fp.cmov(y22, y21, e2); //  31.  y2 = CMOV(y22, y21, e2)  # If g(x2) is square, this is its sqrt
+  tv2 = Fp.square(y1);          //  32. tv2 = y1^2
+  tv2 = Fp.mul(tv2, gxd);       //  33. tv2 = tv2 * gxd
+  let e3 = Fp.equals(tv2, gx1); //  34.  e3 = tv2 == gx1
+  let xn = Fp.cmov(x2n, x1n, e3); //  35.  xn = CMOV(x2n, x1n, e3)  # If e3, x = x1, else x = x2
+  let y = Fp.cmov(y2, y1, e3);  //  36.   y = CMOV(y2, y1, e3)    # If e3, y = y1, else y = y2
+  let e4 = Fp.isOdd(y);         //  37.  e4 = sgn0(y) == 1        # Fix sign of y
+  y = Fp.cmov(y, Fp.negate(y), e3 !== e4); //  38.   y = CMOV(y, -y, e3 XOR e4)
+  return { xMn: xn, xMd: xd, yMn: y, yMd: 1n }; //  39. return (xn, xd, y, 1)
+}
+
+const ELL2_C1_EDWARDS = FpSqrtEven(Fp, Fp.negate(BigInt(486664))); // sgn0(c1) MUST equal 0
+function map_to_curve_elligator2_edwards25519(u: bigint) {
+  const { xMn, xMd, yMn, yMd } = map_to_curve_elligator2_curve25519(u); //  1.  (xMn, xMd, yMn, yMd) = map_to_curve_elligator2_curve25519(u)
+  let xn = Fp.mul(xMn, yMd); //  2.  xn = xMn * yMd
+  xn = Fp.mul(xn, ELL2_C1_EDWARDS); //  3.  xn = xn * c1
+  let xd = Fp.mul(xMd, yMn); //  4.  xd = xMd * yMn    # xn / xd = c1 * xM / yM
+  let yn = Fp.sub(xMn, xMd); //  5.  yn = xMn - xMd
+  let yd = Fp.add(xMn, xMd); //  6.  yd = xMn + xMd    # (n / d - 1) / (n / d + 1) = (n - d) / (n + d)
+  let tv1 = Fp.mul(xd, yd); //  7. tv1 = xd * yd
+  let e = Fp.equals(tv1, Fp.ZERO); //  8.   e = tv1 == 0
+  xn = Fp.cmov(xn, Fp.ZERO, e); //  9.  xn = CMOV(xn, 0, e)
+  xd = Fp.cmov(xd, Fp.ONE, e); //  10. xd = CMOV(xd, 1, e)
+  yn = Fp.cmov(yn, Fp.ONE, e); //  11. yn = CMOV(yn, 1, e)
+  yd = Fp.cmov(yd, Fp.ONE, e); //  12. yd = CMOV(yd, 1, e)
+
+  const inv = Fp.invertBatch([xd, yd]); // batch division
+  return { x: Fp.mul(xn, inv[0]), y: Fp.mul(yn, inv[1]) }; //  13. return (xn, xd, yn, yd)
+}

 const ED25519_DEF = {
  // Param: a
@@ -121,14 +194,10 @@ const ED25519_DEF = {
    p: Fp.ORDER,
    m: 1,
    k: 128,
-    expand: true,
+    expand: 'xmd',
    hash: sha512,
  },
-  mapToCurve: (scalars: bigint[]): { x: bigint; y: bigint } => {
-    throw new Error('Not supported yet');
-    // const { x, y } = calcElligatorRistrettoMap(scalars[0]).toAffine();
-    // return { x, y };
-  },
+  mapToCurve: (scalars: bigint[]) => map_to_curve_elligator2_edwards25519(scalars[0]),
 } as const;

 export const ed25519 = twistedEdwards(ED25519_DEF);
@@ -191,7 +260,7 @@ const invertSqrt = (number: bigint) => uvRatio(_1n, number);

 const MAX_255B = BigInt('0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff');
 const bytes255ToNumberLE = (bytes: Uint8Array) =>
-  ed25519.utils.mod(bytesToNumberLE(bytes) & MAX_255B);
+  ed25519.CURVE.Fp.create(bytesToNumberLE(bytes) & MAX_255B);

 type ExtendedPoint = ExtendedPointType;

@@ -200,7 +269,7 @@ type ExtendedPoint = ExtendedPointType;
 function calcElligatorRistrettoMap(r0: bigint): ExtendedPoint {
  const { d } = ed25519.CURVE;
  const P = ed25519.CURVE.Fp.ORDER;
-  const { mod } = ed25519.utils;
+  const mod = ed25519.CURVE.Fp.create;
  const r = mod(SQRT_M1 * r0 * r0); // 1
  const Ns = mod((r + _1n) * ONE_MINUS_D_SQ); // 2
  let c = BigInt(-1); // 3
@@ -258,7 +327,7 @@ export class RistrettoPoint {
    hex = ensureBytes(hex, 32);
    const { a, d } = ed25519.CURVE;
    const P = ed25519.CURVE.Fp.ORDER;
-    const { mod } = ed25519.utils;
+    const mod = ed25519.CURVE.Fp.create;
    const emsg = 'RistrettoPoint.fromHex: the hex is not valid encoding of RistrettoPoint';
    const s = bytes255ToNumberLE(hex);
    // 1. Check that s_bytes is the canonical encoding of a field element, or else abort.
@@ -288,7 +357,7 @@ export class RistrettoPoint {
  toRawBytes(): Uint8Array {
    let { x, y, z, t } = this.ep;
    const P = ed25519.CURVE.Fp.ORDER;
-    const { mod } = ed25519.utils;
+    const mod = ed25519.CURVE.Fp.create;
    const u1 = mod(mod(z + y) * mod(z - y)); // 1
    const u2 = mod(x * y); // 2
    // Square root always exists
@@ -326,7 +395,7 @@ export class RistrettoPoint {
    assertRstPoint(other);
    const a = this.ep;
    const b = other.ep;
-    const { mod } = ed25519.utils;
+    const mod = ed25519.CURVE.Fp.create;
    // (x1 * y2 == y1 * x2) | (y1 * y2 == x1 * x2)
    const one = mod(a.x * b.y) === mod(a.y * b.x);
    const two = mod(a.y * b.y) === mod(a.x * b.x);
--- a/src/ed448.ts
+++ b/src/ed448.ts
@@ -2,7 +2,7 @@
 import { shake256 } from '@noble/hashes/sha3';
 import { concatBytes, randomBytes, utf8ToBytes, wrapConstructor } from '@noble/hashes/utils';
 import { twistedEdwards } from './abstract/edwards.js';
-import { mod, pow2, Fp } from './abstract/modular.js';
+import { mod, pow2, Fp as Field } from './abstract/modular.js';
 import { montgomery } from './abstract/montgomery.js';

 /**
@@ -52,6 +52,83 @@ function adjustScalarBytes(bytes: Uint8Array): Uint8Array {
  return bytes;
 }

+const Fp = Field(ed448P, 456, true);
+
+// Hash To Curve Elligator2 Map
+const ELL2_C1 = (Fp.ORDER - BigInt(3)) / BigInt(4); // 1. c1 = (q - 3) / 4         # Integer arithmetic
+const ELL2_J = BigInt(156326);
+function map_to_curve_elligator2_curve448(u: bigint) {
+  let tv1 = Fp.square(u); // 1.  tv1 = u^2
+  let e1 = Fp.equals(tv1, Fp.ONE); // 2.   e1 = tv1 == 1
+  tv1 = Fp.cmov(tv1, Fp.ZERO, e1); // 3.  tv1 = CMOV(tv1, 0, e1)  # If Z * u^2 == -1, set tv1 = 0
+  let xd = Fp.sub(Fp.ONE, tv1); // 4.   xd = 1 - tv1
+  let x1n = Fp.negate(ELL2_J); // 5.  x1n = -J
+  let tv2 = Fp.square(xd); // 6.  tv2 = xd^2
+  let gxd = Fp.mul(tv2, xd); // 7.  gxd = tv2 * xd          # gxd = xd^3
+  let gx1 = Fp.mul(tv1, Fp.negate(ELL2_J)); // 8.  gx1 = -J * tv1          # x1n + J * xd
+  gx1 = Fp.mul(gx1, x1n); // 9.  gx1 = gx1 * x1n         # x1n^2 + J * x1n * xd
+  gx1 = Fp.add(gx1, tv2); // 10. gx1 = gx1 + tv2         # x1n^2 + J * x1n * xd + xd^2
+  gx1 = Fp.mul(gx1, x1n); // 11. gx1 = gx1 * x1n         # x1n^3 + J * x1n^2 * xd + x1n * xd^2
+  let tv3 = Fp.square(gxd); // 12. tv3 = gxd^2
+  tv2 = Fp.mul(gx1, gxd); // 13. tv2 = gx1 * gxd         # gx1 * gxd
+  tv3 = Fp.mul(tv3, tv2); // 14. tv3 = tv3 * tv2         # gx1 * gxd^3
+  let y1 = Fp.pow(tv3, ELL2_C1); // 15.  y1 = tv3^c1            # (gx1 * gxd^3)^((p - 3) / 4)
+  y1 = Fp.mul(y1, tv2); // 16.  y1 = y1 * tv2          # gx1 * gxd * (gx1 * gxd^3)^((p - 3) / 4)
+  let x2n = Fp.mul(x1n, Fp.negate(tv1)); // 17. x2n = -tv1 * x1n        # x2 = x2n / xd = -1 * u^2 * x1n / xd
+  let y2 = Fp.mul(y1, u); // 18.  y2 = y1 * u
+  y2 = Fp.cmov(y2, Fp.ZERO, e1); // 19.  y2 = CMOV(y2, 0, e1)
+  tv2 = Fp.square(y1); // 20. tv2 = y1^2
+  tv2 = Fp.mul(tv2, gxd); // 21. tv2 = tv2 * gxd
+  let e2 = Fp.equals(tv2, gx1); // 22.  e2 = tv2 == gx1
+  let xn = Fp.cmov(x2n, x1n, e2); // 23.  xn = CMOV(x2n, x1n, e2)  # If e2, x = x1, else x = x2
+  let y = Fp.cmov(y2, y1, e2); // 24.   y = CMOV(y2, y1, e2)    # If e2, y = y1, else y = y2
+  let e3 = Fp.isOdd(y); // 25.  e3 = sgn0(y) == 1        # Fix sign of y
+  y = Fp.cmov(y, Fp.negate(y), e2 !== e3); // 26.   y = CMOV(y, -y, e2 XOR e3)
+  return { xn, xd, yn: y, yd: Fp.ONE }; // 27. return (xn, xd, y, 1)
+}
+function map_to_curve_elligator2_edwards448(u: bigint) {
+  let { xn, xd, yn, yd } = map_to_curve_elligator2_curve448(u); // 1. (xn, xd, yn, yd) = map_to_curve_elligator2_curve448(u)
+  let xn2 = Fp.square(xn); // 2.  xn2 = xn^2
+  let xd2 = Fp.square(xd); // 3.  xd2 = xd^2
+  let xd4 = Fp.square(xd2); // 4.  xd4 = xd2^2
+  let yn2 = Fp.square(yn); // 5.  yn2 = yn^2
+  let yd2 = Fp.square(yd); // 6.  yd2 = yd^2
+  let xEn = Fp.sub(xn2, xd2); // 7.  xEn = xn2 - xd2
+  let tv2 = Fp.sub(xEn, xd2); // 8.  tv2 = xEn - xd2
+  xEn = Fp.mul(xEn, xd2); // 9.  xEn = xEn * xd2
+  xEn = Fp.mul(xEn, yd); // 10. xEn = xEn * yd
+  xEn = Fp.mul(xEn, yn); // 11. xEn = xEn * yn
+  xEn = Fp.mul(xEn, 4n); // 12. xEn = xEn * 4
+  tv2 = Fp.mul(tv2, xn2); // 13. tv2 = tv2 * xn2
+  tv2 = Fp.mul(tv2, yd2); // 14. tv2 = tv2 * yd2
+  let tv3 = Fp.mul(yn2, 4n); // 15. tv3 = 4 * yn2
+  let tv1 = Fp.add(tv3, yd2); // 16. tv1 = tv3 + yd2
+  tv1 = Fp.mul(tv1, xd4); // 17. tv1 = tv1 * xd4
+  let xEd = Fp.add(tv1, tv2); // 18. xEd = tv1 + tv2
+  tv2 = Fp.mul(tv2, xn); // 19. tv2 = tv2 * xn
+  let tv4 = Fp.mul(xn, xd4); // 20. tv4 = xn * xd4
+  let yEn = Fp.sub(tv3, yd2); // 21. yEn = tv3 - yd2
+  yEn = Fp.mul(yEn, tv4); // 22. yEn = yEn * tv4
+  yEn = Fp.sub(yEn, tv2); // 23. yEn = yEn - tv2
+  tv1 = Fp.add(xn2, xd2); // 24. tv1 = xn2 + xd2
+  tv1 = Fp.mul(tv1, xd2); // 25. tv1 = tv1 * xd2
+  tv1 = Fp.mul(tv1, xd); // 26. tv1 = tv1 * xd
+  tv1 = Fp.mul(tv1, yn2); // 27. tv1 = tv1 * yn2
+  tv1 = Fp.mul(tv1, BigInt(-2)); // 28. tv1 = -2 * tv1
+  let yEd = Fp.add(tv2, tv1); // 29. yEd = tv2 + tv1
+  tv4 = Fp.mul(tv4, yd2); // 30. tv4 = tv4 * yd2
+  yEd = Fp.add(yEd, tv4); // 31. yEd = yEd + tv4
+  tv1 = Fp.mul(xEd, yEd); // 32. tv1 = xEd * yEd
+  let e = Fp.equals(tv1, Fp.ZERO); // 33.   e = tv1 == 0
+  xEn = Fp.cmov(xEn, Fp.ZERO, e); // 34. xEn = CMOV(xEn, 0, e)
+  xEd = Fp.cmov(xEd, Fp.ONE, e); // 35. xEd = CMOV(xEd, 1, e)
+  yEn = Fp.cmov(yEn, Fp.ONE, e); // 36. yEn = CMOV(yEn, 1, e)
+  yEd = Fp.cmov(yEd, Fp.ONE, e); // 37. yEd = CMOV(yEd, 1, e)
+
+  const inv = Fp.invertBatch([xEd, yEd]); // batch division
+  return { x: Fp.mul(xEn, inv[0]), y: Fp.mul(yEn, inv[1]) }; // 38. return (xEn, xEd, yEn, yEd)
+}
+
 const ED448_DEF = {
  // Param: a
  a: BigInt(1),
@@ -60,8 +137,9 @@ const ED448_DEF = {
    '726838724295606890549323807888004534353641360687318060281490199180612328166730772686396383698676545930088884461843637361053498018326358'
  ),
  // Finite field 𝔽p over which we'll do calculations; 2n ** 448n - 2n ** 224n - 1n
-  Fp: Fp(ed448P, 456),
-  // Subgroup order: how many points ed448 has; 2n**446n - 13818066809895115352007386748515426880336692474882178609894547503885n
+  Fp,
+  // Subgroup order: how many points curve has;
+  // 2n**446n - 13818066809895115352007386748515426880336692474882178609894547503885n
  n: BigInt(
    '181709681073901722637330951972001133588410340171829515070372549795146003961539585716195755291692375963310293709091662304773755859649779'
  ),
@@ -111,6 +189,15 @@ const ED448_DEF = {
    // square root exists, and the decoding fails.
    return { isValid: mod(x2 * v, P) === u, value: x };
  },
+  htfDefaults: {
+    DST: 'edwards448_XOF:SHAKE256_ELL2_RO_',
+    p: Fp.ORDER,
+    m: 1,
+    k: 224,
+    expand: 'xof',
+    hash: shake256,
+  },
+  mapToCurve: (scalars: bigint[]) => map_to_curve_elligator2_edwards448(scalars[0]),
 } as const;

 export const ed448 = twistedEdwards(ED448_DEF);
--- a/src/jubjub.ts
+++ b/src/jubjub.ts
@@ -1,5 +1,5 @@
 /*! noble-curves - MIT License (c) 2022 Paul Miller (paulmillr.com) */
-import { sha256 } from '@noble/hashes/sha256';
+import { sha512 } from '@noble/hashes/sha512';
 import { concatBytes, randomBytes, utf8ToBytes } from '@noble/hashes/utils';
 import { twistedEdwards } from './abstract/edwards.js';
 import { blake2s } from '@noble/hashes/blake2s';
@@ -8,6 +8,7 @@ import { Fp } from './abstract/modular.js';
 /**
 * jubjub Twisted Edwards curve.
 * https://neuromancer.sk/std/other/JubJub
+ * jubjub does not use EdDSA, so `hash`/sha512 params are passed because interface expects them.
 */

 export const jubjub = twistedEdwards({
@@ -15,16 +16,16 @@ export const jubjub = twistedEdwards({
  a: BigInt('0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000000'),
  d: BigInt('0x2a9318e74bfa2b48f5fd9207e6bd7fd4292d7f6d37579d2601065fd6d6343eb1'),
  // Finite field 𝔽p over which we'll do calculations
+  // Same value as bls12-381 Fr (not Fp)
  Fp: Fp(BigInt('0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001')),
-  // Subgroup order: how many points ed25519 has
-  // 2n ** 252n + 27742317777372353535851937790883648493n;
+  // Subgroup order: how many points curve has
  n: BigInt('0xe7db4ea6533afa906673b0101343b00a6682093ccc81082d0970e5ed6f72cb7'),
  // Cofactor
  h: BigInt(8),
  // Base point (x, y) aka generator point
  Gx: BigInt('0x11dafe5d23e1218086a365b99fbf3d3be72f6afd7d1f72623e6b071492d1122b'),
  Gy: BigInt('0x1d523cf1ddab1a1793132e78c866c0c33e26ba5cc220fed7cc3f870e59d292aa'),
-  hash: sha256,
+  hash: sha512,
  randomBytes,
 } as const);

--- a/src/p256.ts
+++ b/src/p256.ts
@@ -37,7 +37,7 @@ export const P256 = createCurve(
      p: Fp.ORDER,
      m: 1,
      k: 128,
-      expand: true,
+      expand: 'xmd',
      hash: sha256,
    },
  } as const,
--- a/src/p384.ts
+++ b/src/p384.ts
@@ -41,7 +41,7 @@ export const P384 = createCurve({
      p: Fp.ORDER,
      m: 1,
      k: 192,
-      expand: true,
+      expand: 'xmd',
      hash: sha384,
    },
  } as const,
--- a/src/p521.ts
+++ b/src/p521.ts
@@ -54,7 +54,7 @@ export const P521 = createCurve({
    p: Fp.ORDER,
    m: 1,
    k: 256,
-    expand: true,
+      expand: 'xmd',
    hash: sha512,
  },
 } as const, sha512);
--- a/src/secp256k1.ts
+++ b/src/secp256k1.ts
@@ -15,10 +15,8 @@ import { randomBytes } from '@noble/hashes/utils';
 import { isogenyMap } from './abstract/hash-to-curve.js';

 /**
- * secp256k1 belongs to Koblitz curves: it has
- * efficiently computable Frobenius endomorphism.
- * Endomorphism improves efficiency:
- * Uses 2x less RAM, speeds up precomputation by 2x and ECDH / sign key recovery by 20%.
+ * secp256k1 belongs to Koblitz curves: it has efficiently computable endomorphism.
+ * Endomorphism uses 2x less RAM, speeds up precomputation by 2x and ECDH / key recovery by 20%.
 * Should always be used for Projective's double-and-add multiplication.
 * For affines cached multiplication, it trades off 1/2 init time & 1/3 ram for 20% perf hit.
 * https://gist.github.com/paulmillr/eb670806793e84df628a7c434a873066
@@ -56,7 +54,9 @@ function sqrtMod(y: bigint): bigint {
  const b223 = (pow2(b220, _3n, P) * b3) % P;
  const t1 = (pow2(b223, _23n, P) * b22) % P;
  const t2 = (pow2(t1, _6n, P) * b2) % P;
-  return pow2(t2, _2n, P);
+  const root = pow2(t2, _2n, P);
+  if (!Fp.equals(Fp.square(root), y)) throw new Error('Cannot find square root');
+  return root;
 }

 const Fp = Field(secp256k1P, undefined, undefined, { sqrt: sqrtMod });
@@ -127,7 +127,7 @@ export const secp256k1 = createCurve(
        const b1 = -_1n * BigInt('0xe4437ed6010e88286f547fa90abfe4c3');
        const a2 = BigInt('0x114ca50f7a8e2f3f657c1108d9d44cfd8');
        const b2 = a1;
-        const POW_2_128 = BigInt('0x100000000000000000000000000000000');
+        const POW_2_128 = BigInt('0x100000000000000000000000000000000'); // (2n**128n).toString(16)

        const c1 = divNearest(b2 * k, n);
        const c2 = divNearest(-b1 * k, n);
@@ -152,7 +152,7 @@ export const secp256k1 = createCurve(
      p: Fp.ORDER,
      m: 1,
      k: 128,
-      expand: true,
+      expand: 'xmd',
      hash: sha256,
    },
  },
@@ -173,7 +173,7 @@ function normalizePublicKey(publicKey: Hex | PointType<bigint>): PointType<bigin
  } else {
    const bytes = ensureBytes(publicKey);
    // Schnorr is 32 bytes
-    if (bytes.length === 32) {
+    if (bytes.length !== 32) throw new Error('Schnorr pubkeys must be 32 bytes');
    const x = bytesToNumberBE(bytes);
    if (!isValidFieldElement(x)) throw new Error('Point is not on curve');
    const y2 = secp256k1.utils._weierstrassEquation(x); // y² = x³ + ax + b
@@ -185,9 +185,6 @@ function normalizePublicKey(publicKey: Hex | PointType<bigint>): PointType<bigin
    point.assertValidity();
    return point;
  }
-    // Do we need that in schnorr at all?
-    return secp256k1.Point.fromHex(publicKey);
-  }
 }

 const isWithinCurveOrder = secp256k1.utils._isWithinCurveOrder;
@@ -225,10 +222,13 @@ class SchnorrSignature {
  }
  static fromHex(hex: Hex) {
    const bytes = ensureBytes(hex);
-    if (bytes.length !== 64)
-      throw new TypeError(`SchnorrSignature.fromHex: expected 64 bytes, not ${bytes.length}`);
-    const r = bytesToNumberBE(bytes.subarray(0, 32));
-    const s = bytesToNumberBE(bytes.subarray(32, 64));
+    const len = 32; // group length
+    if (bytes.length !== 2 * len)
+      throw new TypeError(
+        `SchnorrSignature.fromHex: expected ${2 * len} bytes, not ${bytes.length}`
+      );
+    const r = bytesToNumberBE(bytes.subarray(0, len));
+    const s = bytesToNumberBE(bytes.subarray(len, 2 * len));
    return new SchnorrSignature(r, s);
  }
  assertValidity() {
--- a/src/stark.ts
+++ b/src/stark.ts
@@ -97,7 +97,7 @@ function getSharedSecret0x(privKeyA: Hex, pubKeyB: Hex) {
  return starkCurve.getSharedSecret(normalizePrivateKey(privKeyA), pubKeyB);
 }

-function sign0x(msgHash: Hex, privKey: Hex, opts: any) {
+function sign0x(msgHash: Hex, privKey: Hex, opts?: any) {
  if (typeof privKey === 'string') privKey = strip0x(privKey).padStart(64, '0');
  return starkCurve.sign(ensureBytes0x(msgHash), normalizePrivateKey(privKey), opts);
 }
@@ -138,11 +138,12 @@ function hashKeyWithIndex(key: Uint8Array, index: number) {
 export function grindKey(seed: Hex) {
  const _seed = ensureBytes0x(seed);
  const sha256mask = 2n ** 256n;
-  const limit = sha256mask - starkCurve.utils.mod(sha256mask, starkCurve.CURVE.n);
+  const Fn = Fp(CURVE.n);
+  const limit = sha256mask - Fn.create(sha256mask);
  for (let i = 0; ; i++) {
    const key = hashKeyWithIndex(_seed, i);
    // key should be in [0, limit)
-    if (key < limit) return starkCurve.utils.mod(key, starkCurve.CURVE.n).toString(16);
+    if (key < limit) return Fn.create(key).toString(16);
  }
 }

@@ -203,6 +204,8 @@ function pedersenPrecompute(p1: ProjectivePoint, p2: ProjectivePoint): Projectiv
    out.push(p);
    p = p.double();
  }
+  // NOTE: we cannot use wNAF here, because last 4 bits will require full 248 bits multiplication
+  // We can add support for this to wNAF, but it will complicate wNAF.
  p = p2;
  for (let i = 0; i < 4; i++) {
    out.push(p);
--- a/test/basic.test.js
+++ b/test/basic.test.js
@@ -1,7 +1,8 @@
 import { deepStrictEqual, throws } from 'assert';
-import { should } from 'micro-should';
+import { should, describe } from 'micro-should';
 import * as fc from 'fast-check';
 import * as mod from '../lib/esm/abstract/modular.js';
+import { bytesToHex as toHex } from '../lib/esm/abstract/utils.js';
 // Generic tests for all curves in package
 import { secp192r1 } from '../lib/esm/p192.js';
 import { secp224r1 } from '../lib/esm/p224.js';
@@ -15,7 +16,286 @@ import { starkCurve } from '../lib/esm/stark.js';
 import { pallas, vesta } from '../lib/esm/pasta.js';
 import { bn254 } from '../lib/esm/bn.js';
 import { jubjub } from '../lib/esm/jubjub.js';
+import { bls12_381 } from '../lib/esm/bls12-381.js';

+// Fields tests
+const FIELDS = {
+  secp192r1: { Fp: [secp192r1.CURVE.Fp] },
+  secp224r1: { Fp: [secp224r1.CURVE.Fp] },
+  secp256r1: { Fp: [secp256r1.CURVE.Fp] },
+  secp521r1: { Fp: [secp521r1.CURVE.Fp] },
+  secp256k1: { Fp: [secp256k1.CURVE.Fp] },
+  stark: { Fp: [starkCurve.CURVE.Fp] },
+  jubjub: { Fp: [jubjub.CURVE.Fp] },
+  ed25519: { Fp: [ed25519.CURVE.Fp] },
+  ed448: { Fp: [ed448.CURVE.Fp] },
+  bn254: { Fp: [bn254.CURVE.Fp] },
+  pallas: { Fp: [pallas.CURVE.Fp] },
+  vesta: { Fp: [vesta.CURVE.Fp] },
+  bls12: {
+    Fp: [bls12_381.CURVE.Fp],
+    Fp2: [
+      bls12_381.CURVE.Fp2,
+      fc.array(fc.bigInt(1n, bls12_381.CURVE.Fp.ORDER - 1n), {
+        minLength: 2,
+        maxLength: 2,
+      }),
+      (Fp2, num) => Fp2.fromBigTuple([num[0], num[1]]),
+    ],
+    // Fp6: [bls12_381.CURVE.Fp6],
+    Fp12: [
+      bls12_381.CURVE.Fp12,
+      fc.array(fc.bigInt(1n, bls12_381.CURVE.Fp.ORDER - 1n), {
+        minLength: 12,
+        maxLength: 12,
+      }),
+      (Fp12, num) => Fp12.fromBigTwelve(num),
+    ],
+  },
+};
+
+for (const c in FIELDS) {
+  const curve = FIELDS[c];
+  for (const f in curve) {
+    const Fp = curve[f][0];
+    const name = `${c}/${f}:`;
+    const FC_BIGINT = curve[f][1] ? curve[f][1] : fc.bigInt(1n, Fp.ORDER - 1n);
+
+    const create = curve[f][2] ? curve[f][2].bind(null, Fp) : (num) => Fp.create(num);
+    describe(name, () => {
+      should('equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            const b = create(num);
+            deepStrictEqual(Fp.equals(a, b), true);
+            deepStrictEqual(Fp.equals(b, a), true);
+          })
+        );
+      });
+      should('non-equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, (num1, num2) => {
+            const a = create(num1);
+            const b = create(num2);
+            deepStrictEqual(Fp.equals(a, b), num1 === num2);
+            deepStrictEqual(Fp.equals(b, a), num1 === num2);
+          })
+        );
+      });
+      should('add/subtract/commutativity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, (num1, num2) => {
+            const a = create(num1);
+            const b = create(num2);
+            deepStrictEqual(Fp.add(a, b), Fp.add(b, a));
+          })
+        );
+      });
+      should('add/subtract/associativity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, FC_BIGINT, (num1, num2, num3) => {
+            const a = create(num1);
+            const b = create(num2);
+            const c = create(num3);
+            deepStrictEqual(Fp.add(a, Fp.add(b, c)), Fp.add(Fp.add(a, b), c));
+          })
+        );
+      });
+      should('add/subtract/x+0=x', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.add(a, Fp.ZERO), a);
+          })
+        );
+      });
+      should('add/subtract/x-0=x', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.sub(a, Fp.ZERO), a);
+            deepStrictEqual(Fp.sub(a, a), Fp.ZERO);
+          })
+        );
+      });
+      should('add/subtract/negate equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num1) => {
+            const a = create(num1);
+            const b = create(num1);
+            deepStrictEqual(Fp.sub(Fp.ZERO, a), Fp.negate(a));
+            deepStrictEqual(Fp.sub(a, b), Fp.add(a, Fp.negate(b)));
+            deepStrictEqual(Fp.sub(a, b), Fp.add(a, Fp.mul(b, Fp.create(-1n))));
+          })
+        );
+      });
+      should('add/subtract/negate', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.negate(a), Fp.sub(Fp.ZERO, a));
+            deepStrictEqual(Fp.negate(a), Fp.mul(a, Fp.create(-1n)));
+          })
+        );
+      });
+      should('negate(0)', () => {
+        deepStrictEqual(Fp.negate(Fp.ZERO), Fp.ZERO);
+      });
+
+      should('multiply/commutativity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, (num1, num2) => {
+            const a = create(num1);
+            const b = create(num2);
+            deepStrictEqual(Fp.mul(a, b), Fp.mul(b, a));
+          })
+        );
+      });
+      should('multiply/associativity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, FC_BIGINT, (num1, num2, num3) => {
+            const a = create(num1);
+            const b = create(num2);
+            const c = create(num3);
+            deepStrictEqual(Fp.mul(a, Fp.mul(b, c)), Fp.mul(Fp.mul(a, b), c));
+          })
+        );
+      });
+      should('multiply/distributivity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, FC_BIGINT, (num1, num2, num3) => {
+            const a = create(num1);
+            const b = create(num2);
+            const c = create(num3);
+            deepStrictEqual(Fp.mul(a, Fp.add(b, c)), Fp.add(Fp.mul(b, a), Fp.mul(c, a)));
+          })
+        );
+      });
+      should('multiply/add equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.mul(a, 0n), Fp.ZERO);
+            deepStrictEqual(Fp.mul(a, Fp.ZERO), Fp.ZERO);
+            deepStrictEqual(Fp.mul(a, 1n), a);
+            deepStrictEqual(Fp.mul(a, Fp.ONE), a);
+            deepStrictEqual(Fp.mul(a, 2n), Fp.add(a, a));
+            deepStrictEqual(Fp.mul(a, 3n), Fp.add(Fp.add(a, a), a));
+            deepStrictEqual(Fp.mul(a, 4n), Fp.add(Fp.add(Fp.add(a, a), a), a));
+          })
+        );
+      });
+      should('multiply/square equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.square(a), Fp.mul(a, a));
+          })
+        );
+      });
+      should('multiply/pow equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.pow(a, 0n), Fp.ONE);
+            deepStrictEqual(Fp.pow(a, 1n), a);
+            deepStrictEqual(Fp.pow(a, 2n), Fp.mul(a, a));
+            deepStrictEqual(Fp.pow(a, 3n), Fp.mul(Fp.mul(a, a), a));
+          })
+        );
+      });
+
+      should('square(0)', () => {
+        deepStrictEqual(Fp.square(Fp.ZERO), Fp.ZERO);
+        deepStrictEqual(Fp.mul(Fp.ZERO, Fp.ZERO), Fp.ZERO);
+      });
+
+      should('square(1)', () => {
+        deepStrictEqual(Fp.square(Fp.ONE), Fp.ONE);
+        deepStrictEqual(Fp.mul(Fp.ONE, Fp.ONE), Fp.ONE);
+      });
+
+      should('square(-1)', () => {
+        const minus1 = Fp.negate(Fp.ONE);
+        deepStrictEqual(Fp.square(minus1), Fp.ONE);
+        deepStrictEqual(Fp.mul(minus1, minus1), Fp.ONE);
+      });
+
+      const isSquare = mod.FpIsSquare(Fp);
+      // Not implemented
+      if (Fp !== bls12_381.CURVE.Fp12) {
+        should('multiply/sqrt', () => {
+          fc.assert(
+            fc.property(FC_BIGINT, (num) => {
+              const a = create(num);
+              let root;
+              try {
+                root = Fp.sqrt(a);
+              } catch (e) {
+                deepStrictEqual(isSquare(a), false);
+                return;
+              }
+              deepStrictEqual(isSquare(a), true);
+              deepStrictEqual(Fp.equals(Fp.square(root), a), true, 'sqrt(a)^2 == a');
+              deepStrictEqual(Fp.equals(Fp.square(Fp.negate(root)), a), true, '(-sqrt(a))^2 == a');
+            })
+          );
+        });
+
+        should('sqrt(0)', () => {
+          deepStrictEqual(Fp.sqrt(Fp.ZERO), Fp.ZERO);
+          const sqrt1 = Fp.sqrt(Fp.ONE);
+          deepStrictEqual(
+            Fp.equals(sqrt1, Fp.ONE) || Fp.equals(sqrt1, Fp.negate(Fp.ONE)),
+            true,
+            'sqrt(1) = 1 or -1'
+          );
+        });
+      }
+
+      should('div/division by one equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            if (Fp.equals(a, Fp.ZERO)) return; // No division by zero
+            deepStrictEqual(Fp.div(a, Fp.ONE), a);
+            deepStrictEqual(Fp.div(a, a), Fp.ONE);
+          })
+        );
+      });
+      should('zero division equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, (num) => {
+            const a = create(num);
+            deepStrictEqual(Fp.div(Fp.ZERO, a), Fp.ZERO);
+          })
+        );
+      });
+      should('div/division distributivity', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, FC_BIGINT, (num1, num2, num3) => {
+            const a = create(num1);
+            const b = create(num2);
+            const c = create(num3);
+            deepStrictEqual(Fp.div(Fp.add(a, b), c), Fp.add(Fp.div(a, c), Fp.div(b, c)));
+          })
+        );
+      });
+      should('div/division and multiplication equality', () => {
+        fc.assert(
+          fc.property(FC_BIGINT, FC_BIGINT, (num1, num2) => {
+            const a = create(num1);
+            const b = create(num2);
+            deepStrictEqual(Fp.div(a, b), Fp.mul(a, Fp.invert(b)));
+          })
+        );
+      });
+    });
+  }
+}
+
+// Group tests
 // prettier-ignore
 const CURVES = {
  secp192r1, secp224r1, secp256r1, secp384r1, secp521r1,
@@ -29,15 +309,16 @@ const CURVES = {
 };

 const NUM_RUNS = 5;
+
 const getXY = (p) => ({ x: p.x, y: p.y });

 function equal(a, b, comment) {
-  deepStrictEqual(a.equals(b), true, `eq(${comment})`);
+  deepStrictEqual(a.equals(b), true, 'eq(${comment})');
  if (a.toAffine && b.toAffine) {
-    deepStrictEqual(getXY(a.toAffine()), getXY(b.toAffine()), `eqToAffine(${comment})`);
+    deepStrictEqual(getXY(a.toAffine()), getXY(b.toAffine()), 'eqToAffine(${comment})');
  } else if (!a.toAffine && !b.toAffine) {
    // Already affine
-    deepStrictEqual(getXY(a), getXY(b), `eqAffine(${comment})`);
+    deepStrictEqual(getXY(a), getXY(b), 'eqAffine(${comment})');
  } else throw new Error('Different point types');
 }

@@ -62,46 +343,49 @@ for (const name in CURVES) {

    const G = [p.ZERO, p.BASE];
    for (let i = 2; i < 10; i++) G.push(G[1].multiply(i));
+    const title = `${name}/${pointName}`;
+    describe(title, () => {
+      describe('basic group laws', () => {
        // Here we check basic group laws, to verify that points works as group
-    should(`${name}/${pointName}/Basic group laws (zero)`, () => {
+        should('(zero)', () => {
          equal(G[0].double(), G[0], '(0*G).double() = 0');
          equal(G[0].add(G[0]), G[0], '0*G + 0*G = 0');
          equal(G[0].subtract(G[0]), G[0], '0*G - 0*G = 0');
          equal(G[0].negate(), G[0], '-0 = 0');
          for (let i = 0; i < G.length; i++) {
            const p = G[i];
-        equal(p, p.add(G[0]), `${i}*G + 0 = ${i}*G`);
-        equal(G[0].multiply(i + 1), G[0], `${i + 1}*0 = 0`);
+            equal(p, p.add(G[0]), '${i}*G + 0 = ${i}*G');
+            equal(G[0].multiply(i + 1), G[0], '${i + 1}*0 = 0');
          }
        });
-    should(`${name}/${pointName}/Basic group laws (one)`, () => {
+        should('(one)', () => {
          equal(G[1].double(), G[2], '(1*G).double() = 2*G');
          equal(G[1].subtract(G[1]), G[0], '1*G - 1*G = 0');
          equal(G[1].add(G[1]), G[2], '1*G + 1*G = 2*G');
        });
-    should(`${name}/${pointName}/Basic group laws (sanity tests)`, () => {
-      equal(G[2].double(), G[4], `(2*G).double() = 4*G`);
-      equal(G[2].add(G[2]), G[4], `2*G + 2*G = 4*G`);
-      equal(G[7].add(G[3].negate()), G[4], `7*G - 3*G = 4*G`);
+        should('(sanity tests)', () => {
+          equal(G[2].double(), G[4], '(2*G).double() = 4*G');
+          equal(G[2].add(G[2]), G[4], '2*G + 2*G = 4*G');
+          equal(G[7].add(G[3].negate()), G[4], '7*G - 3*G = 4*G');
        });
-    should(`${name}/${pointName}/Basic group laws (addition commutativity)`, () => {
-      equal(G[4].add(G[3]), G[3].add(G[4]), `4*G + 3*G = 3*G + 4*G`);
-      equal(G[4].add(G[3]), G[3].add(G[2]).add(G[2]), `4*G + 3*G = 3*G + 2*G + 2*G`);
+        should('(addition commutativity)', () => {
+          equal(G[4].add(G[3]), G[3].add(G[4]), '4*G + 3*G = 3*G + 4*G');
+          equal(G[4].add(G[3]), G[3].add(G[2]).add(G[2]), '4*G + 3*G = 3*G + 2*G + 2*G');
        });
-    should(`${name}/${pointName}/Basic group laws (double)`, () => {
+        should('(double)', () => {
          equal(G[3].double(), G[6], '(3*G).double() = 6*G');
        });
-    should(`${name}/${pointName}/Basic group laws (multiply)`, () => {
+        should('(multiply)', () => {
          equal(G[2].multiply(3), G[6], '(2*G).multiply(3) = 6*G');
        });
-    should(`${name}/${pointName}/Basic group laws (same point addition)`, () => {
-      equal(G[3].add(G[3]), G[6], `3*G + 3*G = 6*G`);
+        should('(same point addition)', () => {
+          equal(G[3].add(G[3]), G[6], '3*G + 3*G = 6*G');
        });
-    should(`${name}/${pointName}/Basic group laws (same point (negative) addition)`, () => {
+        should('(same point (negative) addition)', () => {
          equal(G[3].add(G[3].negate()), G[0], '3*G + (- 3*G) = 0*G');
          equal(G[3].subtract(G[3]), G[0], '3*G - 3*G = 0*G');
        });
-    should(`${name}/${pointName}/Basic group laws (curve order)`, () => {
+        should('(curve order)', () => {
          equal(G[1].multiply(CURVE_ORDER - 1n).add(G[1]), G[0], '(N-1)*G + G = 0');
          equal(G[1].multiply(CURVE_ORDER - 1n).add(G[2]), G[1], '(N-1)*G + 2*G = 1*G');
          equal(G[1].multiply(CURVE_ORDER - 2n).add(G[2]), G[0], '(N-2)*G + 2*G = 0');
@@ -109,7 +393,7 @@ for (const name in CURVES) {
          const carry = CURVE_ORDER % 2n === 1n ? G[1] : G[0];
          equal(G[1].multiply(half).double().add(carry), G[0], '((N/2) * G).double() = 0');
        });
-    should(`${name}/${pointName}/Basic group laws (inversion)`, () => {
+        should('(inversion)', () => {
          const a = 1234n;
          const b = 5678n;
          const c = a * b;
@@ -117,7 +401,7 @@ for (const name in CURVES) {
          const inv = mod.invert(b, CURVE_ORDER);
          equal(G[1].multiply(c).multiply(inv), G[1].multiply(a), 'c*G * (1/b)*G = a*G');
        });
-    should(`${name}/${pointName}/Basic group laws (multiply, rand)`, () =>
+        should('(multiply, rand)', () =>
          fc.assert(
            fc.property(FC_BIGINT, FC_BIGINT, (a, b) => {
              const c = mod.mod(a + b, CURVE_ORDER);
@@ -125,27 +409,29 @@ for (const name in CURVES) {
              const pA = G[1].multiply(a);
              const pB = G[1].multiply(b);
              const pC = G[1].multiply(c);
-          equal(pA.add(pB), pB.add(pA), `pA + pB = pB + pA`);
-          equal(pA.add(pB), pC, `pA + pB = pC`);
+              equal(pA.add(pB), pB.add(pA), 'pA + pB = pB + pA');
+              equal(pA.add(pB), pC, 'pA + pB = pC');
            }),
            { numRuns: NUM_RUNS }
          )
        );
-    should(`${name}/${pointName}/Basic group laws (multiply2, rand)`, () =>
+        should('(multiply2, rand)', () =>
          fc.assert(
            fc.property(FC_BIGINT, FC_BIGINT, (a, b) => {
              const c = mod.mod(a * b, CURVE_ORDER);
              const pA = G[1].multiply(a);
              const pB = G[1].multiply(b);
-          equal(pA.multiply(b), pB.multiply(a), `b*pA = a*pB`);
-          equal(pA.multiply(b), G[1].multiply(c), `b*pA = c*G`);
+              equal(pA.multiply(b), pB.multiply(a), 'b*pA = a*pB');
+              equal(pA.multiply(b), G[1].multiply(c), 'b*pA = c*G');
            }),
            { numRuns: NUM_RUNS }
          )
        );
+      });

      for (const op of ['add', 'subtract']) {
-      should(`${name}/${pointName}/${op} type check`, () => {
+        describe(op, () => {
+          should('type check', () => {
            throws(() => G[1][op](0), '0');
            throws(() => G[1][op](0n), '0n');
            G[1][op](G[2]);
@@ -153,17 +439,21 @@ for (const name in CURVES) {
            throws(() => G[1][op](123.456), '123.456');
            throws(() => G[1][op](true), 'true');
            throws(() => G[1][op]('1'), "'1'");
-        throws(() => G[1][op]({ x: 1n, y: 1n, z: 1n, t: 1n }), '{ x: 1n, y: 1n, z: 1n, t: 1n }');
+            throws(
+              () => G[1][op]({ x: 1n, y: 1n, z: 1n, t: 1n }),
+              '{ x: 1n, y: 1n, z: 1n, t: 1n }'
+            );
            throws(() => G[1][op](new Uint8Array([])), 'ui8a([])');
            throws(() => G[1][op](new Uint8Array([0])), 'ui8a([0])');
            throws(() => G[1][op](new Uint8Array([1])), 'ui8a([1])');
            throws(() => G[1][op](new Uint8Array(4096).fill(1)), 'ui8a(4096*[1])');
-        if (G[1].toAffine) throws(() => G[1][op](C.Point.BASE), `Point ${op} ${pointName}`);
-        throws(() => G[1][op](o.BASE), `${op}/other curve point`);
+            if (G[1].toAffine) throws(() => G[1][op](C.Point.BASE), 'Point ${op} ${pointName}');
+            throws(() => G[1][op](o.BASE), '${op}/other curve point');
+          });
        });
      }

-    should(`${name}/${pointName}/equals type check`, () => {
+      should('equals type check', () => {
        throws(() => G[1].equals(0), '0');
        throws(() => G[1].equals(0n), '0n');
        deepStrictEqual(G[1].equals(G[2]), false, '1*G != 2*G');
@@ -178,13 +468,14 @@ for (const name in CURVES) {
        throws(() => G[1].equals(new Uint8Array([0])), 'ui8a([0])');
        throws(() => G[1].equals(new Uint8Array([1])), 'ui8a([1])');
        throws(() => G[1].equals(new Uint8Array(4096).fill(1)), 'ui8a(4096*[1])');
-      if (G[1].toAffine) throws(() => G[1].equals(C.Point.BASE), `Point.equals(${pointName})`);
+        if (G[1].toAffine) throws(() => G[1].equals(C.Point.BASE), 'Point.equals(${pointName})');
        throws(() => G[1].equals(o.BASE), 'other curve point');
      });

      for (const op of ['multiply', 'multiplyUnsafe']) {
        if (!p.BASE[op]) continue;
-      should(`${name}/${pointName}/${op} type check`, () => {
+        describe(op, () => {
+          should('type check', () => {
            if (op !== 'multiplyUnsafe') {
              throws(() => G[1][op](0), '0');
              throws(() => G[1][op](0n), '0n');
@@ -203,23 +494,24 @@ for (const name in CURVES) {
            throws(() => G[1][op](new Uint8Array(4096).fill(1)), 'ui8a(4096*[1])');
            throws(() => G[1][op](o.BASE), 'other curve point');
          });
+        });
      }
      // Complex point (Extended/Jacobian/Projective?)
      if (p.BASE.toAffine) {
-      should(`${name}/${pointName}/toAffine()`, () => {
-        equal(p.ZERO.toAffine(), C.Point.ZERO, `0 = 0`);
-        equal(p.BASE.toAffine(), C.Point.BASE, `1 = 1`);
+        should('toAffine()', () => {
+          equal(p.ZERO.toAffine(), C.Point.ZERO, '0 = 0');
+          equal(p.BASE.toAffine(), C.Point.BASE, '1 = 1');
        });
      }
      if (p.fromAffine) {
-      should(`${name}/${pointName}/fromAffine()`, () => {
-        equal(p.ZERO, p.fromAffine(C.Point.ZERO), `0 = 0`);
-        equal(p.BASE, p.fromAffine(C.Point.BASE), `1 = 1`);
+        should('fromAffine()', () => {
+          equal(p.ZERO, p.fromAffine(C.Point.ZERO), '0 = 0');
+          equal(p.BASE, p.fromAffine(C.Point.BASE), '1 = 1');
        });
      }
      // toHex/fromHex (if available)
      if (p.fromHex && p.BASE.toHex) {
-      should(`${name}/${pointName}/fromHex(toHex()) roundtrip`, () => {
+        should('fromHex(toHex()) roundtrip', () => {
          fc.assert(
            fc.property(FC_BIGINT, (x) => {
              const hex = p.BASE.multiply(x).toHex();
@@ -228,9 +520,11 @@ for (const name in CURVES) {
          );
        });
      }
+    });
  }
+  describe(name, () => {
    // Generic complex things (getPublicKey/sign/verify/getSharedSecret)
-  should(`${name}/getPublicKey type check`, () => {
+    should('getPublicKey type check', () => {
      throws(() => C.getPublicKey(0), '0');
      throws(() => C.getPublicKey(0n), '0n');
      throws(() => C.getPublicKey(false), 'false');
@@ -245,23 +539,27 @@ for (const name in CURVES) {
      throws(() => C.getPublicKey(new Uint8Array([1])));
      throws(() => C.getPublicKey(new Uint8Array(4096).fill(1)));
    });
-  should(`${name}.verify()/should verify random signatures`, () =>
+    should('.verify() should verify random signatures', () =>
      fc.assert(
        fc.property(fc.hexaString({ minLength: 64, maxLength: 64 }), (msg) => {
          const priv = C.utils.randomPrivateKey();
          const pub = C.getPublicKey(priv);
          const sig = C.sign(msg, priv);
-        deepStrictEqual(C.verify(sig, msg, pub), true);
+          deepStrictEqual(
+            C.verify(sig, msg, pub),
+            true,
+            'priv=${toHex(priv)},pub=${toHex(pub)},msg=${msg}'
+          );
        }),
        { numRuns: NUM_RUNS }
      )
    );
-  should(`${name}.sign()/edge cases`, () => {
+    should('.sign() edge cases', () => {
      throws(() => C.sign());
      throws(() => C.sign(''));
    });

-  should(`${name}.verify()/should not verify signature with wrong hash`, () => {
+    should('.verify() should not verify signature with wrong hash', () => {
      const MSG = '01'.repeat(32);
      const PRIV_KEY = 0x2n;
      const WRONG_MSG = '11'.repeat(32);
@@ -271,7 +569,7 @@ for (const name in CURVES) {
    });
    // NOTE: fails for ed, because of empty message. Since we convert it to scalar,
    // need to check what other implementations do. Empty message != new Uint8Array([0]), but what scalar should be in that case?
-  // should(`${name}/should not verify signature with wrong message`, () => {
+    // should('should not verify signature with wrong message', () => {
    //   fc.assert(
    //     fc.property(
    //       fc.array(fc.integer({ min: 0x00, max: 0xff })),
@@ -293,7 +591,7 @@ for (const name in CURVES) {
    // });

    if (C.getSharedSecret) {
-    should(`${name}/getSharedSecret() should be commutative`, () => {
+      should('getSharedSecret() should be commutative', () => {
        for (let i = 0; i < NUM_RUNS; i++) {
          const asec = C.utils.randomPrivateKey();
          const apub = C.getPublicKey(asec);
@@ -308,7 +606,24 @@ for (const name in CURVES) {
        }
      });
    }
+  });
 }
+
+should('secp224k1 sqrt bug', () => {
+  const { Fp } = secp224r1.CURVE;
+  const sqrtMinus1 = Fp.sqrt(-1n);
+  // Verified against sage
+  deepStrictEqual(
+    sqrtMinus1,
+    23621584063597419797792593680131996961517196803742576047493035507225n
+  );
+  deepStrictEqual(
+    Fp.negate(sqrtMinus1),
+    3338362603553219996874421406887633712040719456283732096017030791656n
+  );
+  deepStrictEqual(Fp.square(sqrtMinus1), Fp.create(-1n));
+});
+
 // ESM is broken.
 import url from 'url';
 if (import.meta.url === url.pathToFileURL(process.argv[1]).href) {
--- a/test/bls12-381.test.js
+++ b/test/bls12-381.test.js
--- a/test/ed25519.test.js
+++ b/test/ed25519.test.js
@@ -1,7 +1,14 @@
-import { deepStrictEqual, throws } from 'assert';
-import { should } from 'micro-should';
+import { deepEqual, deepStrictEqual, strictEqual, throws } from 'assert';
+import { describe, should } from 'micro-should';
 import * as fc from 'fast-check';
-import { ed25519, ed25519ctx, ed25519ph, x25519, RistrettoPoint } from '../lib/esm/ed25519.js';
+import {
+  ed25519,
+  ed25519ctx,
+  ed25519ph,
+  x25519,
+  RistrettoPoint,
+  ED25519_TORSION_SUBGROUP,
+} from '../lib/esm/ed25519.js';
 import { readFileSync } from 'fs';
 import { default as zip215 } from './ed25519/zip215.json' assert { type: 'json' };
 import { hexToBytes, bytesToHex, randomBytes } from '@noble/hashes/utils';
@@ -10,6 +17,7 @@ import { sha512 } from '@noble/hashes/sha512';
 import { default as ed25519vectors } from './wycheproof/eddsa_test.json' assert { type: 'json' };
 import { default as x25519vectors } from './wycheproof/x25519_test.json' assert { type: 'json' };

+describe('ed25519', () => {
  const ed = ed25519;
  const hex = bytesToHex;

@@ -27,12 +35,12 @@ function utf8ToBytes(str) {

  ed.utils.precompute(8);

-should('ed25519/should not accept >32byte private keys', () => {
+  should('not accept >32byte private keys', () => {
    const invalidPriv =
      100000000000000000000000000000000000009000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000090000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000800073278156000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000n;
    throws(() => ed.getPublicKey(invalidPriv));
  });
-should('ed25519/should verify recent signature', () => {
+  should('verify recent signature', () => {
    fc.assert(
      fc.property(
        fc.hexaString({ minLength: 2, maxLength: 32 }),
@@ -48,7 +56,7 @@ should('ed25519/should verify recent signature', () => {
      { numRuns: 5 }
    );
  });
-should('ed25519/should not verify signature with wrong message', () => {
+  should('not verify signature with wrong message', () => {
    fc.assert(
      fc.property(
        fc.array(fc.integer({ min: 0x00, max: 0xff })),
@@ -72,38 +80,40 @@ should('ed25519/should not verify signature with wrong message', () => {
  const privKey = to32Bytes('a665a45920422f9d417e4867ef');
  const msg = hexToBytes('874f9960c5d2b7a9b5fad383e1ba44719ebb743a');
  const wrongMsg = hexToBytes('589d8c7f1da0a24bc07b7381ad48b1cfc211af1c');
-should('ed25519/basic methods/should sign and verify', () => {
+  describe('basic methods', () => {
+    should('sign and verify', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), true);
    });
-should('ed25519/basic methods/should not verify signature with wrong public key', () => {
+    should('not verify signature with wrong public key', () => {
      const publicKey = ed.getPublicKey(12);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), false);
    });
-should('ed25519/basic methods/should not verify signature with wrong hash', () => {
+    should('not verify signature with wrong hash', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, wrongMsg, publicKey), false);
    });
-
-should('ed25519/sync methods/should sign and verify', () => {
+  });
+  describe('sync methods', () => {
+    should('sign and verify', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), true);
    });
-should('ed25519/sync methods/should not verify signature with wrong public key', () => {
+    should('not verify signature with wrong public key', () => {
      const publicKey = ed.getPublicKey(12);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), false);
    });
-should('ed25519/sync methods/should not verify signature with wrong hash', () => {
+    should('not verify signature with wrong hash', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, wrongMsg, publicKey), false);
    });
-
+  });
  // https://xmr.llcoins.net/addresstests.html
  should(
    'ed25519/BASE_POINT.multiply()/should create right publicKey without SHA-512 hashing TEST 1',
@@ -646,10 +656,23 @@ for (let i = 0; i < VECTORS_RFC8032_PH.length; i++) {

  should('X25519 base point', () => {
    const { y } = ed25519.Point.BASE;
-  const u = ed25519.utils.mod((y + 1n) * ed25519.utils.invert(1n - y, ed25519.CURVE.P));
+    const { Fp } = ed25519.CURVE;
+    const u = Fp.create((y + 1n) * Fp.invert(1n - y));
    deepStrictEqual(hex(numberToBytesLE(u, 32)), x25519.Gu);
  });

+  should('isTorsionFree()', () => {
+    const orig = ed.utils.getExtendedPublicKey(ed.utils.randomPrivateKey()).point;
+    for (const hex of ED25519_TORSION_SUBGROUP.slice(1)) {
+      const dirty = orig.add(ed.Point.fromHex(hex));
+      const cleared = dirty.clearCofactor();
+      strictEqual(orig.isTorsionFree(), true, `orig must be torsionFree: ${hex}`);
+      strictEqual(dirty.isTorsionFree(), false, `dirty must not be torsionFree: ${hex}`);
+      strictEqual(cleared.isTorsionFree(), true, `cleared must be torsionFree: ${hex}`);
+    }
+  });
+});
+
 // ESM is broken.
 import url from 'url';
 if (import.meta.url === url.pathToFileURL(process.argv[1]).href) {
--- a/test/ed448.test.js
+++ b/test/ed448.test.js
@@ -1,5 +1,5 @@
 import { deepStrictEqual, throws } from 'assert';
-import { should } from 'micro-should';
+import { describe, should } from 'micro-should';
 import * as fc from 'fast-check';
 import { ed448, ed448ph, x448 } from '../lib/esm/ed448.js';
 import { hexToBytes, bytesToHex, randomBytes } from '@noble/hashes/utils';
@@ -7,6 +7,7 @@ import { numberToBytesLE } from '../lib/esm/abstract/utils.js';
 import { default as ed448vectors } from './wycheproof/ed448_test.json' assert { type: 'json' };
 import { default as x448vectors } from './wycheproof/x448_test.json' assert { type: 'json' };

+describe('ed448', () => {
  const ed = ed448;
  const hex = bytesToHex;
  ed.utils.precompute(4);
@@ -323,7 +324,7 @@ for (let i = 0; i < VECTORS_RFC8032.length; i++) {
    });
  }

-should('ed448/should not accept >57byte private keys', async () => {
+  should('not accept >57byte private keys', async () => {
    const invalidPriv =
      100000000000000000000000000000000000009000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000090000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000800073278156000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000n;
    throws(() => ed.getPublicKey(invalidPriv));
@@ -334,7 +335,7 @@ function to57Bytes(numOrStr) {
    return hexToBytes(hex.padStart(114, '0'));
  }

-should('ed448/should verify recent signature', () => {
+  should('verify recent signature', () => {
    fc.assert(
      fc.property(
        fc.hexaString({ minLength: 2, maxLength: 57 }),
@@ -350,7 +351,7 @@ should('ed448/should verify recent signature', () => {
      { numRuns: 5 }
    );
  });
-should('ed448/should not verify signature with wrong message', () => {
+  should('not verify signature with wrong message', () => {
    fc.assert(
      fc.property(
        fc.array(fc.integer({ min: 0x00, max: 0xff })),
@@ -374,39 +375,43 @@ should('ed448/should not verify signature with wrong message', () => {
  const privKey = to57Bytes('a665a45920422f9d417e4867ef');
  const msg = hexToBytes('874f9960c5d2b7a9b5fad383e1ba44719ebb743a');
  const wrongMsg = hexToBytes('589d8c7f1da0a24bc07b7381ad48b1cfc211af1c');
-should('ed25519/basic methods/should sign and verify', () => {
+  describe('basic methods', () => {
+    should('sign and verify', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), true);
    });
-should('ed25519/basic methods/should not verify signature with wrong public key', () => {
+    should('not verify signature with wrong public key', () => {
      const publicKey = ed.getPublicKey(12);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), false);
    });
-should('ed25519/basic methods/should not verify signature with wrong hash', () => {
+    should('not verify signature with wrong hash', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, wrongMsg, publicKey), false);
    });
+  });

-should('ed25519/sync methods/should sign and verify', () => {
+  describe('sync methods', () => {
+    should('sign and verify', () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), true);
    });
-should('ed25519/sync methods/should not verify signature with wrong public key', async () => {
+    should('not verify signature with wrong public key', async () => {
      const publicKey = ed.getPublicKey(12);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, msg, publicKey), false);
    });
-should('ed25519/sync methods/should not verify signature with wrong hash', async () => {
+    should('not verify signature with wrong hash', async () => {
      const publicKey = ed.getPublicKey(privKey);
      const signature = ed.sign(msg, privKey);
      deepStrictEqual(ed.verify(signature, wrongMsg, publicKey), false);
    });
+  });

-should('ed25519/BASE_POINT.multiply()/should throw Point#multiply on TEST 5', () => {
+  should('BASE_POINT.multiply() throws in Point#multiply on TEST 5', () => {
    for (const num of [0n, 0, -1n, -1, 1.1]) {
      throws(() => ed.Point.BASE.multiply(num));
    }
@@ -651,11 +656,13 @@ for (let i = 0; i < VECTORS_RFC8032_PH.length; i++) {

  should('X448 base point', () => {
    const { x, y } = ed448.Point.BASE;
-  const { P } = ed448.CURVE;
-  const invX = ed448.utils.invert(x * x, P); // x²
-  const u = ed448.utils.mod(y * y * invX, P); // (y²/x²)
+    const { Fp } = ed448.CURVE;
+    // const invX = Fp.invert(x * x); // x²
+    const u = Fp.div(Fp.create(y * y), Fp.create(x * x)); // (y²/x²)
+    // const u = Fp.create(y * y * invX);
    deepStrictEqual(hex(numberToBytesLE(u, 56)), x448.Gu);
  });
+});

 // ESM is broken.
 import url from 'url';
--- a/test/hash-to-curve.test.js
+++ b/test/hash-to-curve.test.js
@@ -1,19 +1,30 @@
 import { deepStrictEqual } from 'assert';
-import { should } from 'micro-should';
+import { describe, should } from 'micro-should';
 import { bytesToHex } from '@noble/hashes/utils';
 // Generic tests for all curves in package
 import { sha256 } from '@noble/hashes/sha256';
 import { sha512 } from '@noble/hashes/sha512';
+import { shake128, shake256 } from '@noble/hashes/sha3';
 import { secp256r1 } from '../lib/esm/p256.js';
 import { secp384r1 } from '../lib/esm/p384.js';
 import { secp521r1 } from '../lib/esm/p521.js';
+import { ed25519 } from '../lib/esm/ed25519.js';
+import { ed448 } from '../lib/esm/ed448.js';
 import { secp256k1 } from '../lib/esm/secp256k1.js';
 import { bls12_381 } from '../lib/esm/bls12-381.js';
-import { stringToBytes, expand_message_xmd } from '../lib/esm/abstract/hash-to-curve.js';
-
+import {
+  stringToBytes,
+  expand_message_xmd,
+  expand_message_xof,
+} from '../lib/esm/abstract/hash-to-curve.js';
+// XMD
 import { default as xmd_sha256_38 } from './hash-to-curve/expand_message_xmd_SHA256_38.json' assert { type: 'json' };
 import { default as xmd_sha256_256 } from './hash-to-curve/expand_message_xmd_SHA256_256.json' assert { type: 'json' };
 import { default as xmd_sha512_38 } from './hash-to-curve/expand_message_xmd_SHA512_38.json' assert { type: 'json' };
+// XOF
+import { default as xof_shake128_36 } from './hash-to-curve/expand_message_xof_SHAKE128_36.json' assert { type: 'json' };
+import { default as xof_shake128_256 } from './hash-to-curve/expand_message_xof_SHAKE128_256.json' assert { type: 'json' };
+import { default as xof_shake256_36 } from './hash-to-curve/expand_message_xof_SHAKE256_36.json' assert { type: 'json' };
 // P256
 import { default as p256_ro } from './hash-to-curve/P256_XMD:SHA-256_SSWU_RO_.json' assert { type: 'json' };
 import { default as p256_nu } from './hash-to-curve/P256_XMD:SHA-256_SSWU_NU_.json' assert { type: 'json' };
@@ -40,9 +51,10 @@ import { default as ed448_ro } from './hash-to-curve/edwards448_XOF:SHAKE256_ELL
 import { default as ed448_nu } from './hash-to-curve/edwards448_XOF:SHAKE256_ELL2_NU_.json' assert { type: 'json' };

 function testExpandXMD(hash, vectors) {
+  describe(`${vectors.hash}/${vectors.DST.length}`, () => {
    for (let i = 0; i < vectors.tests.length; i++) {
      const t = vectors.tests[i];
-    should(`expand_message_xmd/${vectors.hash}/${vectors.DST.length}/${i}`, () => {
+      should(`${vectors.hash}/${vectors.DST.length}/${i}`, () => {
        const p = expand_message_xmd(
          stringToBytes(t.msg),
          stringToBytes(vectors.DST),
@@ -52,11 +64,38 @@ function testExpandXMD(hash, vectors) {
        deepStrictEqual(bytesToHex(p), t.uniform_bytes);
      });
    }
+  });
 }

+describe('expand_message_xmd', () => {
  testExpandXMD(sha256, xmd_sha256_38);
  testExpandXMD(sha256, xmd_sha256_256);
  testExpandXMD(sha512, xmd_sha512_38);
+});
+
+function testExpandXOF(hash, vectors) {
+  describe(`${vectors.hash}/${vectors.DST.length}`, () => {
+    for (let i = 0; i < vectors.tests.length; i++) {
+      const t = vectors.tests[i];
+      should(`${i}`, () => {
+        const p = expand_message_xof(
+          stringToBytes(t.msg),
+          stringToBytes(vectors.DST),
+          +t.len_in_bytes,
+          vectors.k,
+          hash
+        );
+        deepStrictEqual(bytesToHex(p), t.uniform_bytes);
+      });
+    }
+  });
+}
+
+describe('expand_message_xof', () => {
+  testExpandXOF(shake128, xof_shake128_36);
+  testExpandXOF(shake128, xof_shake128_256);
+  testExpandXOF(shake256, xof_shake256_36);
+});

 function stringToFp(s) {
  // bls-G2 support
@@ -68,9 +107,10 @@ function stringToFp(s) {
 }

 function testCurve(curve, ro, nu) {
+  describe(`${ro.curve}/${ro.ciphersuite}`, () => {
    for (let i = 0; i < ro.vectors.length; i++) {
      const t = ro.vectors[i];
-    should(`${ro.curve}/${ro.ciphersuite}(${i})`, () => {
+      should(`(${i})`, () => {
        const p = curve.Point.hashToCurve(stringToBytes(t.msg), {
          DST: ro.dst,
        });
@@ -78,9 +118,11 @@ function testCurve(curve, ro, nu) {
        deepStrictEqual(p.y, stringToFp(t.P.y), 'Py');
      });
    }
+  });
+  describe(`${nu.curve}/${nu.ciphersuite}`, () => {
    for (let i = 0; i < nu.vectors.length; i++) {
      const t = nu.vectors[i];
-    should(`${nu.curve}/${nu.ciphersuite}(${i})`, () => {
+      should(`(${i})`, () => {
        const p = curve.Point.encodeToCurve(stringToBytes(t.msg), {
          DST: nu.dst,
        });
@@ -88,6 +130,7 @@ function testCurve(curve, ro, nu) {
        deepStrictEqual(p.y, stringToFp(t.P.y), 'Py');
      });
    }
+  });
 }

 testCurve(secp256r1, p256_ro, p256_nu);
@@ -97,8 +140,8 @@ testCurve(secp521r1, p521_ro, p521_nu);
 testCurve(bls12_381.G1, g1_ro, g1_nu);
 testCurve(bls12_381.G2, g2_ro, g2_nu);
 testCurve(secp256k1, secp256k1_ro, secp256k1_nu);
-//testCurve(ed25519, ed25519_ro, ed25519_nu);
-//testCurve(ed448, ed448_ro, ed448_nu);
+testCurve(ed25519, ed25519_ro, ed25519_nu);
+testCurve(ed448, ed448_ro, ed448_nu);

 // ESM is broken.
 import url from 'url';
--- a/test/jubjub.test.js
+++ b/test/jubjub.test.js
@@ -1,5 +1,5 @@
 import { jubjub, findGroupHash } from '../lib/esm/jubjub.js';
-import { should } from 'micro-should';
+import { describe, should } from 'micro-should';
 import { deepStrictEqual, throws } from 'assert';
 import { hexToBytes, bytesToHex } from '@noble/hashes/utils';

@@ -18,6 +18,7 @@ const G_PROOF = new jubjub.ExtendedPoint(

 const getXY = (p) => ({ x: p.x, y: p.y });

+describe('jubjub', () => {
  should('toHex/fromHex', () => {
    // More than field
    throws(() =>
@@ -32,8 +33,8 @@ should('toHex/fromHex', () => {
    throws(() =>
      jubjub.Point.fromHex(
        new Uint8Array([
-        7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-        0,
+          7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0,
        ])
      )
    );
@@ -61,11 +62,18 @@ should('toHex/fromHex', () => {
  });

  should('Find generators', () => {
-  const spend = findGroupHash(new Uint8Array(), new Uint8Array([90, 99, 97, 115, 104, 95, 71, 95]));
-  const proof = findGroupHash(new Uint8Array(), new Uint8Array([90, 99, 97, 115, 104, 95, 72, 95]));
+    const spend = findGroupHash(
+      new Uint8Array(),
+      new Uint8Array([90, 99, 97, 115, 104, 95, 71, 95])
+    );
+    const proof = findGroupHash(
+      new Uint8Array(),
+      new Uint8Array([90, 99, 97, 115, 104, 95, 72, 95])
+    );
    deepStrictEqual(getXY(spend.toAffine()), getXY(G_SPEND.toAffine()));
    deepStrictEqual(getXY(proof.toAffine()), getXY(G_PROOF.toAffine()));
  });
+});

 // ESM is broken.
 import url from 'url';
--- a/test/nist.test.js
+++ b/test/nist.test.js
@@ -1,5 +1,5 @@
 import { deepStrictEqual, throws } from 'assert';
-import { should } from 'micro-should';
+import { describe, should } from 'micro-should';
 import { secp192r1, P192 } from '../lib/esm/p192.js';
 import { secp224r1, P224 } from '../lib/esm/p224.js';
 import { secp256r1, P256 } from '../lib/esm/p256.js';
@@ -344,10 +344,11 @@ function runWycheproof(name, CURVE, group, index) {

 for (const name in WYCHEPROOF_ECDSA) {
  const { curve, hashes } = WYCHEPROOF_ECDSA[name];
+  describe('Wycheproof/WYCHEPROOF_ECDSA', () => {
    for (const hName in hashes) {
      const { hash, tests } = hashes[hName];
      const CURVE = curve.create(hash);
-    should(`Wycheproof/WYCHEPROOF_ECDSA ${name}/${hName}`, () => {
+      should(`${name}/${hName}`, () => {
        for (let i = 0; i < tests.length; i++) {
          const groups = tests[i].testGroups;
          for (let j = 0; j < groups.length; j++) {
@@ -357,6 +358,7 @@ for (const name in WYCHEPROOF_ECDSA) {
        }
      });
    }
+  });
 }

 const hexToBigint = (hex) => BigInt(`0x${hex}`);
--- a/test/secp256k1.test.js
+++ b/test/secp256k1.test.js
@@ -1,12 +1,13 @@
 import * as fc from 'fast-check';
 import { secp256k1, schnorr } from '../lib/esm/secp256k1.js';
+import { Fp } from '../lib/esm/abstract/modular.js';
 import { readFileSync } from 'fs';
 import { default as ecdsa } from './vectors/ecdsa.json' assert { type: 'json' };
 import { default as ecdh } from './vectors/ecdh.json' assert { type: 'json' };
 import { default as privates } from './vectors/privates.json' assert { type: 'json' };
 import { default as points } from './vectors/points.json' assert { type: 'json' };
 import { default as wp } from './vectors/wychenproof.json' assert { type: 'json' };
-import { should } from 'micro-should';
+import { should, describe } from 'micro-should';
 import { deepStrictEqual, throws } from 'assert';
 import { hexToBytes, bytesToHex } from '@noble/hashes/utils';

@@ -16,7 +17,6 @@ const privatesTxt = readFileSync('./test/vectors/privates-2.txt', 'utf-8');
 const schCsv = readFileSync('./test/vectors/schnorr.csv', 'utf-8');

 const FC_BIGINT = fc.bigInt(1n + 1n, secp.CURVE.n - 1n);
-const P = secp.CURVE.Fp.ORDER;
 // prettier-ignore
 const INVALID_ITEMS = ['deadbeef', Math.pow(2, 53), [1], 'xyzxyzxyxyzxyzxyxyzxyzxyxyzxyzxyxyzxyzxyxyzxyzxyxyzxyzxyxyzxyzxy', secp.CURVE.n + 2n];

@@ -30,7 +30,8 @@ function hexToNumber(hex) {
  return BigInt(`0x${hex}`);
 }

-should('secp256k1.getPublicKey()', () => {
+describe('secp256k1', () => {
+  should('getPublicKey()', () => {
    const data = privatesTxt
      .split('\n')
      .filter((line) => line)
@@ -49,12 +50,12 @@ should('secp256k1.getPublicKey()', () => {
      deepStrictEqual(toBEHex(point3.y), y);
    }
  });
-should('secp256k1.getPublicKey() rejects invalid keys', () => {
-  // for (const item of INVALID_ITEMS) {
-  //   throws(() => secp.getPublicKey(item));
-  // }
+  should('getPublicKey() rejects invalid keys', () => {
+    for (const item of INVALID_ITEMS) {
+      throws(() => secp.getPublicKey(item));
+    }
  });
-should('secp256k1.precompute', () => {
+  should('precompute', () => {
    secp.utils.precompute(4);
    const data = privatesTxt
      .split('\n')
@@ -75,7 +76,7 @@ should('secp256k1.precompute', () => {
    }
  });

-should('secp256k1.Point.isValidPoint()', () => {
+  should('Point.isValidPoint()', () => {
    for (const vector of points.valid.isPoint) {
      const { P, expected } = vector;
      if (expected) {
@@ -86,7 +87,7 @@ should('secp256k1.Point.isValidPoint()', () => {
    }
  });

-should('secp256k1.Point.fromPrivateKey()', () => {
+  should('Point.fromPrivateKey()', () => {
    for (const vector of points.valid.pointFromScalar) {
      const { d, expected } = vector;
      let p = secp.Point.fromPrivateKey(d);
@@ -94,7 +95,7 @@ should('secp256k1.Point.fromPrivateKey()', () => {
    }
  });

-should('secp256k1.Point#toHex(compressed)', () => {
+  should('Point#toHex(compressed)', () => {
    for (const vector of points.valid.pointCompress) {
      const { P, compress, expected } = vector;
      let p = secp.Point.fromHex(P);
@@ -102,7 +103,7 @@ should('secp256k1.Point#toHex(compressed)', () => {
    }
  });

-should('secp256k1.Point#toHex() roundtrip (failed case)', () => {
+  should('Point#toHex() roundtrip (failed case)', () => {
    const point1 =
      secp.Point.fromPrivateKey(
        88572218780422190464634044548753414301110513745532121983949500266768436236425n
@@ -111,7 +112,7 @@ should('secp256k1.Point#toHex() roundtrip (failed case)', () => {
    // deepStrictEqual(secp.Point.fromHex(hex).toHex(true), hex);
  });

-should('secp256k1.Point#toHex() roundtrip', () => {
+  should('Point#toHex() roundtrip', () => {
    fc.assert(
      fc.property(FC_BIGINT, (x) => {
        const point1 = secp.Point.fromPrivateKey(x);
@@ -121,7 +122,7 @@ should('secp256k1.Point#toHex() roundtrip', () => {
    );
  });

-should('secp256k1.Point#add(other)', () => {
+  should('Point#add(other)', () => {
    for (const vector of points.valid.pointAdd) {
      const { P, Q, expected } = vector;
      let p = secp.Point.fromHex(P);
@@ -136,7 +137,7 @@ should('secp256k1.Point#add(other)', () => {
    }
  });

-should('secp256k1.Point#multiply(privateKey)', () => {
+  should('Point#multiply(privateKey)', () => {
    for (const vector of points.valid.pointMultiply) {
      const { P, d, expected } = vector;
      const p = secp.Point.fromHex(P);
@@ -175,7 +176,7 @@ should('secp256k1.Point#multiply(privateKey)', () => {
  //   console.log(secp.ProjectivePoint.normalizeZ([p0.multiplyUnsafe(z)])[0])
  // });

-should('secp256k1.Signature.fromCompactHex() roundtrip', () => {
+  should('Signature.fromCompactHex() roundtrip', () => {
    fc.assert(
      fc.property(FC_BIGINT, FC_BIGINT, (r, s) => {
        const sig = new secp.Signature(r, s);
@@ -184,7 +185,7 @@ should('secp256k1.Signature.fromCompactHex() roundtrip', () => {
    );
  });

-should('secp256k1.Signature.fromDERHex() roundtrip', () => {
+  should('Signature.fromDERHex() roundtrip', () => {
    fc.assert(
      fc.property(FC_BIGINT, FC_BIGINT, (r, s) => {
        const sig = new secp.Signature(r, s);
@@ -193,7 +194,7 @@ should('secp256k1.Signature.fromDERHex() roundtrip', () => {
    );
  });

-should('secp256k1.sign()/should create deterministic signatures with RFC 6979', () => {
+  should('sign()/should create deterministic signatures with RFC 6979', () => {
    for (const vector of ecdsa.valid) {
      let usig = secp.sign(vector.m, vector.d);
      let sig = usig.toCompactHex();
@@ -203,18 +204,21 @@ should('secp256k1.sign()/should create deterministic signatures with RFC 6979',
    }
  });

-should('secp256k1.sign()/should not create invalid deterministic signatures with RFC 6979', () => {
+  should(
+    'secp256k1.sign()/should not create invalid deterministic signatures with RFC 6979',
+    () => {
      for (const vector of ecdsa.invalid.sign) {
        throws(() => secp.sign(vector.m, vector.d));
      }
-});
+    }
+  );

-should('secp256k1.sign()/edge cases', () => {
+  should('sign()/edge cases', () => {
    throws(() => secp.sign());
    throws(() => secp.sign(''));
  });

-should('secp256k1.sign()/should create correct DER encoding against libsecp256k1', () => {
+  should('sign()/should create correct DER encoding against libsecp256k1', () => {
    const CASES = [
      [
        'd1a9dc8ed4e46a6a3e5e594615ca351d7d7ef44df1e4c94c1802f3592183794b',
@@ -237,7 +241,7 @@ should('secp256k1.sign()/should create correct DER encoding against libsecp256k1
      deepStrictEqual(secp.Signature.fromCompact(rs).toDERHex(), exp);
    }
  });
-should('secp256k1.sign()/sign ecdsa extraData', () => {
+  should('sign()/sign ecdsa extraData', () => {
    const ent1 = '0000000000000000000000000000000000000000000000000000000000000000';
    const ent2 = '0000000000000000000000000000000000000000000000000000000000000001';
    const ent3 = '6e723d3fd94ed5d2b6bdd4f123364b0f3ca52af829988a63f8afe91d29db1c33';
@@ -258,7 +262,7 @@ should('secp256k1.sign()/sign ecdsa extraData', () => {
    }
  });

-should('secp256k1.verify()/should verify signature', () => {
+  should('verify()/should verify signature', () => {
    const MSG = '01'.repeat(32);
    const PRIV_KEY = 0x2n;
    const signature = secp.sign(MSG, PRIV_KEY);
@@ -266,7 +270,7 @@ should('secp256k1.verify()/should verify signature', () => {
    deepStrictEqual(publicKey.length, 65);
    deepStrictEqual(secp.verify(signature, MSG, publicKey), true);
  });
-should('secp256k1.verify()/should not verify signature with wrong public key', () => {
+  should('verify()/should not verify signature with wrong public key', () => {
    const MSG = '01'.repeat(32);
    const PRIV_KEY = 0x2n;
    const WRONG_PRIV_KEY = 0x22n;
@@ -275,7 +279,7 @@ should('secp256k1.verify()/should not verify signature with wrong public key', (
    deepStrictEqual(publicKey.length, 130);
    deepStrictEqual(secp.verify(signature, MSG, publicKey), false);
  });
-should('secp256k1.verify()/should not verify signature with wrong hash', () => {
+  should('verify()/should not verify signature with wrong hash', () => {
    const MSG = '01'.repeat(32);
    const PRIV_KEY = 0x2n;
    const WRONG_MSG = '11'.repeat(32);
@@ -284,7 +288,7 @@ should('secp256k1.verify()/should not verify signature with wrong hash', () => {
    deepStrictEqual(publicKey.length, 65);
    deepStrictEqual(secp.verify(signature, WRONG_MSG, publicKey), false);
  });
-should('secp256k1.verify()/should verify random signatures', () =>
+  should('verify()/should verify random signatures', () =>
    fc.assert(
      fc.property(FC_BIGINT, fc.hexaString({ minLength: 64, maxLength: 64 }), (privKey, msg) => {
        const pub = secp.getPublicKey(privKey);
@@ -293,10 +297,11 @@ should('secp256k1.verify()/should verify random signatures', () =>
      })
    )
  );
-should('secp256k1.verify()/should not verify signature with invalid r/s', () => {
+  should('verify()/should not verify signature with invalid r/s', () => {
    const msg = new Uint8Array([
-    0xbb, 0x5a, 0x52, 0xf4, 0x2f, 0x9c, 0x92, 0x61, 0xed, 0x43, 0x61, 0xf5, 0x94, 0x22, 0xa1, 0xe3,
-    0x00, 0x36, 0xe7, 0xc3, 0x2b, 0x27, 0x0c, 0x88, 0x07, 0xa4, 0x19, 0xfe, 0xca, 0x60, 0x50, 0x23,
+      0xbb, 0x5a, 0x52, 0xf4, 0x2f, 0x9c, 0x92, 0x61, 0xed, 0x43, 0x61, 0xf5, 0x94, 0x22, 0xa1,
+      0xe3, 0x00, 0x36, 0xe7, 0xc3, 0x2b, 0x27, 0x0c, 0x88, 0x07, 0xa4, 0x19, 0xfe, 0xca, 0x60,
+      0x50, 0x23,
    ]);
    const x = 100260381870027870612475458630405506840396644859280795015145920502443964769584n;
    const y = 41096923727651821103518389640356553930186852801619204169823347832429067794568n;
@@ -312,7 +317,7 @@ should('secp256k1.verify()/should not verify signature with invalid r/s', () =>
    // Verifies, but it shouldn't, because signature S > curve order
    deepStrictEqual(verified, false);
  });
-should('secp256k1.verify()/should not verify msg = curve order', () => {
+  should('verify()/should not verify msg = curve order', () => {
    const msg = 'fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141';
    const x = 55066263022277343669578718895168534326250603453777594175500187360389116729240n;
    const y = 32670510020758816978083085130507043184471273380659243275938904335757337482424n;
@@ -322,7 +327,7 @@ should('secp256k1.verify()/should not verify msg = curve order', () => {
    const sig = new secp.Signature(r, s);
    deepStrictEqual(secp.verify(sig, msg, pub), false);
  });
-should('secp256k1.verify()/should verify non-strict msg bb5a...', () => {
+  should('verify()/should verify non-strict msg bb5a...', () => {
    const msg = 'bb5a52f42f9c9261ed4361f59422a1e30036e7c32b270c8807a419feca605023';
    const x = 3252872872578928810725465493269682203671229454553002637820453004368632726370n;
    const y = 17482644437196207387910659778872952193236850502325156318830589868678978890912n;
@@ -362,7 +367,7 @@ for (let vec of vectors) {
    });
  }

-should('secp256k1.recoverPublicKey()/should recover public key from recovery bit', () => {
+  should('recoverPublicKey()/should recover public key from recovery bit', () => {
    const message = '00000000000000000000000000000000000000000000000000000000deadbeef';
    const privateKey = 123456789n;
    const publicKey = secp.Point.fromHex(secp.getPublicKey(privateKey)).toHex(false);
@@ -370,17 +375,17 @@ should('secp256k1.recoverPublicKey()/should recover public key from recovery bit
    const recoveredPubkey = sig.recoverPublicKey(message);
    // const recoveredPubkey = secp.recoverPublicKey(message, signature, recovery);
    deepStrictEqual(recoveredPubkey !== null, true);
-  deepStrictEqual(recoveredPubkey.toHex(), publicKey);
+    deepStrictEqual(recoveredPubkey.toHex(false), publicKey);
    deepStrictEqual(secp.verify(sig, message, publicKey), true);
  });
-should('secp256k1.recoverPublicKey()/should not recover zero points', () => {
+  should('recoverPublicKey()/should not recover zero points', () => {
    const msgHash = '6b8d2c81b11b2d699528dde488dbdf2f94293d0d33c32e347f255fa4a6c1f0a9';
    const sig =
      '79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f817986b8d2c81b11b2d699528dde488dbdf2f94293d0d33c32e347f255fa4a6c1f0a9';
    const recovery = 0;
    throws(() => secp.recoverPublicKey(msgHash, sig, recovery));
  });
-should('secp256k1.recoverPublicKey()/should handle all-zeros msghash', () => {
+  should('recoverPublicKey()/should handle all-zeros msghash', () => {
    const privKey = secp.utils.randomPrivateKey();
    const pub = secp.getPublicKey(privKey);
    const zeros = '0000000000000000000000000000000000000000000000000000000000000000';
@@ -388,7 +393,7 @@ should('secp256k1.recoverPublicKey()/should handle all-zeros msghash', () => {
    const recoveredKey = sig.recoverPublicKey(zeros);
    deepStrictEqual(recoveredKey.toRawBytes(), pub);
  });
-should('secp256k1.recoverPublicKey()/should handle RFC 6979 vectors', () => {
+  should('recoverPublicKey()/should handle RFC 6979 vectors', () => {
    for (const vector of ecdsa.valid) {
      let usig = secp.sign(vector.m, vector.d);
      let sig = usig.toDERHex();
@@ -402,7 +407,7 @@ should('secp256k1.recoverPublicKey()/should handle RFC 6979 vectors', () => {
  function derToPub(der) {
    return der.slice(46);
  }
-should('secp256k1.getSharedSecret()/should produce correct results', () => {
+  should('getSharedSecret()/should produce correct results', () => {
    // TODO: Once der is there, run all tests.
    for (const vector of ecdh.testGroups[0].tests.slice(0, 230)) {
      if (vector.result === 'invalid' || vector.private.length !== 64) {
@@ -415,7 +420,7 @@ should('secp256k1.getSharedSecret()/should produce correct results', () => {
      }
    }
  });
-should('secp256k1.getSharedSecret()/priv/pub order matters', () => {
+  should('getSharedSecret()/priv/pub order matters', () => {
    for (const vector of ecdh.testGroups[0].tests.slice(0, 100)) {
      if (vector.result === 'valid') {
        let priv = vector.private;
@@ -424,27 +429,36 @@ should('secp256k1.getSharedSecret()/priv/pub order matters', () => {
      }
    }
  });
-should('secp256k1.getSharedSecret()/rejects invalid keys', () => {
+  should('getSharedSecret()/rejects invalid keys', () => {
    throws(() => secp.getSharedSecret('01', '02'));
  });

-should('secp256k1.utils.isValidPrivateKey()', () => {
+  should('utils.isValidPrivateKey()', () => {
    for (const vector of privates.valid.isPrivate) {
      const { d, expected } = vector;
      deepStrictEqual(secp.utils.isValidPrivateKey(d), expected);
    }
  });
+  should('have proper curve equation in assertValidity()', () => {
+    throws(() => {
+      const { Fp } = secp.CURVE;
+      let point = new secp.Point(Fp.create(-2n), Fp.create(-1n));
+      point.assertValidity();
+    });
+  });
+
+  const Fn = Fp(secp.CURVE.n);
  const normal = secp.utils._normalizePrivateKey;
  const tweakUtils = {
    privateAdd: (privateKey, tweak) => {
      const p = normal(privateKey);
      const t = normal(tweak);
-    return secp.utils._bigintToBytes(secp.utils.mod(p + t, secp.CURVE.n));
+      return secp.utils._bigintToBytes(Fn.create(p + t));
    },

    privateNegate: (privateKey) => {
      const p = normal(privateKey);
-    return secp.utils._bigintToBytes(secp.CURVE.n - p);
+      return secp.utils._bigintToBytes(Fn.negate(p));
    },

    pointAddScalar: (p, tweak, isCompressed) => {
@@ -463,45 +477,45 @@ const tweakUtils = {
    },
  };

-should('secp256k1.privateAdd()', () => {
+  should('privateAdd()', () => {
    for (const vector of privates.valid.add) {
      const { a, b, expected } = vector;
      deepStrictEqual(bytesToHex(tweakUtils.privateAdd(a, b)), expected);
    }
  });
-should('secp256k1.privateNegate()', () => {
+  should('privateNegate()', () => {
    for (const vector of privates.valid.negate) {
      const { a, expected } = vector;
      deepStrictEqual(bytesToHex(tweakUtils.privateNegate(a)), expected);
    }
  });
-should('secp256k1.pointAddScalar()', () => {
+  should('pointAddScalar()', () => {
    for (const vector of points.valid.pointAddScalar) {
      const { description, P, d, expected } = vector;
      const compressed = !!expected && expected.length === 66; // compressed === 33 bytes
      deepStrictEqual(bytesToHex(tweakUtils.pointAddScalar(P, d, compressed)), expected);
    }
  });
-should('secp256k1.pointAddScalar() invalid', () => {
+  should('pointAddScalar() invalid', () => {
    for (const vector of points.invalid.pointAddScalar) {
      const { P, d, exception } = vector;
      throws(() => tweakUtils.pointAddScalar(P, d));
    }
  });
-should('secp256k1.pointMultiply()', () => {
+  should('pointMultiply()', () => {
    for (const vector of points.valid.pointMultiply) {
      const { P, d, expected } = vector;
      deepStrictEqual(bytesToHex(tweakUtils.pointMultiply(P, d, true)), expected);
    }
  });
-should('secp256k1.pointMultiply() invalid', () => {
+  should('pointMultiply() invalid', () => {
    for (const vector of points.invalid.pointMultiply) {
      const { P, d, exception } = vector;
      throws(() => tweakUtils.pointMultiply(P, d));
    }
  });

-should('secp256k1.wychenproof vectors', () => {
+  should('wychenproof vectors', () => {
    for (let group of wp.testGroups) {
      const pubKey = secp.Point.fromHex(group.key.uncompressed);
      for (let test of group.tests) {
@@ -528,6 +542,7 @@ should('secp256k1.wychenproof vectors', () => {
      }
    }
  });
+});

 // ESM is broken.
 import url from 'url';
Author	SHA1	Message	Date
Paul Miller	b9482bb17d	Release 0.5.2.	2023-01-13 16:23:52 +01:00
Paul Miller	74475dca68	Fix lint	2023-01-13 16:02:07 +01:00
Paul Miller	f4cf21b9c8	tests: Use describe()	2023-01-13 16:00:13 +01:00
Paul Miller	5312d92b2c	edwards: Fix isTorsionFree()	2023-01-13 15:58:04 +01:00
Paul Miller	d1770c0ac7	Rename test	2023-01-13 01:29:54 +01:00
Paul Miller	2d37edf7d1	Remove utils.mod(), utils.invert()	2023-01-13 01:26:00 +01:00
Paul Miller	36998fede8	Fix sqrt	2023-01-13 01:21:51 +01:00
Paul Miller	83960d445d	Refactor: weierstrass assertValidity and others	2023-01-12 21:18:51 +01:00
Paul Miller	23cc2aa5d1	edwards, montgomery, weierstrass: refactor	2023-01-12 20:40:16 +01:00
Paul Miller	e45d7c2d25	utils: new util; ed448: small adjustment	2023-01-12 20:39:43 +01:00
Paul Miller	bfe929aac3	modular: Tonneli-Shanks refactoring	2023-01-12 20:38:42 +01:00
Paul Miller	069452dbe7	BLS, jubjub refactoring	2023-01-12 20:38:10 +01:00
Paul Miller	2e81f31d2e	ECDSA: signUnhashed(), support for key recovery from bits 2/3	2023-01-08 20:02:04 +01:00
Paul Miller	9f7df0f13b	ECDSA adjustments	2023-01-08 18:46:55 +01:00
Paul Miller	5600629bca	Refactor	2023-01-08 18:02:54 +01:00
Paul Miller	2bd5e9ac16	Release 0.5.1.	2022-12-31 10:31:10 +01:00
Paul Miller	6890c26091	Fix readme toc	2022-12-31 10:29:25 +01:00
Paul Miller	a15e3a93a9	Docs	2022-12-31 10:00:29 +01:00
Paul Miller	910c508da9	hash-to-curve: elligator in 25519, 448. Stark: adjust type	2022-12-31 07:51:29 +01:00
Paul Miller	12da04a2bb	Improve modular math	2022-12-31 07:49:42 +01:00
Paul Miller	cc2c84f040	Improve field tests	2022-12-31 07:49:09 +01:00
Paul Miller	5d42549acc	hash-to-curve: add xmd/xof support	2022-12-31 07:48:13 +01:00