Stage 3: BLS12-381 as kilic/bls12-381 computes it for drand

The field layer, an extension type over BigInt that a later implementation
with fixed limbs can replace alone; the tower Fp2, Fp6 and Fp12 with the
formulas of kilic; G1 and G2 with their compressed encodings and the
verdicts of FromCompressed (flags, the point at infinity, coordinates below
p, the curve and the subgroup, checked in G2 by psi(P) = [x]P); the optimal
ate pairing with the final exponentiation of kilic, GT serialized c1 before
c0 at every level; and the hash to G1 of RFC 9380 with the DST of Quicknet.

tool/bls12381_go_vectors.go writes test/vectors/bls12381_vectors.json with
kilic and kyber-bls12381: the frozen edge cases of datekeys-ts with their Go
verdicts recomputed, and decodings, sums, multiples, pairings, hashes, maps
and BLS signatures drawn from a fixed seed. BigInt is not constant time:
the README says where that matters.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

lib/src/bls12381_curve.dart

@@ -0,0 +1,497 @@
+/// The groups G1 and G2 of BLS12-381 and their compressed encodings (spec
+/// §12.2):
+///
+///     G1: y² = x³ + 4          over Fp
+///     G2: y² = x³ + 4·(1 + u)  over Fp2, the M-type twist
+///
+/// A point is kept in Jacobian coordinates (X, Y, Z), the affine point being
+/// (X/Z², Y/Z³) and Z = 0 the point at infinity. The decoders follow
+/// FromCompressed of kilic/bls12-381 v0.1.0, which drand/kyber-bls12381
+/// calls in UnmarshalBinary: the compression flag set; the point at infinity
+/// only as 0xc0 followed by zeros; coordinates below p; a point on the
+/// curve; and a point of the subgroup of order r. The verdicts are those of
+/// bls12381.ts of datekeys-ts, checked against Go by the vectors of the
+/// tests.
+///
+/// Not constant time: see bls12381_fp.dart.
+library;
+
+import 'dart:typed_data';
+
+import 'bls12381_fp.dart';
+import 'bls12381_tower.dart';
+
+/// The size of a compressed point of G1.
+const g1ByteLength = 48;
+
+/// The size of a compressed point of G2.
+const g2ByteLength = 96;
+
+/// A group of compressed points: G1 (48 bytes) or G2 (96 bytes).
+enum BlsGroup {
+  /// G1, over Fp: the signatures of Quicknet.
+  g1(g1ByteLength),
+
+  /// G2, over Fp2: the public key of Quicknet and U of a tlock stanza.
+  g2(g2ByteLength);
+
+  const BlsGroup(this.byteLength);
+
+  /// The size of a compressed point.
+  final int byteLength;
+}
+
+/// What the decoder of the reference makes of a compressed encoding:
+/// FromCompressed of kilic/bls12-381 and, when it succeeds, whether the
+/// point is the point at infinity.
+enum PointVerdict {
+  /// The canonical encoding of a point of the subgroup other than the point
+  /// at infinity.
+  point,
+
+  /// The canonical encoding of the point at infinity: 0xc0 and zeros.
+  identity,
+
+  /// Anything else (spec §12.2).
+  invalid,
+}
+
+/// The verdict of the decoder of the reference on [bytes] as a compressed
+/// point of [group]: [PointVerdict.point] for the canonical encoding of a
+/// point of the subgroup of order r other than the point at infinity,
+/// [PointVerdict.identity] for the encoding of the point at infinity, and
+/// [PointVerdict.invalid] for anything else, any other length included, as
+/// checkCompressedPoint of datekeys-ts.
+PointVerdict checkCompressedPoint(BlsGroup group, List<int> bytes) {
+  final Object? p = switch (group) {
+    BlsGroup.g1 => G1Point.decode(bytes),
+    BlsGroup.g2 => G2Point.decode(bytes),
+  };
+  return switch (p) {
+    null => PointVerdict.invalid,
+    G1Point(isInfinity: true) ||
+    G2Point(isInfinity: true) => PointVerdict.identity,
+    _ => PointVerdict.point,
+  };
+}
+
+// The flags of the first byte of a compressed point (spec §12.2).
+const _compressionFlag = 0x80;
+const _infinityFlag = 0x40;
+const _signFlag = 0x20;
+
+// The binary digits of r, most significant first, for the subgroup check of
+// G1.
+final String _orderBits = groupOrder.toRadixString(2);
+
+// The binary digits of |x| = 0xd201000000010000, for the subgroup check of
+// G2.
+final String _absXBits = BigInt.parse(
+  'd201000000010000',
+  radix: 16,
+).toRadixString(2);
+
+/// ψx = 1/ξ^((p − 1)/3) and ψy = 1/ξ^((p − 1)/2) of [G2Point.psi], psix and
+/// psiy of kilic. A test computes them again from their definition.
+final Fp2 psiX = Fp2.hex(
+  '0',
+  '1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4897d29650fb85f9b409427eb'
+      '4f49fffd8bfd00000000aaad',
+);
+
+/// See [psiX].
+final Fp2 psiY = Fp2.hex(
+  '135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf'
+      '3e67fa0af1ee7b04121bdea2',
+  '06af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f'
+      '72ec05f4c81084fbede3cc09',
+);
+
+// Whether [bytes] are 0xc0 and zeros, the point at infinity.
+bool _isInfinityEncoding(List<int> bytes) {
+  if (bytes[0] != _compressionFlag | _infinityFlag) return false;
+  for (var i = 1; i < bytes.length; i++) {
+    if (bytes[i] != 0) return false;
+  }
+  return true;
+}
+
+/// A point of E(Fp): y² = x³ + 4, in Jacobian coordinates. The points that
+/// this library builds are in G1, the subgroup of order r, except where a
+/// test builds others.
+final class G1Point {
+  /// The point of the Jacobian coordinates [x], [y], [z].
+  const G1Point.jacobian(this.x, this.y, this.z);
+
+  /// The affine point ([x], [y]), which must be on the curve.
+  G1Point.affine(Fp x, Fp y) : this.jacobian(x, y, Fp.one);
+
+  /// The point at infinity, the identity of the group.
+  static final G1Point infinity = G1Point.jacobian(Fp.one, Fp.one, Fp.zero);
+
+  /// The generator of G1 (spec §63, «Serialización de GT en H2», and
+  /// testdata/vectors/tlock_ibe.json).
+  static final G1Point generator = G1Point.affine(
+    Fp.hex(
+      '17f1d3a73197d7942695638c4fa9ac0fc3688c4f9774b905a14e3a3f171bac586c55e8'
+      '3ff97a1aeffb3af00adb22c6bb',
+    ),
+    Fp.hex(
+      '08b3f481e3aaa0f1a09e30ed741d8ae4fcf5e095d5d00af600db18cb2c04b3edd03cc7'
+      '44a2888ae40caa232946c5e7e1',
+    ),
+  );
+
+  /// The constant b = 4 of the curve.
+  static final Fp b = Fp(BigInt.from(4));
+
+  /// X.
+  final Fp x;
+
+  /// Y.
+  final Fp y;
+
+  /// Z, 0 for the point at infinity.
+  final Fp z;
+
+  /// Whether this is the point at infinity.
+  bool get isInfinity => z.isZero;
+
+  /// Whether the affine point ([x], [y]) is on the curve.
+  static bool isOnCurve(Fp x, Fp y) => y.square().equals(x.square() * x + b);
+
+  /// The affine coordinates, or null for the point at infinity.
+  (Fp, Fp)? toAffine() {
+    if (isInfinity) return null;
+    if (z.isOne) return (x, y);
+    final zi = z.inverse();
+    final zi2 = zi.square();
+    return (x * zi2, y * zi2 * zi);
+  }
+
+  /// Whether this and [o] are the same point.
+  bool equals(G1Point o) {
+    if (isInfinity || o.isInfinity) return isInfinity && o.isInfinity;
+    final z1z1 = z.square();
+    final z2z2 = o.z.square();
+    return (x * z2z2).equals(o.x * z1z1) &&
+        (y * z2z2 * o.z).equals(o.y * z1z1 * z);
+  }
+
+  /// The negation.
+  G1Point operator -() => G1Point.jacobian(x, -y, z);
+
+  /// Twice this: dbl-2009-l, for a = 0.
+  G1Point double() {
+    if (isInfinity) return this;
+    final a = x.square();
+    final yy = y.square();
+    final c = yy.square();
+    final d = ((x + yy).square() - a - c).double();
+    final e = a.double() + a;
+    final f = e.square();
+    final x3 = f - d.double();
+    final c8 = c.double().double().double();
+    return G1Point.jacobian(x3, e * (d - x3) - c8, (y * z).double());
+  }
+
+  /// The sum: add-2007-bl, with its special cases.
+  G1Point operator +(G1Point o) {
+    if (isInfinity) return o;
+    if (o.isInfinity) return this;
+    final z1z1 = z.square();
+    final z2z2 = o.z.square();
+    final u1 = x * z2z2;
+    final u2 = o.x * z1z1;
+    final s1 = y * o.z * z2z2;
+    final s2 = o.y * z * z1z1;
+    final h = u2 - u1;
+    final r = (s2 - s1).double();
+    if (h.isZero) return r.isZero ? double() : infinity;
+    final i = h.double().square();
+    final j = h * i;
+    final v = u1 * i;
+    final x3 = r.square() - j - v.double();
+    final y3 = r * (v - x3) - (s1 * j).double();
+    final z3 = ((z + o.z).square() - z1z1 - z2z2) * h;
+    return G1Point.jacobian(x3, y3, z3);
+  }
+
+  /// The sum with the affine point ([ax], [ay]): madd-2007-bl, with its
+  /// special cases.
+  G1Point addAffine(Fp ax, Fp ay) {
+    if (isInfinity) return G1Point.affine(ax, ay);
+    final z1z1 = z.square();
+    final u2 = ax * z1z1;
+    final s2 = ay * z * z1z1;
+    final h = u2 - x;
+    final r = (s2 - y).double();
+    if (h.isZero) return r.isZero ? double() : infinity;
+    final hh = h.square();
+    final i = hh.double().double();
+    final j = h * i;
+    final v = x * i;
+    final x3 = r.square() - j - v.double();
+    final y3 = r * (v - x3) - (y * j).double();
+    final z3 = (z + h).square() - z1z1 - hh;
+    return G1Point.jacobian(x3, y3, z3);
+  }
+
+  /// [k]·this, k >= 0, by double and add over the binary digits of [k].
+  G1Point multiply(BigInt k) {
+    if (k.isNegative) throw ArgumentError.value(k, 'k', 'negative');
+    return _multiplyBits(k.toRadixString(2));
+  }
+
+  /// Whether this point is in G1, the subgroup of order r: [r]·this is the
+  /// point at infinity.
+  bool get isInSubgroup => _multiplyBits(_orderBits).isInfinity;
+
+  G1Point _multiplyBits(String bits) {
+    final a = toAffine();
+    if (a == null) return infinity;
+    final (ax, ay) = a;
+    var acc = infinity;
+    for (var i = 0; i < bits.length; i++) {
+      acc = acc.double();
+      if (bits.codeUnitAt(i) == 0x31) acc = acc.addAffine(ax, ay);
+    }
+    return acc;
+  }
+
+  /// The 48 bytes of the compressed encoding (spec §12.2), as ToCompressed
+  /// of kilic.
+  Uint8List toBytes() {
+    final a = toAffine();
+    if (a == null) {
+      return Uint8List(g1ByteLength)..[0] = _compressionFlag | _infinityFlag;
+    }
+    final (ax, ay) = a;
+    final out = ax.toBytes();
+    out[0] |= _compressionFlag;
+    if (ay.isLexicographicallyLargest) out[0] |= _signFlag;
+    return out;
+  }
+
+  /// The point of the compressed encoding [bytes], [infinity] for the
+  /// encoding of the point at infinity, or null when FromCompressed of
+  /// kilic/bls12-381 rejects [bytes]: any length other than 48, the
+  /// compression flag clear, the infinity flag with any other bit set, x
+  /// not below p, x of no point of the curve, or a point outside G1.
+  static G1Point? decode(List<int> bytes) {
+    if (bytes.length != g1ByteLength) return null;
+    final flags = bytes[0];
+    if (flags & _compressionFlag == 0) return null;
+    if (flags & _infinityFlag != 0) {
+      return _isInfinityEncoding(bytes) ? infinity : null;
+    }
+    final raw = Uint8List.fromList(bytes);
+    raw[0] &= 0x1f;
+    final x = Fp.fromBytes(raw);
+    if (x == null) return null;
+    var y = (x.square() * x + b).sqrt();
+    if (y == null) return null;
+    if (y.isLexicographicallyLargest != (flags & _signFlag != 0)) y = -y;
+    final p = G1Point.affine(x, y);
+    return p.isInSubgroup ? p : null;
+  }
+}
+
+/// A point of E'(Fp2): y² = x³ + 4·(1 + u), in Jacobian coordinates. The
+/// points that this library builds are in G2, the subgroup of order r,
+/// except where a test builds others.
+final class G2Point {
+  /// The point of the Jacobian coordinates [x], [y], [z].
+  const G2Point.jacobian(this.x, this.y, this.z);
+
+  /// The affine point ([x], [y]), which must be on the curve.
+  G2Point.affine(Fp2 x, Fp2 y) : this.jacobian(x, y, Fp2.one);
+
+  /// The point at infinity, the identity of the group.
+  static final G2Point infinity = G2Point.jacobian(Fp2.one, Fp2.one, Fp2.zero);
+
+  /// The generator of G2 (spec §63, «Serialización de GT en H2», and
+  /// testdata/vectors/tlock_ibe.json).
+  static final G2Point generator = G2Point.affine(
+    Fp2.hex(
+      '024aa2b2f08f0a91260805272dc51051c6e47ad4fa403b02b4510b647ae3d1770bac03'
+          '26a805bbefd48056c8c121bdb8',
+      '13e02b6052719f607dacd3a088274f65596bd0d09920b61ab5da61bbdc7f5049334cf1'
+          '1213945d57e5ac7d055d042b7e',
+    ),
+    Fp2.hex(
+      '0ce5d527727d6e118cc9cdc6da2e351aadfd9baa8cbdd3a76d429a695160d12c923ac9'
+          'cc3baca289e193548608b82801',
+      '0606c4a02ea734cc32acd2b02bc28b99cb3e287e85a763af267492ab572e99ab3f370d'
+          '275cec1da1aaa9075ff05f79be',
+    ),
+  );
+
+  /// The constant b' = 4·(1 + u) of the twist.
+  static final Fp2 b = Fp2(Fp(BigInt.from(4)), Fp(BigInt.from(4)));
+
+  /// X.
+  final Fp2 x;
+
+  /// Y.
+  final Fp2 y;
+
+  /// Z, 0 for the point at infinity.
+  final Fp2 z;
+
+  /// Whether this is the point at infinity.
+  bool get isInfinity => z.isZero;
+
+  /// Whether the affine point ([x], [y]) is on the curve.
+  static bool isOnCurve(Fp2 x, Fp2 y) => y.square().equals(x.square() * x + b);
+
+  /// The affine coordinates, or null for the point at infinity.
+  (Fp2, Fp2)? toAffine() {
+    if (isInfinity) return null;
+    if (z.isOne) return (x, y);
+    final zi = z.inverse();
+    final zi2 = zi.square();
+    return (x * zi2, y * zi2 * zi);
+  }
+
+  /// Whether this and [o] are the same point.
+  bool equals(G2Point o) {
+    if (isInfinity || o.isInfinity) return isInfinity && o.isInfinity;
+    final z1z1 = z.square();
+    final z2z2 = o.z.square();
+    return (x * z2z2).equals(o.x * z1z1) &&
+        (y * z2z2 * o.z).equals(o.y * z1z1 * z);
+  }
+
+  /// The negation.
+  G2Point operator -() => G2Point.jacobian(x, -y, z);
+
+  /// Twice this: dbl-2009-l, for a = 0.
+  G2Point double() {
+    if (isInfinity) return this;
+    final a = x.square();
+    final yy = y.square();
+    final c = yy.square();
+    final d = ((x + yy).square() - a - c).double();
+    final e = a.double() + a;
+    final f = e.square();
+    final x3 = f - d.double();
+    final c8 = c.double().double().double();
+    return G2Point.jacobian(x3, e * (d - x3) - c8, (y * z).double());
+  }
+
+  /// The sum: add-2007-bl, with its special cases.
+  G2Point operator +(G2Point o) {
+    if (isInfinity) return o;
+    if (o.isInfinity) return this;
+    final z1z1 = z.square();
+    final z2z2 = o.z.square();
+    final u1 = x * z2z2;
+    final u2 = o.x * z1z1;
+    final s1 = y * o.z * z2z2;
+    final s2 = o.y * z * z1z1;
+    final h = u2 - u1;
+    final r = (s2 - s1).double();
+    if (h.isZero) return r.isZero ? double() : infinity;
+    final i = h.double().square();
+    final j = h * i;
+    final v = u1 * i;
+    final x3 = r.square() - j - v.double();
+    final y3 = r * (v - x3) - (s1 * j).double();
+    final z3 = ((z + o.z).square() - z1z1 - z2z2) * h;
+    return G2Point.jacobian(x3, y3, z3);
+  }
+
+  /// The sum with the affine point ([ax], [ay]): madd-2007-bl, with its
+  /// special cases.
+  G2Point addAffine(Fp2 ax, Fp2 ay) {
+    if (isInfinity) return G2Point.affine(ax, ay);
+    final z1z1 = z.square();
+    final u2 = ax * z1z1;
+    final s2 = ay * z * z1z1;
+    final h = u2 - x;
+    final r = (s2 - y).double();
+    if (h.isZero) return r.isZero ? double() : infinity;
+    final hh = h.square();
+    final i = hh.double().double();
+    final j = h * i;
+    final v = x * i;
+    final x3 = r.square() - j - v.double();
+    final y3 = r * (v - x3) - (y * j).double();
+    final z3 = (z + h).square() - z1z1 - hh;
+    return G2Point.jacobian(x3, y3, z3);
+  }
+
+  /// [k]·this, k >= 0, by double and add over the binary digits of [k].
+  G2Point multiply(BigInt k) {
+    if (k.isNegative) throw ArgumentError.value(k, 'k', 'negative');
+    return _multiplyBits(k.toRadixString(2));
+  }
+
+  /// ψ, the endomorphism untwist-Frobenius-twist of E': (x, y) ↦
+  /// (x̄·ψx, ȳ·ψy), as psi of kilic. It acts on G2 as the multiplication by
+  /// p, which is x modulo r.
+  G2Point psi() => G2Point.jacobian(
+    x.conjugate() * psiX,
+    y.conjugate() * psiY,
+    z.conjugate(),
+  );
+
+  /// Whether this point is in G2, the subgroup of order r: ψ(P) = [x]·P with
+  /// x = −0xd201000000010000, the test of Scott (eprint 2021/1130) that
+  /// noble uses too. It decides as [r]·P = O, which costs four times more;
+  /// the tests compare the two, and Go, on points outside the subgroup,
+  /// torsion included.
+  bool get isInSubgroup => psi().equals(-_multiplyBits(_absXBits));
+
+  G2Point _multiplyBits(String bits) {
+    final a = toAffine();
+    if (a == null) return infinity;
+    final (ax, ay) = a;
+    var acc = infinity;
+    for (var i = 0; i < bits.length; i++) {
+      acc = acc.double();
+      if (bits.codeUnitAt(i) == 0x31) acc = acc.addAffine(ax, ay);
+    }
+    return acc;
+  }
+
+  /// The 96 bytes of the compressed encoding, x.c1 || x.c0 with the flags
+  /// in the first byte (spec §12.2), as ToCompressed of kilic.
+  Uint8List toBytes() {
+    final a = toAffine();
+    if (a == null) {
+      return Uint8List(g2ByteLength)..[0] = _compressionFlag | _infinityFlag;
+    }
+    final (ax, ay) = a;
+    final out = ax.toBytes();
+    out[0] |= _compressionFlag;
+    if (ay.isLexicographicallyLargest) out[0] |= _signFlag;
+    return out;
+  }
+
+  /// The point of the compressed encoding [bytes], [infinity] for the
+  /// encoding of the point at infinity, or null when FromCompressed of
+  /// kilic/bls12-381 rejects [bytes]: any length other than 96, the
+  /// compression flag clear, the infinity flag with any other bit set, c1 or
+  /// c0 not below p, x of no point of the curve, or a point outside G2.
+  static G2Point? decode(List<int> bytes) {
+    if (bytes.length != g2ByteLength) return null;
+    final flags = bytes[0];
+    if (flags & _compressionFlag == 0) return null;
+    if (flags & _infinityFlag != 0) {
+      return _isInfinityEncoding(bytes) ? infinity : null;
+    }
+    final raw = Uint8List.fromList(bytes);
+    raw[0] &= 0x1f;
+    final c1 = Fp.fromBytes(Uint8List.sublistView(raw, 0, fpByteLength));
+    final c0 = Fp.fromBytes(Uint8List.sublistView(raw, fpByteLength));
+    if (c1 == null || c0 == null) return null;
+    final x = Fp2(c0, c1);
+    var y = (x.square() * x + b).sqrt();
+    if (y == null) return null;
+    if (y.isLexicographicallyLargest != (flags & _signFlag != 0)) y = -y;
+    final p = G2Point.affine(x, y);
+    return p.isInSubgroup ? p : null;
+  }
+}

