Performance: field-op counts + batched group law + smalljac comparison

README.md

@@ -108,6 +108,8 @@ Ramified model:
 | deg2 + deg2 (char2) | GF(2^8) | ~185 ns |
 | deg2 + deg2 (char2) | GF(2^16) | ~600 ns |
 | 2*deg2 (not_char2) | F_65521 | ~194 ns |
+| deg2 + deg2 (not_char2) | F_p (56-bit) | ~741 ns |
+| 2*deg2 (not_char2) | F_p (56-bit) | ~682 ns |
 
 Split model (degree-2 balanced divisors, negative basis):
 
@@ -125,6 +127,63 @@ Run benchmarks with:
 cargo bench
 ```
 
+### Field-operation counts
+
+Wall-clock comparisons across implementations are confounded by field size,
+operation definition, and language, so the explicit-formula literature compares
+**field-operation counts** instead (field-size independent). Counts below are
+for the generic-branch degree-2 `add`/`double` (negative basis), measured by
+running the actual formulas over an instrumented field
+(`cargo test --release g2::split::op_counts -- --nocapture`):
+
+| Operation | M (mul) | S (sqr) | I (inv) | A (add/sub/dbl) |
+|-----------|--------:|--------:|--------:|----------------:|
+| not_char2 — add    | 26 | 2 | 1 | 37 |
+| not_char2 — double | 31 | 3 | 1 | 38 |
+| arbitrary — add    | 30 | 1 | 1 | 36 |
+| arbitrary — double | 38 | 2 | 1 | 44 |
+| char2 — add        | 27 | 1 | 1 | 34 |
+| char2 — double     | 29 | 2 | 1 | 31 |
+
+Each group operation uses exactly **one** field inversion (affine formulas).
+These can be compared directly to the per-formula counts in Lange, Erickson–
+Jacobson–Stein (real genus 2), and Costello–Lauter.
+
+### Wall-clock comparison with smalljac
+
+[smalljac](https://math.mit.edu/~drew/smalljac.html) (Andrew Sutherland) is a
+highly optimized C library whose `hecurve_g2_compose` / `hecurve_g2_square`
+implement the genus-2 **imaginary** (ramified, `deg f = 5`) group law — the same
+model as this crate's `g2::ramified::not_char2`. Built and timed on the same
+machine (smalljac v4.1.3 + ff_poly v1.2.7, ported to arm64), exercising the
+affine path (`ctx = NULL`, one field inversion per op — matching this crate's
+affine formulas):
+
+| field | operation | this crate (ramified nch2) | smalljac |
+|-------|-----------|---------------------------:|---------:|
+| p = 65521 (16-bit) | add / double | 149 / 194 ns | 767 / 873 ns |
+| 56-bit prime (matched width) | add / double | 741 / 682 ns | 1470 / 1619 ns |
+
+**These numbers need context — cross-implementation wall-clock is confounded:**
+
+- **Field width dominates.** ff_poly is compiled for ≤57-bit primes and always
+  does 64-bit-wide Montgomery arithmetic, so smalljac barely changes from 16-bit
+  to 56-bit (767 → 1470 ns); this crate's `PrimeField` uses `u128`-multiply +
+  hardware modulo, much faster at 16-bit but scaling up with the modulus. **Only
+  the matched 56-bit row is a fair wall-clock comparison.**
+- **Batched inversion is disabled.** smalljac's real strength for point counting
+  is amortizing one inversion across many group ops (Montgomery's trick, via its
+  `ctx` state machine). Forcing `ctx = NULL` measures its un-batched affine path
+  — the right comparison for a *single* op, but not how smalljac runs in anger.
+- **Specialized vs general.** This crate's `add`/`double` are degree-2-specialized
+  explicit formulas (≈26 M, 1 I); `hecurve_g2_compose` is a general composition
+  routine that also handles the degenerate-degree cases.
+
+So the field-operation counts above remain the cleaner, field-size-independent
+comparison; the matched-width wall-clock merely confirms the specialized
+explicit formulas are competitive with a mature C implementation. The harness
+and build notes are in [`benches/smalljac-compare/`](benches/smalljac-compare/).
+
 ## Testing
 
 Run all tests:

benches/benchmarks.rs

@@ -477,10 +477,13 @@ fn not_char2_benchmarks(c: &mut Criterion) {
     bench_not_char2_deg2_add::<PrimeField<7>>(c, "F7");
     bench_not_char2_deg2_add::<PrimeField<8191>>(c, "F8191");
     bench_not_char2_deg2_add::<PrimeField<65521>>(c, "F65521");
+    // ~56-bit prime: matched-width comparison vs smalljac (built for 57-bit primes).
+    bench_not_char2_deg2_add::<PrimeField<72057594037927931>>(c, "Fp56");
 
     bench_not_char2_deg2_dbl::<PrimeField<7>>(c, "F7");
     bench_not_char2_deg2_dbl::<PrimeField<8191>>(c, "F8191");
     bench_not_char2_deg2_dbl::<PrimeField<65521>>(c, "F65521");
+    bench_not_char2_deg2_dbl::<PrimeField<72057594037927931>>(c, "Fp56");
 }
 
 fn arbitrary_benchmarks(c: &mut Criterion) {

benches/smalljac-compare/README.md

@@ -0,0 +1,61 @@
+# smalljac wall-clock comparison
+
+[`bench.c`](bench.c) times smalljac's genus-2 imaginary/ramified group law
+(`hecurve_g2_compose` = add, `hecurve_g2_square` = double) on the affine path
+(`ctx = NULL`, one field inversion per op), so it can be put next to this
+crate's `g2::ramified::not_char2` `cargo bench` numbers. See the
+"Wall-clock comparison with smalljac" section of the top-level
+[`README.md`](../../README.md) for the results and the caveats.
+
+smalljac and ff_poly are **not vendored** here (they are GPL; this crate is
+MIT). Download them yourself:
+
+- smalljac v4.1.3 and ff_poly v1.2.7: <https://math.mit.edu/~drew/smalljac.html>
+- GMP (e.g. `brew install gmp` on macOS)
+
+## Building on x86-64 Linux
+
+ff_poly's hand-written x86-64 assembly works as-is there:
+
+```sh
+# build libff_poly.a per its Makefile, then:
+cc -O3 -std=gnu99 -DSMALLJAC_GENUS=2 \
+   -I/path/to/ff_poly -I/path/to/smalljac bench.c \
+   hecurve.c hecurve1.c hecurve2.c mpzutil.c prime.c \
+   /path/to/libff_poly.a -lgmp -o sjbench
+./sjbench
+```
+
+## Building on arm64 (Apple Silicon)
+
+ff_poly and smalljac use x86-64 inline asm (`mulq`, `bsrq`, `bsfq`, …). Replace
+it with portable `__uint128_t` / `__builtin_*` equivalents (semantics identical;
+clang lowers the 128-bit ops to native `mul`/`umulh`/`clz`/`ctz` on aarch64).
+
+**1. Replace `ff_poly/asm.h`** with these portable macros:
+
+```c
+#define _asm_div_q_q(q,r,x,y)        do { __uint128_t _n=((__uint128_t)(unsigned long)(r)<<64)|(unsigned long)(x); (q)=(unsigned long)(_n/(unsigned long)(y)); (r)=(unsigned long)(_n%(unsigned long)(y)); } while(0)
+#define _asm_mult_1_1(z1,z0,x0,y0)   do { __uint128_t _p=(__uint128_t)(unsigned long)(x0)*(unsigned long)(y0); (z0)=(unsigned long)_p; (z1)=(unsigned long)(_p>>64); } while(0)
+#define _asm_mult_2_2_1(z1,z0,x1,x0,y0) do { __uint128_t _p=(__uint128_t)(unsigned long)(x0)*(unsigned long)(y0); (z0)=(unsigned long)_p; (z1)=(unsigned long)(_p>>64)+(unsigned long)(x1)*(unsigned long)(y0); } while(0)
+#define _asm_addto_2_2(z1,z0,x1,x0)  do { __uint128_t _s=(((__uint128_t)(unsigned long)(z1)<<64)|(unsigned long)(z0))+(((__uint128_t)(unsigned long)(x1)<<64)|(unsigned long)(x0)); (z0)=(unsigned long)_s; (z1)=(unsigned long)(_s>>64); } while(0)
+#define _asm_addto_2_1(z1,z0,x0)     do { __uint128_t _s=(((__uint128_t)(unsigned long)(z1)<<64)|(unsigned long)(z0))+(unsigned long)(x0); (z0)=(unsigned long)_s; (z1)=(unsigned long)(_s>>64); } while(0)
+#define _asm_addto_3_3(z2,z1,z0,x2,x1,x0) do { __uint128_t _s=(__uint128_t)(unsigned long)(z0)+(unsigned long)(x0); (z0)=(unsigned long)_s; _s=(__uint128_t)(unsigned long)(z1)+(unsigned long)(x1)+(unsigned long)(_s>>64); (z1)=(unsigned long)_s; (z2)=(unsigned long)(z2)+(unsigned long)(x2)+(unsigned long)(_s>>64); } while(0)
+#define _asm_addto_3_2(z2,z1,z0,x1,x0)    do { __uint128_t _s=(__uint128_t)(unsigned long)(z0)+(unsigned long)(x0); (z0)=(unsigned long)_s; _s=(__uint128_t)(unsigned long)(z1)+(unsigned long)(x1)+(unsigned long)(_s>>64); (z1)=(unsigned long)_s; (z2)=(unsigned long)(z2)+(unsigned long)(_s>>64); } while(0)
+#define _asm_subfrom_2_2(z1,z0,x1,x0) do { __uint128_t _d=(((__uint128_t)(unsigned long)(z1)<<64)|(unsigned long)(z0))-(((__uint128_t)(unsigned long)(x1)<<64)|(unsigned long)(x0)); (z0)=(unsigned long)_d; (z1)=(unsigned long)(_d>>64); } while(0)
+#define _asm_inc_2(z1,z0)            do { if(++(z0)==0UL) (z1)++; } while(0)
+```
+
+(plus the `_asm_mult_3_2_1` / `_asm_mult_3_2_2` / `_asm_square_3_2` helpers,
+which are already written in terms of the macros above — keep them as-is.)
+
+**2. In both `cstd.h` files** (the ff_poly copy *and* the smalljac copy) replace
+the two bit-scan helpers:
+
+```c
+static inline unsigned long _asm_highbit (unsigned long x) { return 63UL - __builtin_clzl(x); }
+static inline unsigned long _asm_lowbit  (unsigned long x) { return __builtin_ctzl(x); }
+```
+
+Then build with the same `cc` line as above. ff_poly is then limited to
+`FF_BITS = 57`, which is why the matched-width comparison uses a 56-bit prime.

benches/smalljac-compare/bench.c

@@ -0,0 +1,90 @@
+/* Micro-benchmark for smalljac's genus-2 Jacobian group law
+ * (hecurve_g2_compose = ADD, hecurve_g2_square = DBL), to compare against the
+ * Rust `divisor-arithmetic` cargo-bench numbers for g2::ramified::not_char2.
+ *
+ * Imaginary/ramified genus-2 model: y^2 = f(x), deg f = 5 (monic).
+ * Times the affine path (ctx = NULL => one field inversion per op), which is
+ * the same affine, one-inversion model the Rust explicit formulas use.
+ *
+ * Build: see README.md in this directory (requires an external smalljac +
+ * ff_poly checkout; on arm64 the x86-64 inline asm must be replaced — the
+ * README lists the two portable patches).
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <gmp.h>
+#include "ff_poly.h"
+#include "mpzutil.h"
+#include "hecurve.h"
+#include "cstd.h"
+
+#define POOL 256
+
+static double now_ns(void) {
+    struct timespec ts;
+    clock_gettime(CLOCK_MONOTONIC, &ts);
+    return ts.tv_sec * 1e9 + ts.tv_nsec;
+}
+
+static void run(unsigned long p, const char *label) {
+    ff_setup_ui(p);
+
+    /* monic degree-5 f with random coefficients f[0..4], f[5]=1, f[6]=0 */
+    ff_t f[HECURVE_DEGREE + 1];
+    for (int i = 0; i < HECURVE_DEGREE + 1; i++) _ff_set_zero(f[i]);
+    for (int i = 0; i < 5; i++) _ff_random(f[i]);
+    _ff_set_one(f[5]);
+
+    /* pool of valid random weight-2 divisors */
+    static ff_t U[POOL][HECURVE_GENUS + 1], V[POOL][HECURVE_GENUS];
+    int n = 0;
+    for (int tries = 0; tries < POOL * 50 && n < POOL; tries++) {
+        ff_t u[HECURVE_GENUS + 1], v[HECURVE_GENUS];
+        hecurve_random(u, v, f);
+        if (_ff_zero(u[2]) || !_ff_one(u[2])) continue; /* keep monic deg-2 u */
+        if (!hecurve_verify(u, v, f)) continue;
+        for (int i = 0; i < HECURVE_GENUS + 1; i++) _ff_set(U[n][i], u[i]);
+        for (int i = 0; i < HECURVE_GENUS; i++) _ff_set(V[n][i], v[i]);
+        n++;
+    }
+    if (n < 8) { printf("%-12s p=%lu: only %d divisors, skipping\n", label, p, n); return; }
+
+    ff_t ru[HECURVE_GENUS + 1], rv[HECURVE_GENUS];
+    unsigned long sink = 0;
+    long iters = 2000000;
+
+    /* ---- ADD (compose distinct divisors) ---- */
+    for (int i = 0; i < n - 1; i++) hecurve_g2_compose(ru, rv, U[i], V[i], U[i+1], V[i+1], f, 0);
+    double t0 = now_ns();
+    for (long k = 0; k < iters; k++) {
+        int i = k % (n - 1);
+        hecurve_g2_compose(ru, rv, U[i], V[i], U[i+1], V[i+1], f, 0);
+        sink += _ff_get_ui(ru[0]);
+    }
+    double add_ns = (now_ns() - t0) / iters;
+
+    /* ---- DBL (square a divisor) ---- */
+    for (int i = 0; i < n; i++) hecurve_g2_square(ru, rv, U[i], V[i], f, 0);
+    t0 = now_ns();
+    for (long k = 0; k < iters; k++) {
+        int i = k % n;
+        hecurve_g2_square(ru, rv, U[i], V[i], f, 0);
+        sink += _ff_get_ui(ru[0]);
+    }
+    double dbl_ns = (now_ns() - t0) / iters;
+
+    printf("%-12s p=%-20lu add=%7.1f ns   double=%7.1f ns   (n=%d, sink=%lu)\n",
+           label, p, add_ns, dbl_ns, n, sink);
+}
+
+int main(void) {
+    run(8191UL, "13-bit");
+    run(65521UL, "16-bit");
+    run(1048583UL, "20-bit");
+    run(1000000007UL, "30-bit");
+    run((1UL << 31) - 1, "31-bit M31");
+    /* 56-bit prime: matched-width comparison vs Rust PrimeField<72057594037927931> */
+    run(72057594037927931UL, "56-bit Fp56");
+    return 0;
+}

src/g2/split/mod.rs

@@ -64,6 +64,8 @@ mod blackbox_tests;
 #[cfg(test)]
 mod char2_tests;
 #[cfg(test)]
+mod op_counts;
+#[cfg(test)]
 mod test_support;
 #[cfg(test)]
 mod wb_vectors;