lib/src/bls12381_fp.dart

@@ -0,0 +1,169 @@
+/// The base field Fp of BLS12-381 and its scalar field order r: the one
+/// layer of the curve code that touches the representation of a field
+/// element.
+///
+/// An element is an extension type over a [BigInt] in 0..p-1, so that it
+/// costs no wrapper at run time and the rest of the code (the tower of
+/// bls12381_tower.dart, the curves, the pairing and the hash to curve) uses
+/// only the operations declared here. A later implementation with fixed
+/// limbs replaces this file alone.
+///
+/// [BigInt] is not constant time: its operations take time and allocate
+/// memory that depend on the values. That is acceptable where the values are
+/// public, as in the verification of a release, whose signature is public
+/// once drand publishes it, and in the decryption of a tlock stanza, whose
+/// inputs are that signature and the stanza. It is not where they are
+/// secret: in the encryption of a tlock stanza, sigma and r leak through the
+/// timing of this code (see ibe.dart).
+///
+/// The [BigInt] of Dart is exact on the VM and compiled to JavaScript alike.
+library;
+
+import 'dart:typed_data';
+
+import 'bytes.dart';
+
+/// The modulus p of the base field (spec §12.2).
+final BigInt fpModulus = BigInt.parse(
+  '1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb1'
+  '53ffffb9feffffffffaaab',
+  radix: 16,
+);
+
+/// The prime order r of G1, G2 and GT, the order of the scalar field (spec
+/// §12.2).
+final BigInt groupOrder = BigInt.parse(
+  '73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001',
+  radix: 16,
+);
+
+final BigInt _p = fpModulus;
+// (p + 1) / 4: the exponent of the square root, as p = 3 mod 4.
+final BigInt _sqrtExponent = (_p + BigInt.one) >> 2;
+// (p - 1) / 2: the largest y that is not lexicographically larger than -y.
+final BigInt _halfP = (_p - BigInt.one) >> 1;
+
+/// The size of an element of Fp in bytes, big-endian.
+const fpByteLength = 48;
+
+/// An element of Fp, the integers modulo p.
+extension type const Fp._(BigInt _v) {
+  /// The element [v], which must be in 0..p-1.
+  factory Fp(BigInt v) {
+    if (v.isNegative || v >= _p) {
+      throw ArgumentError.value(v, 'v', 'not in 0..p-1');
+    }
+    return Fp._(v);
+  }
+
+  /// The element [v] modulo p, for any integer [v].
+  factory Fp.reduce(BigInt v) => Fp._(v % _p);
+
+  /// The element of the hexadecimal [hex], which must be below p: the
+  /// constants of the curves.
+  factory Fp.hex(String hex) => Fp(BigInt.parse(hex, radix: 16));
+
+  /// The element of the 48 big-endian bytes [b], or null when they encode p
+  /// or more: the coordinate decoders of spec §12.2 reject it rather than
+  /// reduce it.
+  static Fp? fromBytes(List<int> b) {
+    if (b.length != fpByteLength) {
+      throw ArgumentError.value(b.length, 'b', 'want $fpByteLength bytes');
+    }
+    final v = _bigFromBytes(b);
+    return v < _p ? Fp._(v) : null;
+  }
+
+  /// The element of any number of big-endian bytes, reduced modulo p: the
+  /// 64 bytes of hash_to_field of RFC 9380.
+  static Fp fromBytesReduced(List<int> b) => Fp._(_bigFromBytes(b) % _p);
+
+  /// 0.
+  static final Fp zero = Fp._(BigInt.zero);
+
+  /// 1.
+  static final Fp one = Fp._(BigInt.one);
+
+  /// 1/2, which is (p + 1)/2.
+  static final Fp half = Fp._((_p + BigInt.one) >> 1);
+
+  /// The integer in 0..p-1 of this element.
+  BigInt toBigInt() => _v;
+
+  /// The 48 big-endian bytes of this element.
+  Uint8List toBytes() =>
+      fromHex(_v.toRadixString(16).padLeft(2 * fpByteLength, '0'));
+
+  /// Whether this is 0.
+  bool get isZero => _v.sign == 0;
+
+  /// Whether this is 1.
+  bool get isOne => _v == BigInt.one;
+
+  /// Whether this and [other] are the same element.
+  bool equals(Fp other) => _v == other._v;
+
+  /// The parity of the integer of this element: sgn0 of RFC 9380 for m = 1.
+  bool get isOdd => _v.isOdd;
+
+  /// Whether this element is larger than its negation as an integer in
+  /// 0..p-1, which sets the sign bit of a compressed point (spec §12.2).
+  bool get isLexicographicallyLargest => _v > _halfP;
+
+  /// The sum.
+  Fp operator +(Fp o) {
+    final s = _v + o._v;
+    return Fp._(s >= _p ? s - _p : s);
+  }
+
+  /// The difference.
+  Fp operator -(Fp o) {
+    final d = _v - o._v;
+    return Fp._(d.isNegative ? d + _p : d);
+  }
+
+  /// The negation.
+  Fp operator -() => _v.sign == 0 ? this : Fp._(_p - _v);
+
+  /// The product.
+  Fp operator *(Fp o) => Fp._((_v * o._v) % _p);
+
+  /// The square.
+  Fp square() => Fp._((_v * _v) % _p);
+
+  /// Twice this.
+  Fp double() => this + this;
+
+  /// The inverse, or 0 for 0, as the inverse of kilic/bls12-381.
+  Fp inverse() => _v.sign == 0 ? this : Fp._(_v.modInverse(_p));
+
+  /// This to the power [e], e >= 0.
+  Fp pow(BigInt e) => Fp._(_v.modPow(e, _p));
+
+  /// A square root, or null when this is not a square. Of the two roots
+  /// it returns this^((p + 1) / 4); when this is not a square, that power is
+  /// a root of its negation instead, which [sqrtOrNegatedRoot] exposes.
+  Fp? sqrt() {
+    final (root, isSquare) = sqrtOrNegatedRoot();
+    return isSquare ? root : null;
+  }
+
+  /// this^((p + 1) / 4) and whether it is a square root of this. When it is
+  /// not, it is a square root of the negation of this, since -1 is not a
+  /// square in Fp: the one exponentiation that the map of RFC 9380 needs.
+  (Fp, bool) sqrtOrNegatedRoot() {
+    final s = _v.modPow(_sqrtExponent, _p);
+    return (Fp._(s), (s * s) % _p == _v);
+  }
+
+  /// (a * b - c * d) mod p, reduced once.
+  static Fp mulSub(Fp a, Fp b, Fp c, Fp d) =>
+      Fp._((a._v * b._v - c._v * d._v) % _p);
+
+  /// (a * b + c * d) mod p, reduced once.
+  static Fp mulAdd(Fp a, Fp b, Fp c, Fp d) =>
+      Fp._((a._v * b._v + c._v * d._v) % _p);
+}
+
+BigInt _bigFromBytes(List<int> b) =>
+    b.isEmpty ? BigInt.zero : BigInt.parse(toHex(b), radix: 16);

lib/src/bls12381_hash.dart

@@ -0,0 +1,340 @@
+/// Hashing to G1 of BLS12-381 per RFC 9380, suite
+/// BLS12381G1_XMD:SHA-256_SSWU_RO_, as HashToCurve of kilic/bls12-381
+/// v0.1.0, which drand/kyber-bls12381 calls: expand_message_xmd with
+/// SHA-256, hash_to_field to two elements of Fp, the simplified SWU map to
+/// the curve E' 11-isogenous to E, the 11-isogeny and the clearing of the
+/// cofactor by h_eff = 0xd201000000010001.
+///
+/// Quicknet (bls-unchained-g1-rfc9380) signs the round on G1 with the DST
+/// [quicknetDst], and tlock hashes the identity of a round to G1 with the
+/// same DST (drand/kyber-bls12381.NewBLS12381Suite). G2 is never hashed to:
+/// no use of the protocol needs it.
+///
+/// As kilic, the two outputs of the map are added on E' before the isogeny,
+/// which is a homomorphism: the result is the RFC's. Where kilic departs from
+/// the RFC, it is in cases that no input reaches: g(x1) = 0, which kilic
+/// takes for a non-square and the RFC for a square, never happens, E'
+/// having no point of order 2; and two outputs of the map that are the same
+/// point or opposite points, with probability about 2^-380 for a hash, which
+/// this code adds on E' as the group law of E' has it, where kilic doubles
+/// with the formula of E. The exceptional u of the map (u = 0 or Z·u² = −1)
+/// give x1 = B/(Z·A) in both, the zInv of kilic being −1/Z, and Go's
+/// vectors cover them.
+///
+/// Not constant time: see bls12381_fp.dart. The inputs it hashes, the
+/// identities of rounds, are public.
+library;
+
+import 'dart:typed_data';
+
+import 'package:crypto/crypto.dart' as crypto;
+
+import 'bls12381_curve.dart';
+import 'bls12381_fp.dart';
+import 'bytes.dart';
+
+/// The DST of RFC 9380 for G1 that Quicknet and tlock use:
+/// BLS_SIG_BLS12381G1_XMD:SHA-256_SSWU_RO_NUL_.
+const quicknetDst = 'BLS_SIG_BLS12381G1_XMD:SHA-256_SSWU_RO_NUL_';
+
+/// expand_message_xmd of RFC 9380, section 5.3.1, with SHA-256: [length]
+/// bytes from [msg] and the domain separation tag [dst]. Throws an
+/// [ArgumentError] for a [dst] longer than 255 bytes or a [length] that
+/// needs more than 255 blocks of SHA-256, as the RFC requires.
+Uint8List expandMessageXmd(List<int> msg, List<int> dst, int length) {
+  const blockBytes = 32;
+  final ell = (length + blockBytes - 1) ~/ blockBytes;
+  if (dst.length > 255) {
+    throw ArgumentError.value(dst.length, 'dst', 'longer than 255 bytes');
+  }
+  if (length < 0 || ell > 255 || length > 65535) {
+    throw ArgumentError.value(length, 'length', 'out of range');
+  }
+  final dstPrime = concatBytes([
+    dst,
+    [dst.length],
+  ]);
+  final b0 = crypto.sha256
+      .convert(
+        concatBytes([
+          Uint8List(64),
+          msg,
+          [length >> 8, length & 0xff, 0],
+          dstPrime,
+        ]),
+      )
+      .bytes;
+  // b_1 is computed even for an empty output, as the RFC does.
+  final out = Uint8List((ell == 0 ? 1 : ell) * blockBytes);
+  var previous = crypto.sha256
+      .convert(
+        concatBytes([
+          b0,
+          [1],
+          dstPrime,
+        ]),
+      )
+      .bytes;
+  out.setRange(0, blockBytes, previous);
+  for (var i = 2; i <= ell; i++) {
+    final mixed = Uint8List(blockBytes);
+    for (var j = 0; j < blockBytes; j++) {
+      mixed[j] = b0[j] ^ previous[j];
+    }
+    previous = crypto.sha256
+        .convert(
+          concatBytes([
+            mixed,
+            [i],
+            dstPrime,
+          ]),
+        )
+        .bytes;
+    out.setRange((i - 1) * blockBytes, i * blockBytes, previous);
+  }
+  return Uint8List.sublistView(out, 0, length);
+}
+
+/// hash_to_field of RFC 9380, section 5.2, for Fp with L = 64: [count]
+/// elements from [msg] and [dst].
+List<Fp> hashToField(List<int> msg, List<int> dst, int count) {
+  final uniform = expandMessageXmd(msg, dst, count * 64);
+  return [
+    for (var i = 0; i < count; i++)
+      Fp.fromBytesReduced(Uint8List.sublistView(uniform, i * 64, i * 64 + 64)),
+  ];
+}
+
+// The curve E': y² = x³ + A'·x + B', 11-isogenous to E, and Z = 11 of the
+// simplified SWU map (RFC 9380, section 8.8.1).
+final Fp _a = Fp.hex(
+  '00144698a3b8e9433d693a02c96d4982b0ea985383ee66a8d8e8981aefd881ac98936f8d'
+  'a0e0f97f5cf428082d584c1d',
+);
+final Fp _b = Fp.hex(
+  '12e2908d11688030018b12e8753eee3b2016c1f0f24f4070a0b9c14fcef35ef55a23215a'
+  '316ceaa5d1cc48e98e172be0',
+);
+final Fp _z = Fp(BigInt.from(11));
+final Fp _minusBOverA = -(_b * _a.inverse());
+final Fp _exceptionalX1 = _b * (_z * _a).inverse();
+// √(−Z): −Z is a square, Z and −1 being none.
+final Fp _sqrtMinusZ = (-_z).sqrt()!;
+
+/// The simplified SWU map of [u] to the curve E', as swuMapG1 of kilic,
+/// with the one exponentiation of RFC 9380, appendix F.2.1.2: when g(x1) is
+/// not a square, the power that would be its root is a root of −g(x1), from
+/// which the root of g(x2) follows. The sign of y is that of [u] (sgn0).
+(Fp, Fp) mapToIsogenousCurve(Fp u) {
+  final tv0 = _z * u.square();
+  final tv1 = tv0.square();
+  final den = (tv0 + tv1).inverse();
+  // The exceptional case of RFC 9380, section 6.6.2: x1 = B/(Z·A).
+  final x1 = den.isZero ? _exceptionalX1 : (den + Fp.one) * _minusBOverA;
+  final gx1 = (x1.square() + _a) * x1 + _b;
+  final (s, isSquare) = gx1.sqrtOrNegatedRoot();
+  final Fp x;
+  var y = s;
+  if (isSquare) {
+    x = x1;
+  } else {
+    // g(x2) = g(x1)·Z³·u⁶ and s² = −g(x1): √g(x2) = s·u³·Z·√(−Z).
+    x = tv0 * x1;
+    y = s * u.square() * u * _z * _sqrtMinusZ;
+  }
+  if (y.isOdd != u.isOdd) y = -y;
+  return (x, y);
+}
+
+/// The 11-isogeny from E' to E of RFC 9380, appendix E.2, as isogenyMapG1
+/// of kilic, with its constants.
+(Fp, Fp) isogenyMap(Fp x, Fp y) {
+  Fp eval(List<Fp> k) {
+    var acc = k.last;
+    for (var i = k.length - 2; i >= 0; i--) {
+      acc = acc * x + k[i];
+    }
+    return acc;
+  }
+
+  final c = _isogenyConstants;
+  return (
+    eval(c[0]) * eval(c[1]).inverse(),
+    y * eval(c[2]) * eval(c[3]).inverse(),
+  );
+}
+
+// The sum of two affine points of E', null for the point at infinity.
+(Fp, Fp)? _addOnIsogenousCurve((Fp, Fp) p, (Fp, Fp) q) {
+  final (x0, y0) = p;
+  final (x1, y1) = q;
+  final Fp lambda;
+  if (!x0.equals(x1)) {
+    lambda = (y1 - y0) * (x1 - x0).inverse();
+  } else if (y0.equals(y1) && !y0.isZero) {
+    final xx = x0.square();
+    lambda = (xx.double() + xx + _a) * y0.double().inverse();
+  } else {
+    return null;
+  }
+  final x2 = lambda.square() - x0 - x1;
+  return (x2, lambda * (x0 - x2) - y0);
+}
+
+/// h_eff of RFC 9380 for G1, 1 − x = 0xd201000000010001.
+final BigInt g1CofactorEff = BigInt.parse('d201000000010001', radix: 16);
+
+/// hash_to_curve of RFC 9380 to G1 with the suite
+/// BLS12381G1_XMD:SHA-256_SSWU_RO_ and the domain separation tag [dst], as
+/// HashToCurve of kilic.
+G1Point hashToG1(List<int> msg, String dst) {
+  final [u0, u1] = hashToField(msg, utf8Bytes(dst), 2);
+  final q = _addOnIsogenousCurve(
+    mapToIsogenousCurve(u0),
+    mapToIsogenousCurve(u1),
+  );
+  if (q == null) return G1Point.infinity;
+  final (x, y) = isogenyMap(q.$1, q.$2);
+  return G1Point.affine(x, y).multiply(g1CofactorEff);
+}
+
+/// The map of one element [u] to G1, as MapToCurve of kilic: the simplified
+/// SWU map, the isogeny and the clearing of the cofactor. Only the tests use
+/// it, to reach the exceptional values of [u] that no hash reaches.
+G1Point mapToG1(Fp u) {
+  final (x0, y0) = mapToIsogenousCurve(u);
+  final (x, y) = isogenyMap(x0, y0);
+  return G1Point.affine(x, y).multiply(g1CofactorEff);
+}
+
+/// The coefficients of the 11-isogeny, lowest degree first: x_num, x_den,
+/// y_num and y_den, the denominators monic (isogenyConstansG1 of kilic,
+/// out of its Montgomery form).
+List<List<Fp>> get isogenyConstants => _isogenyConstants;
+
+final List<List<Fp>> _isogenyConstants = [
+  for (final k in [_isoXNum, _isoXDen, _isoYNum, _isoYDen])
+    [for (final h in k) Fp.hex(h)],
+];
+
+const _isoXNum = [
+  '11a05f2b1e833340b809101dd99815856b303e88a2d7005ff2627b56cdb4e2c85610c2d5'
+      'f2e62d6eaeac1662734649b7',
+  '17294ed3e943ab2f0588bab22147a81c7c17e75b2f6a8417f565e33c70d1e86b4838f2a6'
+      'f318c356e834eef1b3cb83bb',
+  '0d54005db97678ec1d1048c5d10a9a1bce032473295983e56878e501ec68e25c958c3e3d'
+      '2a09729fe0179f9dac9edcb0',
+  '1778e7166fcc6db74e0609d307e55412d7f5e4656a8dbf25f1b33289f1b330835336e25c'
+      'e3107193c5b388641d9b6861',
+  '0e99726a3199f4436642b4b3e4118e5499db995a1257fb3f086eeb65982fac18985a286f'
+      '301e77c451154ce9ac8895d9',
+  '1630c3250d7313ff01d1201bf7a74ab5db3cb17dd952799b9ed3ab9097e68f90a0870d2d'
+      'cae73d19cd13c1c66f652983',
+  '0d6ed6553fe44d296a3726c38ae652bfb11586264f0f8ce19008e218f9c86b2a8da25128'
+      'c1052ecaddd7f225a139ed84',
+  '17b81e7701abdbe2e8743884d1117e53356de5ab275b4db1a682c62ef0f2753339b7c8f8'
+      'c8f475af9ccb5618e3f0c88e',
+  '080d3cf1f9a78fc47b90b33563be990dc43b756ce79f5574a2c596c928c5d1de4fa295f2'
+      '96b74e956d71986a8497e317',
+  '169b1f8e1bcfa7c42e0c37515d138f22dd2ecb803a0c5c99676314baf4bb1b7fa3190b2e'
+      'dc0327797f241067be390c9e',
+  '10321da079ce07e272d8ec09d2565b0dfa7dccdde6787f96d50af36003b14866f69b771f'
+      '8c285decca67df3f1605fb7b',
+  '06e08c248e260e70bd1e962381edee3d31d79d7e22c837bc23c0bf1bc24c6b68c24b1b80'
+      'b64d391fa9c8ba2e8ba2d229',
+];
+
+const _isoXDen = [
+  '08ca8d548cff19ae18b2e62f4bd3fa6f01d5ef4ba35b48ba9c9588617fc8ac62b558d681'
+      'be343df8993cf9fa40d21b1c',
+  '12561a5deb559c4348b4711298e536367041e8ca0cf0800c0126c2588c48bf5713daa884'
+      '6cb026e9e5c8276ec82b3bff',
+  '0b2962fe57a3225e8137e629bff2991f6f89416f5a718cd1fca64e00b11aceacd6a3d096'
+      '7c94fedcfcc239ba5cb83e19',
+  '03425581a58ae2fec83aafef7c40eb545b08243f16b1655154cca8abc28d6fd04976d524'
+      '3eecf5c4130de8938dc62cd8',
+  '13a8e162022914a80a6f1d5f43e7a07dffdfc759a12062bb8d6b44e833b306da9bd29ba8'
+      '1f35781d539d395b3532a21e',
+  '0e7355f8e4e667b955390f7f0506c6e9395735e9ce9cad4d0a43bcef24b8982f7400d24b'
+      'c4228f11c02df9a29f6304a5',
+  '0772caacf16936190f3e0c63e0596721570f5799af53a1894e2e073062aede9cea73b353'
+      '8f0de06cec2574496ee84a3a',
+  '14a7ac2a9d64a8b230b3f5b074cf01996e7f63c21bca68a81996e1cdf9822c580fa5b948'
+      '9d11e2d311f7d99bbdcc5a5e',
+  '0a10ecf6ada54f825e920b3dafc7a3cce07f8d1d7161366b74100da67f39883503826692'
+      'abba43704776ec3a79a1d641',
+  '095fc13ab9e92ad4476d6e3eb3a56680f682b4ee96f7d03776df533978f31c1593174e4b'
+      '4b7865002d6384d168ecdd0a',
+  '000000000000000000000000000000000000000000000000000000000000000000000000'
+      '000000000000000000000001',
+];
+
+const _isoYNum = [
+  '090d97c81ba24ee0259d1f094980dcfa11ad138e48a869522b52af6c956543d3cd0c7aee'
+      '9b3ba3c2be9845719707bb33',
+  '134996a104ee5811d51036d776fb46831223e96c254f383d0f906343eb67ad34d6c56711'
+      '962fa8bfe097e75a2e41c696',
+  '00cc786baa966e66f4a384c86a3b49942552e2d658a31ce2c344be4b91400da7d26d5216'
+      '28b00523b8dfe240c72de1f6',
+  '01f86376e8981c217898751ad8746757d42aa7b90eeb791c09e4a3ec03251cf9de405aba'
+      '9ec61deca6355c77b0e5f4cb',
+  '08cc03fdefe0ff135caf4fe2a21529c4195536fbe3ce50b879833fd221351adc2ee7f8dc'
+      '099040a841b6daecf2e8fedb',
+  '16603fca40634b6a2211e11db8f0a6a074a7d0d4afadb7bd76505c3d3ad5544e203f6326'
+      'c95a807299b23ab13633a5f0',
+  '04ab0b9bcfac1bbcb2c977d027796b3ce75bb8ca2be184cb5231413c4d634f3747a87ac2'
+      '460f415ec961f8855fe9d6f2',
+  '0987c8d5333ab86fde9926bd2ca6c674170a05bfe3bdd81ffd038da6c26c842642f64550'
+      'fedfe935a15e4ca31870fb29',
+  '09fc4018bd96684be88c9e221e4da1bb8f3abd16679dc26c1e8b6e6a1f20cabe69d65201'
+      'c78607a360370e577bdba587',
+  '0e1bba7a1186bdb5223abde7ada14a23c42a0ca7915af6fe06985e7ed1e4d43b9b3f7055'
+      'dd4eba6f2bafaaebca731c30',
+  '19713e47937cd1be0dfd0b8f1d43fb93cd2fcbcb6caf493fd1183e416389e61031bf3a5c'
+      'ce3fbafce813711ad011c132',
+  '18b46a908f36f6deb918c143fed2edcc523559b8aaf0c2462e6bfe7f911f643249d9cdf4'
+      '1b44d606ce07c8a4d0074d8e',
+  '0b182cac101b9399d155096004f53f447aa7b12a3426b08ec02710e807b4633f06c851c1'
+      '919211f20d4c04f00b971ef8',
+  '0245a394ad1eca9b72fc00ae7be315dc757b3b080d4c158013e6632d3c40659cc6cf90ad'
+      '1c232a6442d9d3f5db980133',
+  '05c129645e44cf1102a159f748c4a3fc5e673d81d7e86568d9ab0f5d396a7ce46ba1049b'
+      '6579afb7866b1e715475224b',
+  '15e6be4e990f03ce4ea50b3b42df2eb5cb181d8f84965a3957add4fa95af01b2b665027e'
+      'fec01c7704b456be69c8b604',
+];
+
+const _isoYDen = [
+  '16112c4c3a9c98b252181140fad0eae9601a6de578980be6eec3232b5be72e7a07f3688e'
+      'f60c206d01479253b03663c1',
+  '1962d75c2381201e1a0cbd6c43c348b885c84ff731c4d59ca4a10356f453e01f78a42607'
+      '63529e3532f6102c2e49a03d',
+  '058df3306640da276faaae7d6e8eb15778c4855551ae7f310c35a5dd279cd2eca6757cd6'
+      '36f96f891e2538b53dbf67f2',
+  '16b7d288798e5395f20d23bf89edb4d1d115c5dbddbcd30e123da489e726af41727364f2'
+      'c28297ada8d26d98445f5416',
+  '0be0e079545f43e4b00cc912f8228ddcc6d19c9f0f69bbb0542eda0fc9dec916a20b15dc'
+      '0fd2ededda39142311a5001d',
+  '08d9e5297186db2d9fb266eaac783182b70152c65550d881c5ecd87b6f0f5a6449f38db9'
+      'dfa9cce202c6477faaf9b7ac',
+  '166007c08a99db2fc3ba8734ace9824b5eecfdfa8d0cf8ef5dd365bc400a0051d5fa9c01'
+      'a58b1fb93d1a1399126a775c',
+  '16a3ef08be3ea7ea03bcddfabba6ff6ee5a4375efa1f4fd7feb34fd206357132b920f5b0'
+      '0801dee460ee415a15812ed9',
+  '1866c8ed336c61231a1be54fd1d74cc4f9fb0ce4c6af5920abc5750c4bf39b4852cfe2f7'
+      'bb9248836b233d9d55535d4a',
+  '167a55cda70a6e1cea820597d94a84903216f763e13d87bb5308592e7ea7d4fbc7385ea3'
+      'd529b35e346ef48bb8913f55',
+  '04d2f259eea405bd48f010a01ad2911d9c6dd039bb61a6290e591b36e636a5c871a5c29f'
+      '4f83060400f8b49cba8f6aa8',
+  '0accbb67481d033ff5852c1e48c50c477f94ff8aefce42d28c0f9a88cea7913516f96898'
+      '6f7ebbea9684b529e2561092',
+  '0ad6b9514c767fe3c3613144b45f1496543346d98adf02267d5ceef9a00d9b8693000763'
+      'e3b90ac11e99b138573345cc',
+  '02660400eb2e4f3b628bdd0d53cd76f2bf565b94e72927c1cb748df27942480e420517bd'
+      '8714cc80d1fadc1326ed06f7',
+  '0e0fa1d816ddc03e6b24255e0d7819c171c40f65e273b853324efcd6356caa205ca2f570'
+      'f13497804415473a1d634b8f',
+  '000000000000000000000000000000000000000000000000000000000000000000000000'
+      '000000000000000000000001',
+];

lib/src/bls12381_pairing.dart

@@ -0,0 +1,159 @@
+/// The optimal ate pairing of BLS12-381, e: G1 × G2 → GT, as the Engine of
+/// kilic/bls12-381 v0.1.0 computes it for drand/kyber-bls12381: the Miller
+/// loop over |x| = 0xd201000000010000 with the lines of its doubling and
+/// addition steps, the conjugation for x < 0, and its final exponentiation,
+/// whose hard part raises to 3·(p⁴ − p² + 1)/r (eprint 2016/130). Another
+/// final exponent gives another, equally bilinear, value: the tlock IBE
+/// hashes this one (spec §63, «Serialización de GT en H2»).
+///
+/// Not constant time: see bls12381_fp.dart.
+library;
+
+import 'bls12381_curve.dart';
+import 'bls12381_fp.dart';
+import 'bls12381_tower.dart';
+
+/// The binary digits of |x| = 0xd201000000010000 below its top bit, most
+/// significant first: the Miller loop runs from bit 62 down to bit 0.
+final String millerLoopBits = BigInt.parse(
+  'd201000000010000',
+  radix: 16,
+).toRadixString(2).substring(1);
+
+// The coefficients of the line of a step of the Miller loop.
+typedef _Line = (Fp2, Fp2, Fp2);
+
+// The lines of the Miller loop of an affine point of G2: per bit, the line
+// of its doubling step and, for a bit set, of its addition step. The
+// precompute of kilic.
+List<_Line> _lines(Fp2 qx, Fp2 qy) {
+  final out = <_Line>[];
+  var rx = qx;
+  var ry = qy;
+  var rz = Fp2.one;
+  for (var i = 0; i < millerLoopBits.length; i++) {
+    // doublingStep of kilic.
+    final t0 = (rx * ry).mulFp(Fp.half);
+    final t1 = ry.square();
+    final t2 = rz.square();
+    final t3 = G2Point.b * (t2.double() + t2);
+    final t4 = t3.double() + t3;
+    final t5 = (t1 + t4).mulFp(Fp.half);
+    final t6 = (ry + rz).square() - (t2 + t1);
+    final x2 = rx.square();
+    out.add((t3 - t1, x2.double() + x2, -t6));
+    final t3s = t3.square();
+    rx = (t1 - t4) * t0;
+    ry = t5.square() - (t3s.double() + t3s);
+    rz = t1 * t6;
+    if (millerLoopBits.codeUnitAt(i) == 0x31) {
+      // additionStep of kilic.
+      final a0 = ry - qy * rz;
+      final a1 = rx - qx * rz;
+      final a3 = a1.square();
+      final a4 = a1 * a3;
+      final a2 = rz * a0.square();
+      final a3x = a3 * rx;
+      final a5 = a4 - a3x.double() + a2;
+      final nx = a1 * a5;
+      final ny = (a3x - a5) * a0 - ry * a4;
+      final nz = rz * a4;
+      out.add((a0 * qx - a1 * qy, -a0, a1));
+      rx = nx;
+      ry = ny;
+      rz = nz;
+    }
+  }
+  return out;
+}
+
+/// The Miller loop of the pairs ([P], [Q]) of affine points, as kilic: the
+/// product of the line functions of every pair, conjugated because x < 0.
+/// A pair with a point at infinity contributes nothing, as in kilic's
+/// AddPair.
+Fp12 millerLoop(List<(G1Point, G2Point)> pairs) {
+  final points = <(Fp, Fp, List<_Line>)>[];
+  for (final (p, q) in pairs) {
+    final pa = p.toAffine();
+    final qa = q.toAffine();
+    if (pa == null || qa == null) continue;
+    final (qx, qy) = qa;
+    points.add((pa.$1, pa.$2, _lines(qx, qy)));
+  }
+  var f = Fp12.one;
+  if (points.isEmpty) return f;
+  var j = 0;
+  for (var i = 0; i < millerLoopBits.length; i++) {
+    if (i != 0) f = f.square();
+    for (final (px, py, lines) in points) {
+      final (c0, c1, c2) = lines[j];
+      f = f.mul014(c0, c1.mulFp(px), c2.mulFp(py));
+    }
+    if (millerLoopBits.codeUnitAt(i) == 0x31) {
+      j++;
+      for (final (px, py, lines) in points) {
+        final (c0, c1, c2) = lines[j];
+        f = f.mul014(c0, c1.mulFp(px), c2.mulFp(py));
+      }
+    }
+    j++;
+  }
+  return f.conjugate();
+}
+
+// a^x of an element of the cyclotomic subgroup, x = -0xd201000000010000:
+// the addition chain of kilic (from blst), then the conjugation for x < 0.
+Fp12 _expByX(Fp12 a) {
+  var c = a.cyclotomicSquare();
+  for (final n in const [2, 3, 9, 32, 16]) {
+    c = c * a;
+    for (var i = 0; i < n; i++) {
+      c = c.cyclotomicSquare();
+    }
+  }
+  return c.conjugate();
+}
+
+/// The final exponentiation of kilic: f^((p⁶ − 1)(p² + 1)), then the hard
+/// part f^(λ0 + λ1·p + λ2·p² + λ3·p³) = f^(3·(p⁴ − p² + 1)/r) of eprint
+/// 2016/130, section 3.
+Fp12 finalExponentiation(Fp12 f) {
+  // The easy part.
+  var t2 = f.conjugate() * f.inverse();
+  t2 = t2.frobenius(2) * t2;
+  // The hard part, step by step as kilic.
+  var t1 = t2.conjugate().cyclotomicSquare();
+  var t3 = _expByX(t2);
+  var t4 = t3.cyclotomicSquare();
+  var t5 = t1 * t3;
+  t1 = _expByX(t5);
+  final t0 = _expByX(t1);
+  var t6 = _expByX(t0);
+  t6 = t6 * t4;
+  t4 = _expByX(t6);
+  t5 = t5.conjugate();
+  t4 = t4 * t5 * t2;
+  t5 = (t2.conjugate() * t6).frobenius(1);
+  t3 = (t3 * t0).frobenius(2);
+  t1 = (t1 * t2).frobenius(3);
+  return t3 * t1 * t5 * t4;
+}
+
+/// e([p], [q]), 1 when either is the point at infinity, as Pair of
+/// drand/kyber-bls12381.
+Fp12 pairing(G1Point p, G2Point q) {
+  if (p.isInfinity || q.isInfinity) return Fp12.one;
+  return finalExponentiation(millerLoop([(p, q)]));
+}
+
+/// Whether the product of e(P, Q) over [pairs] is 1: one Miller loop over
+/// all the pairs and one final exponentiation, as Check of kilic's Engine.
+/// The pairs with a point at infinity are left out, as kilic's AddPair does.
+bool pairingCheck(List<(G1Point, G2Point)> pairs) {
+  final live = [
+    for (final (p, q) in pairs)
+      if (!p.isInfinity && !q.isInfinity) (p, q),
+  ];
+  if (live.isEmpty) return true;
+  return finalExponentiation(millerLoop(live)).isOne;
+}

lib/src/bls12381_tower.dart

@@ -0,0 +1,358 @@
+/// The tower of extensions of Fp of BLS12-381, as kilic/bls12-381 builds it
+/// (spec §63, «Serialización de GT en H2»):
+///
+///     Fp2  = Fp[u]  / (u² + 1)
+///     Fp6  = Fp2[v] / (v³ − ξ),  ξ = u + 1
+///     Fp12 = Fp6[w] / (w² − v)
+///
+/// The formulas are those of kilic/bls12-381 v0.1.0, which drand/kyber
+/// pairs with: an element of GT is the value its pairing computes, and H2 of
+/// the tlock IBE hashes it in the order of [Fp12.toBytes], c1 before c0 at
+/// every level. Everything is written over the operations of the field layer
+/// (bls12381_fp.dart), and like it, is not constant time.
+library;
+
+import 'dart:typed_data';
+
+import 'bls12381_fp.dart';
+import 'bytes.dart';
+
+/// An element c0 + c1·u of Fp2.
+final class Fp2 {
+  /// The element [c0] + [c1]·u.
+  const Fp2(this.c0, this.c1);
+
+  /// The element of two hexadecimal coordinates: the constants of the
+  /// curves.
+  Fp2.hex(String c0, String c1) : this(Fp.hex(c0), Fp.hex(c1));
+
+  /// 0.
+  static final Fp2 zero = Fp2(Fp.zero, Fp.zero);
+
+  /// 1.
+  static final Fp2 one = Fp2(Fp.one, Fp.zero);
+
+  /// The coefficient of 1.
+  final Fp c0;
+
+  /// The coefficient of u.
+  final Fp c1;
+
+  /// Whether this is 0.
+  bool get isZero => c0.isZero && c1.isZero;
+
+  /// Whether this is 1.
+  bool get isOne => c0.isOne && c1.isZero;
+
+  /// Whether this and [other] are the same element.
+  bool equals(Fp2 other) => c0.equals(other.c0) && c1.equals(other.c1);
+
+  /// The sum.
+  Fp2 operator +(Fp2 o) => Fp2(c0 + o.c0, c1 + o.c1);
+
+  /// The difference.
+  Fp2 operator -(Fp2 o) => Fp2(c0 - o.c0, c1 - o.c1);
+
+  /// The negation.
+  Fp2 operator -() => Fp2(-c0, -c1);
+
+  /// The product: (a0·b0 − a1·b1) + (a0·b1 + a1·b0)·u.
+  Fp2 operator *(Fp2 o) =>
+      Fp2(Fp.mulSub(c0, o.c0, c1, o.c1), Fp.mulAdd(c0, o.c1, c1, o.c0));
+
+  /// The square: (c0 + c1)(c0 − c1) + 2·c0·c1·u.
+  Fp2 square() => Fp2((c0 + c1) * (c0 - c1), c0.double() * c1);
+
+  /// Twice this.
+  Fp2 double() => Fp2(c0.double(), c1.double());
+
+  /// The product by an element of Fp.
+  Fp2 mulFp(Fp k) => Fp2(c0 * k, c1 * k);
+
+  /// The product by ξ = u + 1, the non-residue of Fp6.
+  Fp2 mulByNonResidue() => Fp2(c0 - c1, c0 + c1);
+
+  /// The conjugate c0 − c1·u, which is also this^p.
+  Fp2 conjugate() => Fp2(c0, -c1);
+
+  /// The inverse, 0 for 0.
+  Fp2 inverse() {
+    final t = Fp.mulAdd(c0, c0, c1, c1).inverse();
+    return Fp2(c0 * t, -(c1 * t));
+  }
+
+  /// This to the power [e], e >= 0.
+  Fp2 pow(BigInt e) {
+    var z = Fp2.one;
+    final bits = e.toRadixString(2);
+    for (var i = 0; i < bits.length; i++) {
+      z = z.square();
+      if (bits.codeUnitAt(i) == 0x31) z = z * this;
+    }
+    return z;
+  }
+
+  /// A square root, or null when this is not a square. Which of the two
+  /// roots it returns is not specified: the decoders choose the sign from
+  /// the flags of the encoding.
+  ///
+  /// The complex method: with N = c0² + c1² the norm, this is a square in
+  /// Fp2 if and only if N is a square in Fp; then one of (c0 ± √N)/2 is a
+  /// square x0² in Fp (their product, −c1²/4, is not, −1 being no square),
+  /// and x0 + c1/(2·x0)·u is a root. The result is checked by squaring.
+  Fp2? sqrt() {
+    if (c1.isZero) {
+      final s = c0.sqrt();
+      if (s != null) return Fp2(s, Fp.zero);
+      // c0 is not a square, so −c0 is: (t·u)² = −t² = c0.
+      final t = (-c0).sqrt();
+      return t == null ? null : Fp2(Fp.zero, t);
+    }
+    final alpha = Fp.mulAdd(c0, c0, c1, c1).sqrt();
+    if (alpha == null) return null;
+    var x0 = ((c0 + alpha) * Fp.half).sqrt();
+    x0 ??= ((c0 - alpha) * Fp.half).sqrt();
+    if (x0 == null || x0.isZero) return null;
+    final y = Fp2(x0, c1 * x0.double().inverse());
+    return y.square().equals(this) ? y : null;
+  }
+
+  /// Whether this is larger than its negation in the order of spec §12.2,
+  /// which sets the sign bit of a compressed point of G2: c1 decides, and c0
+  /// when c1 is 0.
+  bool get isLexicographicallyLargest =>
+      c1.isZero ? c0.isLexicographicallyLargest : c1.isLexicographicallyLargest;
+
+  /// The 96 bytes c1 || c0, the order of kilic/bls12-381 and of the
+  /// compressed points of G2 (spec §12.2).
+  Uint8List toBytes() => concatBytes([c1.toBytes(), c0.toBytes()]);
+}
+
+/// An element c0 + c1·v + c2·v² of Fp6.
+final class Fp6 {
+  /// The element [c0] + [c1]·v + [c2]·v².
+  const Fp6(this.c0, this.c1, this.c2);
+
+  /// 0.
+  static final Fp6 zero = Fp6(Fp2.zero, Fp2.zero, Fp2.zero);
+
+  /// 1.
+  static final Fp6 one = Fp6(Fp2.one, Fp2.zero, Fp2.zero);
+
+  /// The coefficient of 1.
+  final Fp2 c0;
+
+  /// The coefficient of v.
+  final Fp2 c1;
+
+  /// The coefficient of v².
+  final Fp2 c2;
+
+  /// Whether this is 0.
+  bool get isZero => c0.isZero && c1.isZero && c2.isZero;
+
+  /// Whether this is 1.
+  bool get isOne => c0.isOne && c1.isZero && c2.isZero;
+
+  /// Whether this and [other] are the same element.
+  bool equals(Fp6 other) =>
+      c0.equals(other.c0) && c1.equals(other.c1) && c2.equals(other.c2);
+
+  /// The sum.
+  Fp6 operator +(Fp6 o) => Fp6(c0 + o.c0, c1 + o.c1, c2 + o.c2);
+
+  /// The difference.
+  Fp6 operator -(Fp6 o) => Fp6(c0 - o.c0, c1 - o.c1, c2 - o.c2);
+
+  /// The negation.
+  Fp6 operator -() => Fp6(-c0, -c1, -c2);
+
+  /// The product, Algorithm 5.21 of the Guide to Pairing-Based
+  /// Cryptography, as kilic.
+  Fp6 operator *(Fp6 b) {
+    final v0 = c0 * b.c0;
+    final v1 = c1 * b.c1;
+    final v2 = c2 * b.c2;
+    return Fp6(
+      ((c1 + c2) * (b.c1 + b.c2) - v1 - v2).mulByNonResidue() + v0,
+      (c0 + c1) * (b.c0 + b.c1) - v0 - v1 + v2.mulByNonResidue(),
+      (c0 + c2) * (b.c0 + b.c2) - v0 - v2 + v1,
+    );
+  }
+
+  /// The square, CH-SQR2 of eprint 2006/471, as kilic.
+  Fp6 square() {
+    final s0 = c0.square();
+    final s1 = (c0 * c1).double();
+    final s2 = (c0 - c1 + c2).square();
+    final s3 = (c1 * c2).double();
+    final s4 = c2.square();
+    return Fp6(
+      s0 + s3.mulByNonResidue(),
+      s1 + s4.mulByNonResidue(),
+      s1 + s2 + s3 - s0 - s4,
+    );
+  }
+
+  /// The product by v, the non-residue of Fp12: (ξ·c2, c0, c1).
+  Fp6 mulByNonResidue() => Fp6(c2.mulByNonResidue(), c0, c1);
+
+  /// The product by b0 + b1·v, as mul01 of kilic.
+  Fp6 mul01(Fp2 b0, Fp2 b1) {
+    final v0 = c0 * b0;
+    final v1 = c1 * b1;
+    return Fp6(
+      ((c1 + c2) * b1 - v1).mulByNonResidue() + v0,
+      (c0 + c1) * (b0 + b1) - v0 - v1,
+      (c0 + c2) * b0 - v0 + v1,
+    );
+  }
+
+  /// The product by b1·v, as mul1 of kilic.
+  Fp6 mul1(Fp2 b1) => Fp6((c2 * b1).mulByNonResidue(), c0 * b1, c1 * b1);
+
+  /// The inverse, Algorithm 5.23 of the Guide to Pairing-Based
+  /// Cryptography, as kilic.
+  Fp6 inverse() {
+    final a = c0.square() - (c1 * c2).mulByNonResidue();
+    final b = c2.square().mulByNonResidue() - c0 * c1;
+    final c = c1.square() - c0 * c2;
+    final f = ((c2 * b + c1 * c).mulByNonResidue() + c0 * a).inverse();
+    return Fp6(a * f, b * f, c * f);
+  }
+
+  /// The 288 bytes c2 || c1 || c0, the order of kilic/bls12-381.
+  Uint8List toBytes() =>
+      concatBytes([c2.toBytes(), c1.toBytes(), c0.toBytes()]);
+}
+
+/// An element c0 + c1·w of Fp12; GT, the target group of the pairing, is its
+/// subgroup of order r.
+final class Fp12 {
+  /// The element [c0] + [c1]·w.
+  const Fp12(this.c0, this.c1);
+
+  /// 1.
+  static final Fp12 one = Fp12(Fp6.one, Fp6.zero);
+
+  /// The coefficient of 1.
+  final Fp6 c0;
+
+  /// The coefficient of w.
+  final Fp6 c1;
+
+  /// Whether this is 1.
+  bool get isOne => c0.isOne && c1.isZero;
+
+  /// Whether this and [other] are the same element.
+  bool equals(Fp12 other) => c0.equals(other.c0) && c1.equals(other.c1);
+
+  /// The product, Karatsuba over Fp6, as kilic.
+  Fp12 operator *(Fp12 b) {
+    final v0 = c0 * b.c0;
+    final v1 = c1 * b.c1;
+    return Fp12(v0 + v1.mulByNonResidue(), (c0 + c1) * (b.c0 + b.c1) - v0 - v1);
+  }
+
+  /// The square, the complex method of eprint 2006/471, as kilic.
+  Fp12 square() {
+    final t2 = c0 * c1;
+    final t = (c0 + c1) * (c1.mulByNonResidue() + c0);
+    return Fp12(t - t2 - t2.mulByNonResidue(), t2 + t2);
+  }
+
+  /// The conjugate c0 − c1·w, which is also this^(p⁶): the inverse of an
+  /// element of the cyclotomic subgroup, where GT lies.
+  Fp12 conjugate() => Fp12(c0, -c1);
+
+  /// The inverse, as kilic.
+  Fp12 inverse() {
+    final t = (c0.square() - c1.square().mulByNonResidue()).inverse();
+    return Fp12(c0 * t, -(c1 * t));
+  }
+
+  /// The product by the sparse b0 + b1·v + b4·v·w of a line of the Miller
+  /// loop, as mul014 of kilic.
+  Fp12 mul014(Fp2 b0, Fp2 b1, Fp2 b4) {
+    final t0 = c0.mul01(b0, b1);
+    final t1 = c1.mul1(b4);
+    return Fp12(
+      t1.mulByNonResidue() + t0,
+      (c1 + c0).mul01(b0, b1 + b4) - t0 - t1,
+    );
+  }
+
+  /// The square of an element of the cyclotomic subgroup, Granger and
+  /// Scott, as cyclotomicSquare of kilic. Wrong for any other element.
+  Fp12 cyclotomicSquare() {
+    final (t3, t4) = _fp4Square(c0.c0, c1.c1);
+    final r00 = (t3 - c0.c0).double() + t3;
+    final r11 = (t4 + c1.c1).double() + t4;
+    final (u3, u4) = _fp4Square(c1.c0, c0.c2);
+    final (u5, u6) = _fp4Square(c0.c1, c1.c2);
+    final r01 = (u3 - c0.c1).double() + u3;
+    final r12 = (u4 + c1.c2).double() + u4;
+    final w3 = u6.mulByNonResidue();
+    final r10 = (w3 + c1.c0).double() + w3;
+    final r02 = (u5 - c0.c2).double() + u5;
+    return Fp12(Fp6(r00, r01, r02), Fp6(r10, r11, r12));
+  }
+
+  /// This^(p^k), k in 1..3: the coefficient of w^i is conjugated k times and
+  /// multiplied by γ_k^i, with γ_k = ξ^((p^k − 1)/6).
+  Fp12 frobenius(int k) {
+    final g = _frobeniusPowers[k - 1];
+    Fp2 c(Fp2 a) => k.isOdd ? a.conjugate() : a;
+    return Fp12(
+      Fp6(c(c0.c0), c(c0.c1) * g[1], c(c0.c2) * g[3]),
+      Fp6(c(c1.c0) * g[0], c(c1.c1) * g[2], c(c1.c2) * g[4]),
+    );
+  }
+
+  /// The 576 bytes c1 || c0 in the order of kilic/bls12-381, which
+  /// drand/kyber marshals and H2 of the tlock IBE hashes (spec §63,
+  /// «Serialización de GT en H2»).
+  Uint8List toBytes() => concatBytes([c1.toBytes(), c0.toBytes()]);
+}
+
+// The square of a0 + a1·t in Fp4 = Fp2[t]/(t² − ξ), as fp4Square of kilic.
+(Fp2, Fp2) _fp4Square(Fp2 a0, Fp2 a1) {
+  final t0 = a0.square();
+  final t1 = a1.square();
+  return (t1.mulByNonResidue() + t0, (a0 + a1).square() - t0 - t1);
+}
+
+/// γ_k = ξ^((p^k − 1)/6) for k = 1, 2, 3, the coefficients of the Frobenius
+/// maps of Fp12 (frobeniusCoeffs12 of kilic). A test computes them again
+/// from their definition.
+final List<Fp2> frobeniusGammas = [
+  Fp2.hex(
+    '1904d3bf02bb0667c231beb4202c0d1f0fd603fd3cbd5f4f7b2443d784bab9c4f67ea53d'
+        '63e7813d8d0775ed92235fb8',
+    '00fc3e2b36c4e03288e9e902231f9fb854a14787b6c7b36fec0c8ec971f63c5f282d5ac1'
+        '4d6c7ec22cf78a126ddc4af3',
+  ),
+  Fp2.hex(
+    '00000000000000005f19672fdf76ce51ba69c6076a0f77eaddb3a93be6f89688de17d813'
+        '620a00022e01fffffffeffff',
+    '0',
+  ),
+  Fp2.hex(
+    '135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60ef396489f61eb45e304466cf'
+        '3e67fa0af1ee7b04121bdea2',
+    '06af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e77f76e17009241c5ee67992f'
+        '72ec05f4c81084fbede3cc09',
+  ),
+];
+
+// γ_k^1 .. γ_k^5 for k = 1, 2, 3.
+final List<List<Fp2>> _frobeniusPowers = [
+  for (final g in frobeniusGammas)
+    () {
+      final out = [g];
+      for (var i = 1; i < 5; i++) {
+        out.add(out[i - 1] * g);
+      }
+      return out;
+    }(),
+];

test/bls12381_constants.dart

@@ -0,0 +1,36 @@
+// Go reference values of BLS12-381 that the tests running on Node.js use:
+// they read no file. Each is copied from the JSON that Go wrote,
+// testdata/vectors/tlock_ibe.json or test/vectors/bls12381_vectors.json, and
+// bls12381_constants_test.dart checks it against that JSON on the VM.
+library;
+
+/// The generator of G1, compressed (tlock_ibe.json).
+const generatorG1 =
+    '97f1d3a73197d7942695638c4fa9ac0fc3688c4f9774b905a14e3a3f171bac586c55e83f'
+    'f97a1aeffb3af00adb22c6bb';
+
+/// The generator of G2, compressed (tlock_ibe.json).
+const generatorG2 =
+    '93e02b6052719f607dacd3a088274f65596bd0d09920b61ab5da61bbdc7f5049334cf112'
+    '13945d57e5ac7d055d042b7e024aa2b2f08f0a91260805272dc51051c6e47ad4fa403b02'
+    'b4510b647ae3d1770bac0326a805bbefd48056c8c121bdb8';
+
+/// H2(e(G1, G2)) truncated to 16 bytes (tlock_ibe.json).
+const h2OfGenerators = 'cb87319f24560b5231579a09ad79f12e';
+
+/// The identity of round 1000 hashed to G1 with the DST of Quicknet
+/// (bls12381_vectors.json).
+const hashOfRound1000 =
+    '8f5a32d53837b00fbc0ee31ce9966435a41c5188a80ce9934d3c80588b6ad6f643ebda1b'
+    '83ef89e44da9ced6205cdecf';
+
+/// hash_to_curve of "abc" with the DST of RFC 9380, appendix J.9.1: x
+/// (bls12381_vectors.json).
+const hashOfAbcX =
+    '03567bc5ef9c690c2ab2ecdf6a96ef1c139cc0b2f284dca0a9a7943388a49a3aee664ba5'
+    '379a7655d3c68900be2f6903';
+
+/// The same: y.
+const hashOfAbcY =
+    '0b9c15f3fe6e5cf4211f346271d7b01c8f3b28be689c8429c85b67af215533311f0b8dfa'
+    'aa154fa6b88176c229f2885d';

test/bls12381_constants_test.dart

@@ -0,0 +1,46 @@
+// The Go values that the tests on Node.js carry in bls12381_constants.dart
+// are those of the JSON that Go wrote.
+@TestOn('vm')
+library;
+
+import 'dart:convert';
+import 'dart:io';
+
+import 'package:test/test.dart';
+
+import 'bls12381_constants.dart';
+
+typedef Json = Map<String, Object?>;
+
+Json readJson(String path) => jsonDecode(File(path).readAsStringSync()) as Json;
+
+List<Json> section(Json file, String name) =>
+    (file[name]! as List).cast<Json>();
+
+void main() {
+  test('the generators and H2 are those of tlock_ibe.json', () {
+    final v = section(
+      readJson('testdata/vectors/tlock_ibe.json'),
+      'vectors',
+    ).first;
+    expect(
+      [generatorG1, generatorG2, h2OfGenerators],
+      [v['g1'], v['g2'], v['h2']],
+    );
+  });
+
+  test('the hashes to G1 are those of bls12381_vectors.json', () {
+    final hashes = section(
+      readJson('test/vectors/bls12381_vectors.json'),
+      'hash_to_g1',
+    );
+    final round = hashes.firstWhere(
+      (v) => v['label'] == 'the identity of round 1000',
+    );
+    expect(hashOfRound1000, round['point']);
+    final abc = hashes.firstWhere(
+      (v) => (v['label']! as String).contains('msg "abc"'),
+    );
+    expect([hashOfAbcX, hashOfAbcY], [abc['x'], abc['y']]);
+  });
+}

test/bls12381_test.dart

@@ -0,0 +1,632 @@
+// The arithmetic of BLS12-381 by its properties, on the VM and compiled to
+// JavaScript: the field layer against plain BigInt arithmetic, the tower,
+// the constants against their definitions, the group laws and the encodings,
+// the subgroup checks, the pairing and the hash to G1. They read no file:
+// the few Go values they use are in bls12381_constants.dart. The vectors of Go
+// are in bls12381_vectors_test.dart.
+//
+// On Node.js the BigInt of dart2js is about 20 times slower than on the VM,
+// so the loops draw fewer cases there ([web]).
+library;
+
+import 'dart:math';
+import 'dart:typed_data';
+
+import 'package:crypto/crypto.dart';
+import 'package:datekeys/datekeys.dart' show fromHex, toHex;
+import 'package:datekeys/src/bls12381_curve.dart';
+import 'package:datekeys/src/bls12381_fp.dart';
+import 'package:datekeys/src/bls12381_hash.dart';
+import 'package:datekeys/src/bls12381_pairing.dart';
+import 'package:datekeys/src/bls12381_tower.dart';
+import 'package:test/test.dart';
+
+import 'bls12381_constants.dart';
+
+/// Whether the tests run compiled to JavaScript.
+final bool web = identical(0, 0.0);
+
+/// [vm] cases on the VM, [js] on the web.
+int cases(int vm, int js) => web ? js : vm;
+
+final BigInt p = fpModulus;
+final BigInt r = groupOrder;
+
+// A fixed-seed source of field elements, scalars and bytes.
+final class Draw {
+  Draw(int seed) : _r = Random(seed);
+
+  final Random _r;
+
+  Uint8List bytes(int n) =>
+      Uint8List.fromList([for (var i = 0; i < n; i++) _r.nextInt(256)]);
+
+  BigInt big(int n) => BigInt.parse(toHex(bytes(n)), radix: 16);
+
+  Fp fp() => Fp.fromBytesReduced(bytes(64));
+
+  Fp2 fp2() => Fp2(fp(), fp());
+
+  Fp6 fp6() => Fp6(fp2(), fp2(), fp2());
+
+  Fp12 fp12() => Fp12(fp6(), fp6());
+
+  /// A scalar in 1..r-1.
+  BigInt scalar() => big(40) % (r - BigInt.one) + BigInt.one;
+
+  /// A point of E(Fp) from a random x: almost never in G1.
+  G1Point g1OnCurve() {
+    for (;;) {
+      final x = fp();
+      final y = (x.square() * x + G1Point.b).sqrt();
+      if (y != null) return G1Point.affine(x, y);
+    }
+  }
+
+  /// A point of E'(Fp2) from a random x: almost never in G2.
+  G2Point g2OnCurve() {
+    for (;;) {
+      final x = fp2();
+      final y = (x.square() * x + G2Point.b).sqrt();
+      if (y != null) return G2Point.affine(x, y);
+    }
+  }
+}
+
+BigInt big(Fp a) => a.toBigInt();
+
+Fp fp(int v) => Fp(BigInt.from(v));
+
+// The product in Fp2 by the definition u² = −1.
+Fp2 mulFp2(Fp2 a, Fp2 b) =>
+    Fp2(a.c0 * b.c0 - a.c1 * b.c1, a.c0 * b.c1 + a.c1 * b.c0);
+
+final Fp2 xi = Fp2(Fp.one, Fp.one);
+
+// The product in Fp6 by the definition v³ = ξ, schoolbook.
+Fp6 mulFp6(Fp6 a, Fp6 b) {
+  final c = List.filled(5, Fp2.zero);
+  final x = [a.c0, a.c1, a.c2];
+  final y = [b.c0, b.c1, b.c2];
+  for (var i = 0; i < 3; i++) {
+    for (var j = 0; j < 3; j++) {
+      c[i + j] = c[i + j] + mulFp2(x[i], y[j]);
+    }
+  }
+  return Fp6(c[0] + mulFp2(c[3], xi), c[1] + mulFp2(c[4], xi), c[2]);
+}
+
+final Fp6 v6 = Fp6(Fp2.zero, Fp2.one, Fp2.zero);
+
+// The product in Fp12 by the definition w² = v, schoolbook.
+Fp12 mulFp12(Fp12 a, Fp12 b) => Fp12(
+  mulFp6(a.c0, b.c0) + mulFp6(mulFp6(a.c1, b.c1), v6),
+  mulFp6(a.c0, b.c1) + mulFp6(a.c1, b.c0),
+);
+
+void main() {
+  group('Fp, the field layer', () {
+    test('p and r are those of spec §12.2, and p = 3 mod 4', () {
+      expect(
+        p.toRadixString(16),
+        '1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241ea'
+        'bfffeb153ffffb9feffffffffaaab',
+      );
+      expect(
+        r.toRadixString(16),
+        '73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001',
+      );
+      expect(p % BigInt.from(4), BigInt.from(3));
+    });
+
+    test('computes as BigInt arithmetic modulo p', () {
+      final d = Draw(1);
+      for (var i = 0; i < cases(200, 20); i++) {
+        final a = d.fp();
+        final b = d.fp();
+        final c = d.fp();
+        final e = d.fp();
+        expect(big(a + b), (big(a) + big(b)) % p);
+        expect(big(a - b), (big(a) - big(b)) % p);
+        expect(big(-a), (-big(a)) % p);
+        expect(big(a * b), (big(a) * big(b)) % p);
+        expect(big(a.square()), (big(a) * big(a)) % p);
+        expect(big(a.double()), (big(a) * BigInt.two) % p);
+        expect(
+          big(Fp.mulSub(a, b, c, e)),
+          (big(a) * big(b) - big(c) * big(e)) % p,
+        );
+        expect(
+          big(Fp.mulAdd(a, b, c, e)),
+          (big(a) * big(b) + big(c) * big(e)) % p,
+        );
+        expect((a * a.inverse()).isOne, isTrue);
+        expect(a.pow(BigInt.from(5)).equals(a * a * a * a * a), isTrue);
+      }
+      expect(Fp.zero.inverse().isZero, isTrue);
+      expect((-Fp.zero).isZero, isTrue);
+      expect((Fp.half + Fp.half).isOne, isTrue);
+    });
+
+    test('takes square roots exactly of the squares', () {
+      final d = Draw(2);
+      var squares = 0;
+      for (var i = 0; i < cases(60, 8); i++) {
+        final a = d.fp();
+        final legendre = big(a).modPow((p - BigInt.one) >> 1, p);
+        final s = a.sqrt();
+        expect(s != null, legendre == BigInt.one, reason: '$i');
+        if (s != null) {
+          expect(s.square().equals(a), isTrue);
+          squares++;
+        } else {
+          // The power is then a root of −a.
+          final (t, isSquare) = a.sqrtOrNegatedRoot();
+          expect(isSquare, isFalse);
+          expect(t.square().equals(-a), isTrue);
+        }
+      }
+      expect(squares, inInclusiveRange(1, cases(59, 7)));
+      expect(Fp.zero.sqrt()!.isZero, isTrue);
+    });
+
+    test('reads and writes 48 big-endian bytes, never reducing', () {
+      final d = Draw(3);
+      for (var i = 0; i < 20; i++) {
+        final a = d.fp();
+        expect(Fp.fromBytes(a.toBytes())!.equals(a), isTrue);
+      }
+      Uint8List be(BigInt v) => fromHex(v.toRadixString(16).padLeft(96, '0'));
+      expect(Fp.fromBytes(be(p - BigInt.one))!.toBigInt(), p - BigInt.one);
+      expect(Fp.fromBytes(be(p)), isNull);
+      expect(Fp.fromBytes(be(p + BigInt.one)), isNull);
+      expect(Fp.fromBytes(Uint8List(48)..fillRange(0, 48, 0xff)), isNull);
+      expect(() => Fp.fromBytes(Uint8List(47)), throwsArgumentError);
+      expect(Fp.fromBytesReduced(be(p + BigInt.two)).toBigInt(), BigInt.two);
+      expect(() => Fp(p), throwsArgumentError);
+      expect(() => Fp(-BigInt.one), throwsArgumentError);
+    });
+
+    test('the sign of spec §12.2 turns at (p − 1)/2', () {
+      final half = (p - BigInt.one) >> 1;
+      expect(Fp(half).isLexicographicallyLargest, isFalse);
+      expect(Fp(half + BigInt.one).isLexicographicallyLargest, isTrue);
+      expect(Fp.zero.isLexicographicallyLargest, isFalse);
+      expect(fp(3).isOdd, isTrue);
+      expect(fp(4).isOdd, isFalse);
+    });
+  });
+
+  group('the tower', () {
+    test('Fp2, Fp6 and Fp12 multiply as their definitions', () {
+      final d = Draw(4);
+      for (var i = 0; i < cases(20, 3); i++) {
+        final a2 = d.fp2();
+        final b2 = d.fp2();
+        expect((a2 * b2).equals(mulFp2(a2, b2)), isTrue);
+        expect(a2.square().equals(a2 * a2), isTrue);
+        expect((a2 * a2.inverse()).isOne, isTrue);
+        expect(a2.mulByNonResidue().equals(a2 * xi), isTrue);
+        final a6 = d.fp6();
+        final b6 = d.fp6();
+        expect((a6 * b6).equals(mulFp6(a6, b6)), isTrue);
+        expect(a6.square().equals(a6 * a6), isTrue);
+        expect((a6 * a6.inverse()).isOne, isTrue);
+        expect(a6.mulByNonResidue().equals(a6 * v6), isTrue);
+        final x0 = d.fp2();
+        final x1 = d.fp2();
+        expect(a6.mul01(x0, x1).equals(a6 * Fp6(x0, x1, Fp2.zero)), isTrue);
+        expect(a6.mul1(x1).equals(a6 * Fp6(Fp2.zero, x1, Fp2.zero)), isTrue);
+        final a12 = d.fp12();
+        final b12 = d.fp12();
+        expect((a12 * b12).equals(mulFp12(a12, b12)), isTrue);
+        expect(a12.square().equals(a12 * a12), isTrue);
+        expect((a12 * a12.inverse()).isOne, isTrue);
+        final x4 = d.fp2();
+        final sparse = Fp12(Fp6(x0, x1, Fp2.zero), Fp6(Fp2.zero, x4, Fp2.zero));
+        expect(a12.mul014(x0, x1, x4).equals(a12 * sparse), isTrue);
+      }
+    });
+
+    test('takes square roots in Fp2 exactly of the squares', () {
+      final d = Draw(5);
+      for (var i = 0; i < cases(30, 4); i++) {
+        final a = d.fp2();
+        final sq = a.square();
+        final s = sq.sqrt()!;
+        expect(s.square().equals(sq), isTrue);
+        // ξ = 1 + u is not a square: neither is ξ·a².
+        expect((sq * xi).sqrt(), isNull);
+      }
+      // Elements of Fp: a square, and a non-square whose root is in u·Fp.
+      final four = Fp2(fp(4), Fp.zero);
+      expect(four.sqrt()!.square().equals(four), isTrue);
+      final minusFour = Fp2(-fp(4), Fp.zero);
+      final root = minusFour.sqrt()!;
+      expect(root.c0.isZero, isTrue);
+      expect(root.square().equals(minusFour), isTrue);
+      expect(Fp2.zero.sqrt()!.isZero, isTrue);
+    });
+
+    test('the Frobenius coefficients are ξ^((p^k − 1)/6)', () {
+      for (var k = 1; k <= 3; k++) {
+        final e = (p.pow(k) - BigInt.one) ~/ BigInt.from(6);
+        expect(xi.pow(e).equals(frobeniusGammas[k - 1]), isTrue, reason: '$k');
+      }
+    });
+
+    test('the Frobenius map of Fp12 is a ring morphism of order 12', () {
+      final d = Draw(6);
+      final a = d.fp12();
+      final b = d.fp12();
+      for (var k = 1; k <= 3; k++) {
+        expect(
+          (a * b).frobenius(k).equals(a.frobenius(k) * b.frobenius(k)),
+          isTrue,
+          reason: '$k',
+        );
+      }
+      expect(a.frobenius(1).frobenius(1).equals(a.frobenius(2)), isTrue);
+      expect(a.frobenius(2).frobenius(1).equals(a.frobenius(3)), isTrue);
+      var x = a;
+      for (var i = 0; i < 12; i++) {
+        x = x.frobenius(1);
+      }
+      expect(x.equals(a), isTrue);
+      // On an element of Fp2, it is the conjugation.
+      final c = d.fp2();
+      final f = Fp12(Fp6(c, Fp2.zero, Fp2.zero), Fp6.zero);
+      expect(f.frobenius(1).c0.c0.equals(c.conjugate()), isTrue);
+      // And it is the p-th power: (a^p)·a^-p = 1 needs a^p, which the
+      // pairing vectors of Go check through the final exponentiation.
+    });
+
+    test('squares in the cyclotomic subgroup as in the field', () {
+      final d = Draw(7);
+      for (var i = 0; i < cases(5, 2); i++) {
+        // The easy part of the final exponentiation lands in the subgroup.
+        final f = d.fp12();
+        var g = f.conjugate() * f.inverse();
+        g = g.frobenius(2) * g;
+        expect(g.cyclotomicSquare().equals(g.square()), isTrue);
+        expect((g * g.conjugate()).isOne, isTrue);
+      }
+    });
+
+    test('GT is written c1 before c0 at every level', () {
+      final e = [for (var i = 1; i <= 12; i++) fp(i)];
+      final f = Fp12(
+        Fp6(Fp2(e[0], e[1]), Fp2(e[2], e[3]), Fp2(e[4], e[5])),
+        Fp6(Fp2(e[6], e[7]), Fp2(e[8], e[9]), Fp2(e[10], e[11])),
+      );
+      final b = f.toBytes();
+      expect(b, hasLength(576));
+      // c1 of Fp12: c2, c1, c0 of Fp6, each c1 then c0 of Fp2.
+      expect(
+        [for (var i = 0; i < 12; i++) b[48 * i + 47]],
+        [
+          12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, //
+        ],
+      );
+    });
+  });
+
+  group('the curves', () {
+    test('the generators are those of Go, on their curves and in their '
+        'subgroups', () {
+      expect(toHex(G1Point.generator.toBytes()), generatorG1);
+      expect(toHex(G2Point.generator.toBytes()), generatorG2);
+      expect(
+        G1Point.decode(fromHex(generatorG1))!.equals(G1Point.generator),
+        isTrue,
+      );
+      expect(
+        G2Point.decode(fromHex(generatorG2))!.equals(G2Point.generator),
+        isTrue,
+      );
+      final (x1, y1) = G1Point.generator.toAffine()!;
+      final (x2, y2) = G2Point.generator.toAffine()!;
+      expect(G1Point.isOnCurve(x1, y1), isTrue);
+      expect(G2Point.isOnCurve(x2, y2), isTrue);
+      expect(G1Point.generator.multiply(r).isInfinity, isTrue);
+      expect(G2Point.generator.multiply(r).isInfinity, isTrue);
+    });
+
+    test('add, double and multiply as a group of order r', () {
+      final d = Draw(8);
+      for (var i = 0; i < cases(4, 1); i++) {
+        final a = d.scalar();
+        final b = d.scalar();
+        final g1 = G1Point.generator;
+        final g2 = G2Point.generator;
+        expect(
+          (g1.multiply(a) + g1.multiply(b)).equals(g1.multiply((a + b) % r)),
+          isTrue,
+        );
+        expect(
+          (g2.multiply(a) + g2.multiply(b)).equals(g2.multiply((a + b) % r)),
+          isTrue,
+        );
+        final p1 = g1.multiply(a);
+        final p2 = g2.multiply(b);
+        expect((p1 + p1).equals(p1.double()), isTrue);
+        expect((p2 + p2).equals(p2.double()), isTrue);
+        expect((p1 + -p1).isInfinity, isTrue);
+        expect((p2 + -p2).isInfinity, isTrue);
+        expect(p1.multiply(r + BigInt.one).equals(p1), isTrue);
+        expect(p2.multiply(BigInt.zero).isInfinity, isTrue);
+        final (ax, ay) = p1.toAffine()!;
+        expect(g1.addAffine(ax, ay).equals(g1 + p1), isTrue);
+        final (bx, by) = p2.toAffine()!;
+        expect(g2.addAffine(bx, by).equals(g2 + p2), isTrue);
+      }
+      expect(
+        (G1Point.infinity + G1Point.generator).equals(G1Point.generator),
+        isTrue,
+      );
+      expect(G1Point.infinity.double().isInfinity, isTrue);
+      expect(
+        () => G1Point.generator.multiply(-BigInt.one),
+        throwsArgumentError,
+      );
+    });
+
+    test('encode and decode the canonical encodings only', () {
+      final d = Draw(9);
+      final p1 = G1Point.generator.multiply(d.scalar());
+      final p2 = G2Point.generator.multiply(d.scalar());
+      for (final (group, b, point) in [
+        (BlsGroup.g1, p1.toBytes(), p1 as Object),
+        (BlsGroup.g2, p2.toBytes(), p2 as Object),
+      ]) {
+        expect(checkCompressedPoint(group, b), PointVerdict.point);
+        final decoded = group == BlsGroup.g1
+            ? G1Point.decode(b)
+            : G2Point.decode(b);
+        expect(
+          decoded is G1Point
+              ? decoded.equals(point as G1Point)
+              : (decoded! as G2Point).equals(point as G2Point),
+          isTrue,
+        );
+        // The negation flips the sign bit only.
+        final neg = Uint8List.fromList(b)..[0] ^= 0x20;
+        final negated = group == BlsGroup.g1
+            ? (-(point as G1Point)).toBytes()
+            : (-(point as G2Point)).toBytes();
+        expect(negated, neg);
+        // Not canonical: the compression flag cleared, the infinity flag
+        // set, another length, x + p in the first coordinate.
+        expect(
+          checkCompressedPoint(group, Uint8List.fromList(b)..[0] ^= 0x80),
+          PointVerdict.invalid,
+        );
+        expect(
+          checkCompressedPoint(group, Uint8List.fromList(b)..[0] |= 0x40),
+          PointVerdict.invalid,
+        );
+        expect(checkCompressedPoint(group, b.sublist(1)), PointVerdict.invalid);
+        expect(checkCompressedPoint(group, [...b, 0]), PointVerdict.invalid);
+        final raw = Uint8List.fromList(b)..[0] &= 0x1f;
+        final x = BigInt.parse(toHex(raw.sublist(0, 48)), radix: 16) + p;
+        if (x.bitLength <= 381) {
+          final plusP = fromHex(x.toRadixString(16).padLeft(96, '0'));
+          plusP[0] |= b[0] & 0xe0;
+          expect(
+            checkCompressedPoint(group, [...plusP, ...b.sublist(48)]),
+            PointVerdict.invalid,
+          );
+        }
+      }
+      for (final group in BlsGroup.values) {
+        final n = group.byteLength;
+        final infinity = Uint8List(n)..[0] = 0xc0;
+        expect(checkCompressedPoint(group, infinity), PointVerdict.identity);
+        expect(
+          checkCompressedPoint(group, Uint8List.fromList(infinity)..[0] = 0xe0),
+          PointVerdict.invalid,
+        );
+        expect(
+          checkCompressedPoint(
+            group,
+            Uint8List.fromList(infinity)..[n - 1] = 1,
+          ),
+          PointVerdict.invalid,
+        );
+        expect(checkCompressedPoint(group, Uint8List(n)), PointVerdict.invalid);
+      }
+      expect(toHex(G1Point.infinity.toBytes()), 'c0${'00' * 47}');
+      expect(toHex(G2Point.infinity.toBytes()), 'c0${'00' * 95}');
+    });
+
+    test('reject the points of the curve outside the subgroup, torsion '
+        'added to a point of the subgroup included', () {
+      final d = Draw(10);
+      for (var i = 0; i < cases(3, 1); i++) {
+        final q1 = d.g1OnCurve();
+        expect(q1.isInSubgroup, isFalse);
+        expect(
+          checkCompressedPoint(BlsGroup.g1, q1.toBytes()),
+          PointVerdict.invalid,
+        );
+        // [r]·Q is a point of order dividing the cofactor.
+        final t1 = q1.multiply(r);
+        expect(t1.isInfinity, isFalse);
+        final s1 = G1Point.generator.multiply(d.scalar()) + t1;
+        expect(s1.isInSubgroup, isFalse);
+        expect(
+          checkCompressedPoint(BlsGroup.g1, s1.toBytes()),
+          PointVerdict.invalid,
+        );
+        final q2 = d.g2OnCurve();
+        expect(q2.isInSubgroup, isFalse);
+        expect(
+          checkCompressedPoint(BlsGroup.g2, q2.toBytes()),
+          PointVerdict.invalid,
+        );
+        final t2 = q2.multiply(r);
+        final s2 = G2Point.generator.multiply(d.scalar()) + t2;
+        expect(s2.isInSubgroup, isFalse);
+        expect(
+          checkCompressedPoint(BlsGroup.g2, s2.toBytes()),
+          PointVerdict.invalid,
+        );
+      }
+    });
+
+    test('the subgroup check of G2 by ψ decides as [r]·P = O, torsion of '
+        'every small order of the cofactor included', () {
+      final xi = Fp2(Fp.one, Fp.one);
+      expect(
+        xi.pow((p - BigInt.one) ~/ BigInt.from(3)).inverse().equals(psiX),
+        isTrue,
+      );
+      expect(
+        xi.pow((p - BigInt.one) ~/ BigInt.two).inverse().equals(psiY),
+        isTrue,
+      );
+      final g = G2Point.generator;
+      final absX = BigInt.parse('d201000000010000', radix: 16);
+      expect(g.psi().equals(-g.multiply(absX)), isTrue);
+      // The cofactor of G2 (cofactorG2 of kilic), and its prime factors
+      // below 2^21: 13², 23², 2713, 11953 and 262069.
+      final h2 = BigInt.parse(
+        '5d543a95414e7f1091d50792876a202cd91de4547085abaa68a205b2e5a7ddfa628'
+        'f1cb4d9e82ef21537e293a6691ae1616ec6e786f0c70cf1c38e31c7238e5',
+        radix: 16,
+      );
+      final small = [13, 23, 2713, 11953, 262069];
+      for (final f in small) {
+        expect(h2 % BigInt.from(f), BigInt.zero);
+      }
+      final d = Draw(13);
+      final points = <G2Point>[];
+      for (var i = 0; i < cases(3, 1); i++) {
+        final q = d.g2OnCurve();
+        // r·h2 is the order of E'(Fp2).
+        expect(q.multiply(r * h2).isInfinity, isTrue);
+        final s = g.multiply(d.scalar());
+        points.addAll([q, s, s + q.multiply(r)]);
+        // A point of each small order of the cofactor, and the point of the
+        // subgroup plus it.
+        for (final f in web ? small.take(2) : small) {
+          final t = q.multiply(r * h2 ~/ BigInt.from(f));
+          if (!t.isInfinity) points.addAll([t, s + t]);
+        }
+      }
+      var outside = 0;
+      for (final q in points) {
+        final naive = q.multiply(r).isInfinity;
+        expect(q.isInSubgroup, naive);
+        if (!naive) outside++;
+      }
+      expect(outside, greaterThanOrEqualTo(points.length * 2 ~/ 3));
+    });
+  });
+
+  group('the pairing', () {
+    test('e(G1, G2) gives the H2 of tlock_ibe.json', () {
+      // H2 of the IBE of tlock: SHA-256("IBE-H2" || GT), truncated.
+      final gt = pairing(G1Point.generator, G2Point.generator);
+      final h2 = sha256.convert([...'IBE-H2'.codeUnits, ...gt.toBytes()]);
+      expect(toHex(h2.bytes.sublist(0, 16)), h2OfGenerators);
+    });
+
+    test('is bilinear and non-degenerate', () {
+      final d = Draw(11);
+      final p1 = G1Point.generator.multiply(d.scalar());
+      final q2 = G2Point.generator.multiply(d.scalar());
+      final e = pairing(p1, q2);
+      expect(e.isOne, isFalse);
+      expect(pairing(p1.double(), q2).equals(e.square()), isTrue);
+      expect(pairing(p1, q2.double() + q2).equals(e.square() * e), isTrue);
+      expect(pairing(-p1, q2).equals(e.conjugate()), isTrue);
+      expect(pairing(G1Point.infinity, q2).isOne, isTrue);
+      expect(pairing(p1, G2Point.infinity).isOne, isTrue);
+      expect(pairingCheck([(p1, q2), (-p1, q2)]), isTrue);
+      expect(pairingCheck([(p1.double(), q2), (-p1, q2.double())]), isTrue);
+      expect(pairingCheck([(p1, q2), (p1, q2)]), isFalse);
+      expect(pairingCheck([(G1Point.infinity, q2)]), isTrue);
+    });
+
+    test('lands in the subgroup of order r of the cyclotomic subgroup', () {
+      final e = pairing(G1Point.generator, G2Point.generator);
+      var x = Fp12.one;
+      final bits = r.toRadixString(2);
+      for (var i = 0; i < bits.length; i++) {
+        x = x.cyclotomicSquare();
+        if (bits[i] == '1') x = x * e;
+      }
+      expect(x.isOne, isTrue);
+    });
+
+    test('runs the Miller loop over |x| = 0xd201000000010000', () {
+      expect(
+        BigInt.parse('1$millerLoopBits', radix: 2).toRadixString(16),
+        'd201000000010000',
+      );
+    });
+  });
+
+  group('the hash to G1', () {
+    test('maps to E′, which the isogeny maps to E, and lands in G1', () {
+      final d = Draw(12);
+      final a = Fp.hex(
+        '00144698a3b8e9433d693a02c96d4982b0ea985383ee66a8d8e8981aefd881ac9893'
+        '6f8da0e0f97f5cf428082d584c1d',
+      );
+      final b = Fp.hex(
+        '12e2908d11688030018b12e8753eee3b2016c1f0f24f4070a0b9c14fcef35ef55a23'
+        '215a316ceaa5d1cc48e98e172be0',
+      );
+      for (var i = 0; i < cases(10, 2); i++) {
+        final u = d.fp();
+        final (x, y) = mapToIsogenousCurve(u);
+        expect(y.square().equals((x.square() + a) * x + b), isTrue);
+        expect(y.isOdd, u.isOdd);
+        final (ix, iy) = isogenyMap(x, y);
+        expect(G1Point.isOnCurve(ix, iy), isTrue);
+      }
+      // The exceptional u = 0: x1 = B/(Z·A).
+      final (x0, _) = mapToIsogenousCurve(Fp.zero);
+      expect(x0.equals(b * (fp(11) * a).inverse()), isTrue);
+      expect(isogenyConstants.map((k) => k.length), [12, 11, 16, 16]);
+      expect(isogenyConstants[1].last.isOne, isTrue);
+      expect(isogenyConstants[3].last.isOne, isTrue);
+    });
+
+    test('hashes as Go, for the DST of Quicknet and of RFC 9380', () {
+      // The identity of round 1000: SHA-256 of its 8 big-endian bytes.
+      final id = sha256.convert([0, 0, 0, 0, 0, 0, 0x03, 0xe8]).bytes;
+      final p1 = hashToG1(id, quicknetDst);
+      expect(toHex(p1.toBytes()), hashOfRound1000);
+      expect(p1.isInSubgroup, isTrue);
+      final (x, y) = hashToG1(
+        'abc'.codeUnits,
+        'QUUX-V01-CS02-with-BLS12381G1_XMD:SHA-256_SSWU_RO_',
+      ).toAffine()!;
+      expect(
+        [toHex(x.toBytes()), toHex(y.toBytes())],
+        [hashOfAbcX, hashOfAbcY],
+      );
+    });
+
+    test('expand_message_xmd checks its lengths', () {
+      expect(expandMessageXmd([], 'DST'.codeUnits, 0), isEmpty);
+      expect(expandMessageXmd([], 'DST'.codeUnits, 8160), hasLength(8160));
+      expect(
+        () => expandMessageXmd([], 'DST'.codeUnits, 8161),
+        throwsArgumentError,
+      );
+      expect(
+        () => expandMessageXmd([], List.filled(256, 0x41), 32),
+        throwsArgumentError,
+      );
+      // The length is an input of b_0: no output is a prefix of a longer
+      // one.
+      final a = expandMessageXmd([1, 2], 'DST'.codeUnits, 32);
+      final b = expandMessageXmd([1, 2], 'DST'.codeUnits, 64);
+      expect(a, isNot(b.sublist(0, 32)));
+    });
+  });
+}

test/bls12381_vectors_test.dart

@@ -0,0 +1,219 @@
+// The BLS12-381 code against the Go reference: test/vectors/
+// bls12381_vectors.json, written by tool/bls12381_go_vectors.go with
+// kilic/bls12-381 and drand/kyber-bls12381, the libraries of drand and
+// tlock. Decoding verdicts, sums, multiples, pairings, hashes to G1, the map
+// of one element and BLS signatures on G1.
+@TestOn('vm')
+library;
+
+import 'dart:convert';
+import 'dart:io';
+
+import 'package:datekeys/datekeys.dart' show fromHex, toHex;
+import 'package:datekeys/src/bls12381_curve.dart';
+import 'package:datekeys/src/bls12381_fp.dart';
+import 'package:datekeys/src/bls12381_hash.dart';
+import 'package:datekeys/src/bls12381_pairing.dart';
+import 'package:datekeys/src/bls12381_tower.dart';
+import 'package:test/test.dart';
+
+typedef Json = Map<String, Object?>;
+
+final Json vectors = jsonDecode(
+  File('test/vectors/bls12381_vectors.json').readAsStringSync(),
+) as Json;
+
+List<Json> section(String name) => (vectors[name]! as List).cast<Json>();
+
+String s(Json v, String key) => v[key]! as String;
+
+BlsGroup group(Json v) => s(v, 'group') == 'G1' ? BlsGroup.g1 : BlsGroup.g2;
+
+PointVerdict verdictOf(String go) => PointVerdict.values.byName(go);
+
+G1Point g1(String hex) => G1Point.decode(fromHex(hex))!;
+
+G2Point g2(String hex) => G2Point.decode(fromHex(hex))!;
+
+void main() {
+  test('the vectors come from the Go libraries of the reference', () {
+    expect(vectors['generator'], 'tool/bls12381_go_vectors.go');
+    expect(
+      vectors['libraries'],
+      allOf(
+        contains('github.com/kilic/bls12-381 v0.1.0'),
+        contains('github.com/drand/kyber-bls12381 v0.3.4'),
+      ),
+    );
+    expect(section('points'), hasLength(157));
+  });
+
+  test('decodes every frozen edge case of datekeys-ts as Go does', () {
+    final seen = <String>{};
+    for (final v in section('points')) {
+      final got = checkCompressedPoint(group(v), fromHex(s(v, 'hex')));
+      expect(got, verdictOf(s(v, 'go')), reason: s(v, 'label'));
+      seen.add('${s(v, 'group')} ${s(v, 'go')} ${v['class'] ?? ''}');
+    }
+    // Every class of verdict and of failure, in both groups.
+    expect(
+      seen,
+      containsAll([
+        'G1 point ',
+        'G1 identity ',
+        'G1 invalid format',
+        'G1 invalid curve',
+        'G1 invalid subgroup',
+        'G2 point ',
+        'G2 identity ',
+        'G2 invalid format',
+        'G2 invalid curve',
+        'G2 invalid subgroup',
+      ]),
+    );
+  });
+
+  test('decodes the encodings drawn from the seed as Go does, failing where Go '
+      'fails', () {
+    final classes = <String, int>{};
+    for (final v in section('decode')) {
+      final b = fromHex(s(v, 'hex'));
+      final label = s(v, 'label');
+      expect(
+        checkCompressedPoint(group(v), b),
+        verdictOf(s(v, 'go')),
+        reason: label,
+      );
+      final cls = v['class'] as String?;
+      final key = '${s(v, 'group')} ${cls ?? s(v, 'go')}';
+      classes[key] = (classes[key] ?? 0) + 1;
+      // Where Go finds a point of the curve outside the subgroup, so does
+      // this code: the subgroup check is what rejects it.
+      if (cls == 'subgroup' || cls == 'curve') {
+        expect(
+          onCurveOutsideSubgroup(group(v), b),
+          cls == 'subgroup',
+          reason: label,
+        );
+      }
+      // A point re-encodes to its bytes.
+      if (s(v, 'go') == 'point') {
+        final encoded = group(v) == BlsGroup.g1
+            ? g1(s(v, 'hex')).toBytes()
+            : g2(s(v, 'hex')).toBytes();
+        expect(toHex(encoded), s(v, 'hex'), reason: label);
+      }
+    }
+    expect(classes['G1 subgroup'], greaterThanOrEqualTo(16));
+    expect(classes['G2 subgroup'], greaterThanOrEqualTo(8));
+    expect(classes['G1 curve'], greaterThanOrEqualTo(8));
+    expect(classes['G2 curve'], greaterThanOrEqualTo(4));
+  });
+
+  test('adds as kilic', () {
+    for (final v in section('add')) {
+      final label = s(v, 'label');
+      if (group(v) == BlsGroup.g1) {
+        final sum = g1(s(v, 'a')) + g1(s(v, 'b'));
+        expect(toHex(sum.toBytes()), s(v, 'sum'), reason: label);
+        // The mixed addition agrees.
+        final b = g1(s(v, 'b')).toAffine();
+        if (b != null) {
+          expect(
+            toHex(g1(s(v, 'a')).addAffine(b.$1, b.$2).toBytes()),
+            s(v, 'sum'),
+            reason: label,
+          );
+        }
+      } else {
+        final sum = g2(s(v, 'a')) + g2(s(v, 'b'));
+        expect(toHex(sum.toBytes()), s(v, 'sum'), reason: label);
+        final b = g2(s(v, 'b')).toAffine();
+        if (b != null) {
+          expect(
+            toHex(g2(s(v, 'a')).addAffine(b.$1, b.$2).toBytes()),
+            s(v, 'sum'),
+            reason: label,
+          );
+        }
+      }
+    }
+  });
+
+  test('multiplies as kilic, also by 0, r and scalars above r', () {
+    for (final v in section('multiply')) {
+      final k = BigInt.parse(s(v, 'scalar'), radix: 16);
+      final product = group(v) == BlsGroup.g1
+          ? g1(s(v, 'point')).multiply(k).toBytes()
+          : g2(s(v, 'point')).multiply(k).toBytes();
+      expect(toHex(product), s(v, 'product'), reason: s(v, 'label'));
+    }
+  });
+
+  test('pairs as kilic, serialized as kyber-bls12381 marshals GT', () {
+    for (final v in section('pairing')) {
+      final gt = pairing(g1(s(v, 'g1')), g2(s(v, 'g2')));
+      expect(toHex(gt.toBytes()), s(v, 'gt'), reason: s(v, 'label'));
+    }
+  });
+
+  test('hashes to G1 as kilic, for the DST of Quicknet and of RFC 9380', () {
+    final dsts = <String>{};
+    for (final v in section('hash_to_g1')) {
+      final p = hashToG1(fromHex(s(v, 'msg')), s(v, 'dst'));
+      final label = s(v, 'label');
+      expect(toHex(p.toBytes()), s(v, 'point'), reason: label);
+      final (x, y) = p.toAffine()!;
+      expect(
+        [toHex(x.toBytes()), toHex(y.toBytes())],
+        [s(v, 'x'), s(v, 'y')],
+        reason: label,
+      );
+      dsts.add(s(v, 'dst'));
+    }
+    expect(dsts, {
+      quicknetDst,
+      'QUUX-V01-CS02-with-BLS12381G1_XMD:SHA-256_SSWU_RO_',
+    });
+  });
+
+  test('maps one element to G1 as kilic, the exceptional ones included', () {
+    for (final v in section('map_to_g1')) {
+      final u = Fp.fromBytes(fromHex(s(v, 'u')))!;
+      expect(toHex(mapToG1(u).toBytes()), s(v, 'point'), reason: s(v, 'label'));
+    }
+  });
+
+  test('verifies BLS signatures on G1 as kyber sign/bls', () {
+    for (final v in section('signatures')) {
+      final key = g2(s(v, 'public_key'));
+      final sig = G1Point.decode(fromHex(s(v, 'signature')));
+      final ok =
+          sig != null &&
+          !sig.isInfinity &&
+          pairingCheck([
+            (hashToG1(fromHex(s(v, 'msg')), quicknetDst), key),
+            (-sig, G2Point.generator),
+          ]);
+      expect(ok, v['go'], reason: s(v, 'label'));
+    }
+  });
+}
+
+// Whether the compressed [b] has an x of a point of the curve, which then
+// lies outside the subgroup: what kilic reports as "not on correct
+// subgroup" rather than "not on curve".
+bool onCurveOutsideSubgroup(BlsGroup group, List<int> b) {
+  final raw = [b[0] & 0x1f, ...b.sublist(1)];
+  if (group == BlsGroup.g1) {
+    final x = Fp.fromBytes(raw)!;
+    final y = (x.square() * x + G1Point.b).sqrt();
+    return y != null && !G1Point.affine(x, y).isInSubgroup;
+  }
+  final x = Fp2(
+    Fp.fromBytes(raw.sublist(48))!,
+    Fp.fromBytes(raw.sublist(0, 48))!,
+  );
+  final y = (x.square() * x + G2Point.b).sqrt();
+  return y != null && !G2Point.affine(x, y).isInSubgroup;
+}

tool/bls12381_go_vectors.go

@@ -0,0 +1,546 @@
+//go:build ignore
+
+// Prints test/vectors/bls12381_vectors.json: the Go reference values for the
+// BLS12-381 code of lib/src/bls12381_*.dart, from the libraries that drand
+// and tlock use for Quicknet: github.com/kilic/bls12-381 and
+// github.com/drand/kyber-bls12381 over it.
+//
+//   - points: the frozen edge-case encodings of datekeys-ts
+//     (src/lib/dkc/testing/bls12381-vectors.json; valid points, sign-bit
+//     flips, identity encodings with stray flags or payload, missing
+//     compression flag, wrong lengths, uncompressed forms, x + p, points on
+//     the curve outside the subgroup, cofactor torsion), each with the
+//     verdict of the reference recomputed here: the KeyGroup of the drand
+//     crypto schemes (G1 for pedersen-bls-unchained, G2 for
+//     bls-unchained-g1-rfc9380), whose UnmarshalBinary is FromCompressed of
+//     kilic, and Equal(Null()) for the identity, as profile.Validate uses it.
+//   - decode: encodings drawn from a fixed seed: multiples of the generators,
+//     their negations, one flipped bit, random x with the compression flag,
+//     and random x of points on the curve outside the subgroup; each with the
+//     verdict, and for an invalid one the class of the error of kilic
+//     ("format", "curve" or "subgroup").
+//   - add, multiply: sums and multiples of points drawn from the seed, the
+//     special cases included (P + P, P + (-P), the point at infinity, the
+//     scalars 0, 1, r - 1, r, r + 1 and 2^256 - 1), computed by kilic.
+//   - pairing: e(P, Q) for points drawn from the seed, serialized by
+//     kyber-bls12381 (the order of kilic: c1 before c0 at every level).
+//   - hash_to_g1: HashToCurve of kilic for the DST of Quicknet and of tlock
+//     (BLS_SIG_BLS12381G1_XMD:SHA-256_SSWU_RO_NUL_) on messages drawn from
+//     the seed and on the identities of some rounds, and for the DST of the
+//     test vectors of RFC 9380, appendix J.9.1, on their messages, with the
+//     affine coordinates.
+//   - map_to_g1: MapToCurve of kilic (the simplified SWU map, the isogeny and
+//     the clearing of the cofactor) on elements u drawn from the seed and on
+//     the exceptional ones, u = 0 and Z·u² = -1, which no hash reaches.
+//   - signatures: BLS signatures on G1 of kyber's sign/bls (NewSchemeOnG1,
+//     the scheme of Quicknet) under keys drawn from the seed, and the verdict
+//     of its Verify on each and on edited copies.
+//
+// The seed is fixed: the output is the same on every run. Run it from a
+// scratch module that requires the reference implementation at spec-v0.11
+// (a module scratch whose go.mod has require g.activething.com/go/DateKeys
+// v0.0.0 and replace g.activething.com/go/DateKeys => an export of that tag,
+// with its go.sum, and GOFLAGS=-mod=mod), copied into it, passing the frozen
+// file of datekeys-ts:
+//
+//	go run bls12381_go_vectors.go path/to/bls12381-vectors.json > bls12381_vectors.json
+package main
+
+import (
+	"bytes"
+	"crypto/sha256"
+	"encoding/binary"
+	"encoding/hex"
+	"encoding/json"
+	"fmt"
+	"math/big"
+	"os"
+	"runtime/debug"
+	"sort"
+	"strings"
+
+	"github.com/drand/drand/v2/common"
+	"github.com/drand/drand/v2/crypto"
+	"github.com/drand/kyber"
+	bls "github.com/drand/kyber-bls12381"
+	signbls "github.com/drand/kyber/sign/bls"
+	bls12381 "github.com/kilic/bls12-381"
+)
+
+const quicknetDST = "BLS_SIG_BLS12381G1_XMD:SHA-256_SSWU_RO_NUL_"
+
+var (
+	p, _     = new(big.Int).SetString("1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab", 16)
+	order, _ = new(big.Int).SetString("73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001", 16)
+	g1       = bls12381.NewG1()
+	g2       = bls12381.NewG2()
+)
+
+func must[T any](v T, err error) T {
+	if err != nil {
+		panic(err)
+	}
+	return v
+}
+
+// stream is a deterministic byte source: SHA-256 of the label, a counter
+// and a block index.
+type stream struct {
+	label string
+	n     uint32
+}
+
+func (s *stream) bytes(n int) []byte {
+	s.n++
+	var out []byte
+	for i := uint32(0); len(out) < n; i++ {
+		h := sha256.Sum256(binary.BigEndian.AppendUint32(binary.BigEndian.AppendUint32([]byte("DateKeys BLS12-381 vectors "+s.label), s.n), i))
+		out = append(out, h[:]...)
+	}
+	return out[:n]
+}
+
+// scalar is a nonzero scalar below r.
+func (s *stream) scalar() *big.Int {
+	k := new(big.Int).SetBytes(s.bytes(64))
+	k.Mod(k, new(big.Int).Sub(order, big.NewInt(1)))
+	return k.Add(k, big.NewInt(1))
+}
+
+func (s *stream) fp() *big.Int {
+	x := new(big.Int).SetBytes(s.bytes(64))
+	return x.Mod(x, p)
+}
+
+func fp48(x *big.Int) []byte { return x.FillBytes(make([]byte, 48)) }
+
+func hx(b []byte) string { return hex.EncodeToString(b) }
+
+func g1Mul(k *big.Int) *bls12381.PointG1 {
+	return g1.MulScalarBig(g1.New(), g1.One(), k)
+}
+
+func g2Mul(k *big.Int) *bls12381.PointG2 {
+	return g2.MulScalarBig(g2.New(), g2.One(), k)
+}
+
+func g1Enc(q *bls12381.PointG1) []byte { return g1.ToCompressed(g1.New().Set(q)) }
+
+func g2Enc(q *bls12381.PointG2) []byte { return g2.ToCompressed(g2.New().Set(q)) }
+
+// errorClass names the check of FromCompressed of kilic that failed.
+func errorClass(err error) string {
+	switch {
+	case strings.Contains(err.Error(), "not on curve"):
+		return "curve"
+	case strings.Contains(err.Error(), "correct subgroup"):
+		return "subgroup"
+	default:
+		return "format"
+	}
+}
+
+type decodeVector struct {
+	Label string `json:"label"`
+	Group string `json:"group"`
+	Hex   string `json:"hex"`
+	Go    string `json:"go"`
+	Class string `json:"class,omitempty"`
+}
+
+func verdict(group string, b []byte) decodeVector {
+	v := decodeVector{Group: group, Hex: hx(b)}
+	var err error
+	var zero bool
+	if group == "G1" {
+		var q *bls12381.PointG1
+		if q, err = g1.FromCompressed(b); err == nil {
+			zero = g1.IsZero(q)
+		}
+	} else {
+		var q *bls12381.PointG2
+		if q, err = g2.FromCompressed(b); err == nil {
+			zero = g2.IsZero(q)
+		}
+	}
+	switch {
+	case err != nil:
+		v.Go, v.Class = "invalid", errorClass(err)
+	case zero:
+		v.Go = "identity"
+	default:
+		v.Go = "point"
+	}
+	return v
+}
+
+// schemeVerdict is the verdict of the KeyGroup of the drand scheme, as
+// scripts/bls12381-go-verdicts.go of datekeys-ts computes it.
+func schemeVerdict(group string, b []byte) string {
+	id := crypto.UnchainedSchemeID
+	if group == "G2" {
+		id = crypto.SigsOnG1ID
+	}
+	s := must(crypto.SchemeFromName(id))
+	k := s.KeyGroup.Point()
+	switch {
+	case k.UnmarshalBinary(b) != nil:
+		return "invalid"
+	case k.Equal(k.Null()):
+		return "identity"
+	default:
+		return "point"
+	}
+}
+
+type addVector struct {
+	Label string `json:"label"`
+	Group string `json:"group"`
+	A     string `json:"a"`
+	B     string `json:"b"`
+	Sum   string `json:"sum"`
+}
+
+type mulVector struct {
+	Label   string `json:"label"`
+	Group   string `json:"group"`
+	Point   string `json:"point"`
+	Scalar  string `json:"scalar"`
+	Product string `json:"product"`
+}
+
+type pairingVector struct {
+	Label string `json:"label"`
+	G1    string `json:"g1"`
+	G2    string `json:"g2"`
+	GT    string `json:"gt"`
+}
+
+type hashVector struct {
+	Label string `json:"label"`
+	DST   string `json:"dst"`
+	Msg   string `json:"msg"`
+	Point string `json:"point"`
+	X     string `json:"x"`
+	Y     string `json:"y"`
+}
+
+type mapVector struct {
+	Label string `json:"label"`
+	U     string `json:"u"`
+	Point string `json:"point"`
+}
+
+type signatureVector struct {
+	Label     string `json:"label"`
+	PublicKey string `json:"public_key"`
+	Msg       string `json:"msg"`
+	Signature string `json:"signature"`
+	Go        bool   `json:"go"`
+}
+
+func main() {
+	out := struct {
+		Description string            `json:"description"`
+		Generator   string            `json:"generator"`
+		Libraries   string            `json:"libraries"`
+		PointsFrom  string            `json:"points_from"`
+		Points      []decodeVector    `json:"points"`
+		Decode      []decodeVector    `json:"decode"`
+		Add         []addVector       `json:"add"`
+		Multiply    []mulVector       `json:"multiply"`
+		Pairing     []pairingVector   `json:"pairing"`
+		HashToG1    []hashVector      `json:"hash_to_g1"`
+		MapToG1     []mapVector       `json:"map_to_g1"`
+		Signatures  []signatureVector `json:"signatures"`
+	}{
+		Description: "Go reference values for the BLS12-381 code of lib/src/bls12381_*.dart: decoding verdicts, " +
+			"sums, multiples, pairings, hashes to G1 and BLS signatures on G1; see tool/bls12381_go_vectors.go.",
+		Generator: "tool/bls12381_go_vectors.go",
+		Libraries: libraries(),
+		PointsFrom: "the encodings of src/lib/dkc/testing/bls12381-vectors.json of datekeys-ts at 289fe71, " +
+			"with the verdicts recomputed by this generator",
+	}
+
+	// The frozen edge cases of datekeys-ts.
+	var frozen struct {
+		Vectors []struct{ Label, Group, Hex, Go string }
+	}
+	if err := json.Unmarshal(must(os.ReadFile(os.Args[1])), &frozen); err != nil {
+		panic(err)
+	}
+	for _, f := range frozen.Vectors {
+		b := must(hex.DecodeString(f.Hex))
+		v := verdict(f.Group, b)
+		if (len(b) == 48 && f.Group == "G1") || (len(b) == 96 && f.Group == "G2") {
+			if s := schemeVerdict(f.Group, b); s != v.Go {
+				panic(f.Label + ": kilic and the drand scheme disagree")
+			}
+		} else {
+			v.Go, v.Class = schemeVerdict(f.Group, b), "format"
+		}
+		if v.Go != f.Go {
+			panic(f.Label + ": the verdict of datekeys-ts is not the one of Go")
+		}
+		v.Label = f.Label
+		out.Points = append(out.Points, v)
+	}
+
+	// Decoding.
+	s := &stream{label: "decode"}
+	for i := 0; i < 24; i++ {
+		b := g1Enc(g1Mul(s.scalar()))
+		neg := bytes.Clone(b)
+		neg[0] ^= 0x20
+		flip := bytes.Clone(b)
+		bit := 3 + int(binary.BigEndian.Uint16(s.bytes(2)))%(48*8-3)
+		flip[bit/8] ^= 0x80 >> (bit % 8)
+		for j, c := range [][]byte{b, neg, flip} {
+			v := verdict("G1", c)
+			v.Label = fmt.Sprintf("G1 multiple %d, %s", i, []string{"as encoded", "negated", fmt.Sprintf("bit %d flipped", bit)}[j])
+			out.Decode = append(out.Decode, v)
+		}
+	}
+	for i := 0; i < 8; i++ {
+		b := g2Enc(g2Mul(s.scalar()))
+		neg := bytes.Clone(b)
+		neg[0] ^= 0x20
+		flip := bytes.Clone(b)
+		bit := 3 + int(binary.BigEndian.Uint16(s.bytes(2)))%(96*8-3)
+		flip[bit/8] ^= 0x80 >> (bit % 8)
+		for j, c := range [][]byte{b, neg, flip} {
+			v := verdict("G2", c)
+			v.Label = fmt.Sprintf("G2 multiple %d, %s", i, []string{"as encoded", "negated", fmt.Sprintf("bit %d flipped", bit)}[j])
+			out.Decode = append(out.Decode, v)
+		}
+	}
+	// Random x, with the compression flag and either sign: on the curve or
+	// not, and then outside the subgroup.
+	classes := map[string]int{}
+	for i := 0; classes["G1 subgroup"] < 16 || classes["G1 curve"] < 8; i++ {
+		b := fp48(s.fp())
+		b[0] |= 0x80 | s.bytes(1)[0]&0x20
+		v := verdict("G1", b)
+		key := "G1 " + v.Class
+		if classes[key] >= 16 {
+			continue
+		}
+		classes[key]++
+		v.Label = fmt.Sprintf("G1 random x %d", i)
+		out.Decode = append(out.Decode, v)
+	}
+	for i := 0; classes["G2 subgroup"] < 8 || classes["G2 curve"] < 4; i++ {
+		b := append(fp48(s.fp()), fp48(s.fp())...)
+		b[0] |= 0x80 | s.bytes(1)[0]&0x20
+		v := verdict("G2", b)
+		key := "G2 " + v.Class
+		if classes[key] >= 8 {
+			continue
+		}
+		classes[key]++
+		v.Label = fmt.Sprintf("G2 random x %d", i)
+		out.Decode = append(out.Decode, v)
+	}
+
+	// Sums.
+	s = &stream{label: "add"}
+	inf1, inf2 := g1.Zero(), g2.Zero()
+	for i := 0; i < 16; i++ {
+		a, b := g1Mul(s.scalar()), g1Mul(s.scalar())
+		out.Add = append(out.Add, addVector{fmt.Sprintf("G1 sum %d", i), "G1", hx(g1Enc(a)), hx(g1Enc(b)), hx(g1Enc(g1.Add(g1.New(), a, b)))})
+	}
+	a1 := g1Mul(s.scalar())
+	for _, c := range []struct {
+		label string
+		a, b  *bls12381.PointG1
+	}{
+		{"G1 P + P", a1, a1},
+		{"G1 P + (-P)", a1, g1.Neg(g1.New(), a1)},
+		{"G1 P + infinity", a1, inf1},
+		{"G1 infinity + P", inf1, a1},
+		{"G1 infinity + infinity", inf1, inf1},
+	} {
+		out.Add = append(out.Add, addVector{c.label, "G1", hx(g1Enc(c.a)), hx(g1Enc(c.b)), hx(g1Enc(g1.Add(g1.New(), c.a, c.b)))})
+	}
+	for i := 0; i < 8; i++ {
+		a, b := g2Mul(s.scalar()), g2Mul(s.scalar())
+		out.Add = append(out.Add, addVector{fmt.Sprintf("G2 sum %d", i), "G2", hx(g2Enc(a)), hx(g2Enc(b)), hx(g2Enc(g2.Add(g2.New(), a, b)))})
+	}
+	a2 := g2Mul(s.scalar())
+	for _, c := range []struct {
+		label string
+		a, b  *bls12381.PointG2
+	}{
+		{"G2 P + P", a2, a2},
+		{"G2 P + (-P)", a2, g2.Neg(g2.New(), a2)},
+		{"G2 P + infinity", a2, inf2},
+		{"G2 infinity + P", inf2, a2},
+	} {
+		out.Add = append(out.Add, addVector{c.label, "G2", hx(g2Enc(c.a)), hx(g2Enc(c.b)), hx(g2Enc(g2.Add(g2.New(), c.a, c.b)))})
+	}
+
+	// Multiples.
+	s = &stream{label: "multiply"}
+	max256 := new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 256), big.NewInt(1))
+	special := []struct {
+		label string
+		k     *big.Int
+	}{
+		{"0", big.NewInt(0)},
+		{"1", big.NewInt(1)},
+		{"2", big.NewInt(2)},
+		{"r - 1", new(big.Int).Sub(order, big.NewInt(1))},
+		{"r", new(big.Int).Set(order)},
+		{"r + 1", new(big.Int).Add(order, big.NewInt(1))},
+		{"2^256 - 1", max256},
+	}
+	for i := 0; i < 12; i++ {
+		q, k := g1Mul(s.scalar()), new(big.Int).SetBytes(s.bytes(32))
+		out.Multiply = append(out.Multiply, mulVector{fmt.Sprintf("G1 multiple %d", i), "G1", hx(g1Enc(q)), k.Text(16), hx(g1Enc(g1.MulScalarBig(g1.New(), q, k)))})
+	}
+	q1 := g1Mul(s.scalar())
+	for _, c := range special {
+		out.Multiply = append(out.Multiply, mulVector{"G1 by " + c.label, "G1", hx(g1Enc(q1)), c.k.Text(16), hx(g1Enc(g1.MulScalarBig(g1.New(), q1, c.k)))})
+	}
+	for i := 0; i < 6; i++ {
+		q, k := g2Mul(s.scalar()), new(big.Int).SetBytes(s.bytes(32))
+		out.Multiply = append(out.Multiply, mulVector{fmt.Sprintf("G2 multiple %d", i), "G2", hx(g2Enc(q)), k.Text(16), hx(g2Enc(g2.MulScalarBig(g2.New(), q, k)))})
+	}
+	q2 := g2Mul(s.scalar())
+	for _, c := range special {
+		out.Multiply = append(out.Multiply, mulVector{"G2 by " + c.label, "G2", hx(g2Enc(q2)), c.k.Text(16), hx(g2Enc(g2.MulScalarBig(g2.New(), q2, c.k)))})
+	}
+
+	// Pairings, through kyber-bls12381.
+	s = &stream{label: "pairing"}
+	suite := bls.NewBLS12381Suite()
+	for i := 0; i < 6; i++ {
+		a := suite.G1().Point().Mul(scalarOf(s.scalar()), nil)
+		b := suite.G2().Point().Mul(scalarOf(s.scalar()), nil)
+		out.Pairing = append(out.Pairing, pairingVector{fmt.Sprintf("e(P, Q) %d", i), marshal(a), marshal(b), marshal(suite.Pair(a, b))})
+	}
+	{
+		a := suite.G1().Point().Mul(scalarOf(s.scalar()), nil)
+		b := suite.G2().Point().Mul(scalarOf(s.scalar()), nil)
+		out.Pairing = append(out.Pairing,
+			pairingVector{"e(infinity, Q)", marshal(suite.G1().Point().Null()), marshal(b), marshal(suite.Pair(suite.G1().Point().Null(), b))},
+			pairingVector{"e(P, infinity)", marshal(a), marshal(suite.G2().Point().Null()), marshal(suite.Pair(a, suite.G2().Point().Null()))})
+	}
+
+	// Hashes to G1.
+	s = &stream{label: "hash"}
+	hashOne := func(label, dst string, msg []byte) hashVector {
+		q := must(g1.HashToCurve(msg, []byte(dst)))
+		u := g1.ToUncompressed(g1.New().Set(q))
+		return hashVector{label, dst, hx(msg), hx(g1Enc(q)), hx(u[:48]), hx(u[48:])}
+	}
+	for i := 0; i < 24; i++ {
+		n := int(s.bytes(1)[0]) % 80
+		out.HashToG1 = append(out.HashToG1, hashOne(fmt.Sprintf("a message of %d bytes", n), quicknetDST, s.bytes(n)))
+	}
+	sch := must(crypto.SchemeFromName(crypto.SigsOnG1ID))
+	for _, round := range []uint64{1, 1000, 1001, 1004, 2000, 83903165811, 1<<53 - 1} {
+		id := sch.DigestBeacon(&common.Beacon{Round: round})
+		out.HashToG1 = append(out.HashToG1, hashOne(fmt.Sprintf("the identity of round %d", round), quicknetDST, id))
+	}
+	const rfcDST = "QUUX-V01-CS02-with-BLS12381G1_XMD:SHA-256_SSWU_RO_"
+	for _, m := range []string{"", "abc", "abcdef0123456789", "q128_" + strings.Repeat("q", 128), "a512_" + strings.Repeat("a", 512)} {
+		out.HashToG1 = append(out.HashToG1, hashOne(fmt.Sprintf("RFC 9380 J.9.1, msg %q", shorten(m)), rfcDST, []byte(m)))
+	}
+
+	// The map to G1 of one element.
+	s = &stream{label: "map"}
+	z := big.NewInt(11)
+	// Z·u² = -1: u² = -1/Z, a square since -Z is one.
+	minusInvZ := new(big.Int).Sub(p, new(big.Int).ModInverse(z, p))
+	root := new(big.Int).ModSqrt(minusInvZ, p)
+	if root == nil {
+		panic("-1/Z is not a square")
+	}
+	us := []struct {
+		label string
+		u     *big.Int
+	}{
+		{"u = 0, exceptional", big.NewInt(0)},
+		{"Z·u² = -1, exceptional", root},
+		{"Z·u² = -1, the other root, exceptional", new(big.Int).Sub(p, root)},
+		{"u = 1", big.NewInt(1)},
+		{"u = p - 1", new(big.Int).Sub(p, big.NewInt(1))},
+	}
+	for i := 0; i < 8; i++ {
+		us = append(us, struct {
+			label string
+			u     *big.Int
+		}{fmt.Sprintf("u drawn from the seed %d", i), s.fp()})
+	}
+	for _, c := range us {
+		q := must(g1.MapToCurve(fp48(c.u)))
+		out.MapToG1 = append(out.MapToG1, mapVector{c.label, hx(fp48(c.u)), hx(g1Enc(q))})
+	}
+
+	// BLS signatures on G1, the scheme of Quicknet.
+	s = &stream{label: "signatures"}
+	scheme := signbls.NewSchemeOnG1(suite)
+	for i := 0; i < 6; i++ {
+		sk := scalarOf(s.scalar())
+		pk := suite.G2().Point().Mul(sk, nil)
+		msg := s.bytes(32)
+		sig := must(scheme.Sign(sk, msg))
+		other := s.bytes(32)
+		for _, c := range []struct {
+			label    string
+			msg, sig []byte
+		}{
+			{fmt.Sprintf("signature %d", i), msg, sig},
+			{fmt.Sprintf("signature %d of another message", i), other, sig},
+			{fmt.Sprintf("signature %d negated", i), msg, negate(sig)},
+		} {
+			out.Signatures = append(out.Signatures, signatureVector{c.label, marshal(pk), hx(c.msg), hx(c.sig), scheme.Verify(pk, c.msg, c.sig) == nil})
+		}
+	}
+
+	enc := json.NewEncoder(os.Stdout)
+	enc.SetIndent("", " ")
+	enc.SetEscapeHTML(false)
+	if err := enc.Encode(out); err != nil {
+		panic(err)
+	}
+}
+
+func scalarOf(k *big.Int) kyber.Scalar {
+	return bls.NewKyberScalar().SetBytes(k.FillBytes(make([]byte, 32)))
+}
+
+func marshal(m interface{ MarshalBinary() ([]byte, error) }) string {
+	return hex.EncodeToString(must(m.MarshalBinary()))
+}
+
+func negate(sig []byte) []byte {
+	c := bytes.Clone(sig)
+	c[0] ^= 0x20
+	return c
+}
+
+func shorten(m string) string {
+	if len(m) > 20 {
+		return m[:20] + "…"
+	}
+	return m
+}
+
+// libraries names the versions of the libraries this program ran with.
+func libraries() string {
+	info, ok := debug.ReadBuildInfo()
+	if !ok {
+		panic("no build info")
+	}
+	var out []string
+	for _, d := range info.Deps {
+		switch d.Path {
+		case "github.com/kilic/bls12-381", "github.com/drand/kyber", "github.com/drand/kyber-bls12381", "github.com/drand/drand/v2":
+			out = append(out, d.Path+" "+d.Version)
+		}
+	}
+	sort.Strings(out)
+	return strings.Join(out, ", ")
+}