diff --git a/README.md b/README.md index f2d1e2e..841ede2 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,11 @@ This crate provides highly optimized implementations of divisor addition and dou - All formulas are cross-checked against a generic Cantor reference implementation (`generic::split`) +- **Batched group law** (ramified `not_char2`): `add_batch` / `double_batch` + amortize the single field inversion across a whole batch of independent + operations via Montgomery's trick (`field::batch_invert`) — the same strategy + smalljac uses for generic-group order computations + - **Field Implementations**: - `PrimeField

` - Prime fields F_p for small primes - `BinaryExtField` - Binary extension fields GF(2^k) for k ≤ 24 @@ -108,8 +113,6 @@ 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): @@ -149,40 +152,58 @@ 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 +### Wall-clock comparison with smalljac (scalar and batched) [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/). +model as this crate's `g2::ramified::not_char2`. Both are compared two ways: + +- **scalar** — one field inversion per group operation (this crate's plain + `add`/`double`; smalljac with `ctx = NULL`); +- **batched** — one field inversion shared across a batch of `N = 1024` + independent operations via Montgomery's trick (this crate's + [`add_batch`/`double_batch`] + [`field::batch_invert`]; smalljac's + `hecurve_ctx_t` state machine + `ff_parallel_invert`). This is the throughput + metric that matters for the generic-group order computations smalljac targets. + +All numbers below were measured **on the same machine** (Apple Silicon, single +core; smalljac v4.1.3 + ff_poly v1.2.7 ported to arm64), on the same depressed +monic quintic (`f₄ = 0` — required, else smalljac falls back to generic Cantor), +in ns per group operation: + +| field | op | this crate, scalar | smalljac, scalar | this crate, batched | smalljac, batched | +|-------|----|------------------:|-----------------:|-------------------:|------------------:| +| p = 65521 (16-bit) | add | 222 | 100 | 99 | 47 | +| p = 65521 (16-bit) | double | 245 | 105 | 115 | 54 | +| 56-bit prime | add | 673 | 190 | 348 | 48 | +| 56-bit prime | double | 691 | 210 | 384 | 54 | + +Batching helps both implementations (it removes the inversion): this crate's +56-bit add drops 673 → 348 ns (1.9×), smalljac's 190 → 48 ns (4×). But smalljac +is faster in every cell, and the reason is the **field-arithmetic layer**, not +the genus-2 formulas (the [operation counts](#field-operation-counts) — ~26 M, +1 I — are essentially the same). Measured field ops at the 56-bit prime: + +| field operation | this crate (`PrimeField`) | smalljac (ff_poly) | +|-----------------|--------------------------:|-------------------:| +| multiply | 8.2 ns | 4.5 ns | +| inversion | 529 ns | 157 ns | +| amortized inversion (batch of 1024) | ~25 ns | ~10 ns | + +smalljac uses single-word **Montgomery** arithmetic (no division in the +multiply; a tuned binary-GCD inversion), whereas this crate's `PrimeField` uses +schoolbook `u128`-multiply + hardware modulo and an extended-Euclidean inversion +— ~1.8× slower to multiply and ~3.4× slower to invert. Since the scalar group op +is inversion-dominated (529 of 673 ns at 56-bit), that gap is what the batched +path removes, and it is what a Montgomery `PrimeField` would close. The takeaway: +**the explicit formulas and the batched group law are sound and competitive in +operation count; closing the wall-clock gap to smalljac is a field-arithmetic +optimization (Montgomery reduction), not a formula one.** + +The harness (scalar + batched + raw field ops) and arm64 build notes are in +[`benches/smalljac-compare/`](benches/smalljac-compare/); run this crate's side +with `cargo bench -- not_char2`. ## Testing diff --git a/benches/benchmarks.rs b/benches/benchmarks.rs index c21121f..11286ca 100644 --- a/benches/benchmarks.rs +++ b/benches/benchmarks.rs @@ -1,6 +1,6 @@ //! Benchmarks for genus 2 ramified divisor arithmetic operations. -use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion}; +use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput}; use divisor_arithmetic::field::{BinaryExtField, Field, PrimeField}; use divisor_arithmetic::g2::ramified::{arbitrary, char2, not_char2}; @@ -506,10 +506,64 @@ fn char2_benchmarks(c: &mut Criterion) { fn field_benchmarks(c: &mut Criterion) { bench_field_ops::>(c, "F65521"); + bench_field_ops::>(c, "Fp56"); bench_field_ops::>(c, "GF256"); bench_field_ops::>(c, "GF65536"); } +// ============================================================================= +// Batched group law (Montgomery simultaneous inversion) +// ============================================================================= + +fn rand_deg2(rng: &mut impl rand::Rng) -> not_char2::DivisorCoords { + not_char2::DivisorCoords::deg2( + F::random(rng), + F::random(rng), + F::random(rng), + F::random(rng), + ) +} + +/// Amortized throughput of the batched group law: one field inversion is shared +/// across a batch of `n` independent ops via `batch_invert`. This is the metric +/// that matters for smalljac-style generic-group algorithms (and matches its +/// `ctx` + `ff_parallel_invert` path). Criterion reports elements/sec; per-op +/// time = batch_time / n. +fn bench_not_char2_batched(c: &mut Criterion, name: &str, n: usize) { + let mut rng = rand::thread_rng(); + let cc = not_char2::CurveConstants { + f3: F::random(&mut rng), + f2: F::random(&mut rng), + f1: F::random(&mut rng), + f0: F::random(&mut rng), + }; + let pairs: Vec<_> = (0..n) + .map(|_| (rand_deg2::(&mut rng), rand_deg2::(&mut rng))) + .collect(); + let singles: Vec<_> = (0..n).map(|_| rand_deg2::(&mut rng)).collect(); + + let mut g = c.benchmark_group("not_char2_batched"); + g.throughput(Throughput::Elements(n as u64)); + g.bench_with_input( + BenchmarkId::new("add_batch", format!("{name}/N={n}")), + &pairs, + |b, pairs| b.iter(|| not_char2::add_batch(black_box(pairs), black_box(&cc))), + ); + g.bench_with_input( + BenchmarkId::new("dbl_batch", format!("{name}/N={n}")), + &singles, + |b, singles| b.iter(|| not_char2::double_batch(black_box(singles), black_box(&cc))), + ); + g.finish(); +} + +fn batched_benchmarks(c: &mut Criterion) { + // Matched 56-bit width vs smalljac's batched (ctx + ff_parallel_invert) path, + // plus 16-bit for the field-width contrast. + bench_not_char2_batched::>(c, "Fp56", 1024); + bench_not_char2_batched::>(c, "F65521", 1024); +} + criterion_group!( benches, field_benchmarks, @@ -517,5 +571,6 @@ criterion_group!( arbitrary_benchmarks, char2_benchmarks, split_benchmarks, + batched_benchmarks, ); criterion_main!(benches); diff --git a/benches/smalljac-compare/README.md b/benches/smalljac-compare/README.md index 7a87a34..74c0ac1 100644 --- a/benches/smalljac-compare/README.md +++ b/benches/smalljac-compare/README.md @@ -1,11 +1,23 @@ # 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. +(`hecurve_g2_compose` = add, `hecurve_g2_square` = double) so it can be put next +to this crate's `g2::ramified::not_char2` `cargo bench` numbers. It reports: + +- **raw field ops** — multiply, inversion, and amortized (parallel) inversion; +- **scalar** group law — `ctx = NULL`, one inversion per op; +- **batched** group law — the `hecurve_ctx_t` state machine + `ff_parallel_invert` + (Montgomery's trick), one inversion shared across `N = 1024` ops, matching this + crate's `add_batch`/`double_batch`. + +See the "Wall-clock comparison with smalljac" section of the top-level +[`README.md`](../../README.md) for the results and analysis. + +> **Important:** the curve must use a *depressed* quintic (`f[4] = 0`, as the +> harness does). With `f[4] ≠ 0`, `hecurve_g2_compose` silently reverts to the +> slow generic **Cantor** path and the measurement is ~5–8× too slow and +> meaningless. This crate's `not_char2` model is `y² = x⁵ + f₃x³ + …`, i.e. +> `f₄ = 0`, so matching it is also the correct comparison. smalljac and ff_poly are **not vendored** here (they are GPL; this crate is MIT). Download them yourself: diff --git a/benches/smalljac-compare/bench.c b/benches/smalljac-compare/bench.c index 8128101..e626b20 100644 --- a/benches/smalljac-compare/bench.c +++ b/benches/smalljac-compare/bench.c @@ -1,14 +1,18 @@ -/* 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. +/* Micro-benchmark for smalljac's genus-2 Jacobian group law, scalar and + * batched, to compare against the Rust `divisor-arithmetic` cargo-bench numbers + * for g2::ramified::not_char2. See ../../README.md ("smalljac comparison"). * - * 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. + * Imaginary/ramified genus-2 model: y^2 = f(x), deg f = 5 (monic, depressed: + * f[4] = 0 — REQUIRED, else hecurve reverts to the slow generic Cantor path and + * the comparison is meaningless; the Rust not_char2 model also uses f4 = 0). * - * 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). + * - scalar: ctx = NULL => one field inversion per op. + * - batched: ctx state machine + ff_parallel_invert => one inversion per batch + * (Montgomery's trick) — the throughput metric that matters for the + * generic-group algorithms smalljac is built for. + * + * Build: see README.md in this directory (needs an external smalljac + ff_poly + * checkout; on arm64 apply the two portability patches listed there). */ #include #include @@ -19,7 +23,8 @@ #include "hecurve.h" #include "cstd.h" -#define POOL 256 +#define POOL 1100 +#define BATCH 1024 static double now_ns(void) { struct timespec ts; @@ -27,64 +32,130 @@ static double now_ns(void) { return ts.tv_sec * 1e9 + ts.tv_nsec; } -static void run(unsigned long p, const char *label) { - ff_setup_ui(p); +static ff_t U[POOL][HECURVE_GENUS + 1], V[POOL][HECURVE_GENUS]; +static ff_t f[HECURVE_DEGREE + 1]; +static hecurve_ctx_t ctx[BATCH]; +static ff_t RU[BATCH][HECURVE_GENUS + 1], RV[BATCH][HECURVE_GENUS]; +static ff_t invs[BATCH]; +static int idx[BATCH]; - /* monic degree-5 f with random coefficients f[0..4], f[5]=1, f[6]=0 */ - ff_t f[HECURVE_DEGREE + 1]; +static int build_pool(unsigned long p) { + ff_setup_ui(p); 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]; + for (int i = 0; i < 4; i++) _ff_random(f[i]); /* f0..f3 random, f4 = 0 */ + _ff_set_one(f[5]); /* monic degree 5 */ 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 (_ff_zero(u[2]) || !_ff_one(u[2])) continue; 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; } + return n; +} + +static void run_fieldops(unsigned long p, const char *label) { + ff_setup_ui(p); + ff_t a, b, acc, t; + _ff_set_ui(a, 12345 % p); _ff_set_ui(b, 67891 % p); _ff_set(acc, a); + long iters = 50000000; + double t0 = now_ns(); + for (long k = 0; k < iters; k++) ff_mult(acc, acc, b); + double m_ns = (now_ns() - t0) / iters; + unsigned long sink = _ff_get_ui(acc); + _ff_set_ui(acc, 12345 % p); + long it2 = 5000000; + t0 = now_ns(); + for (long k = 0; k < it2; k++) { ff_invert(t, acc); _ff_mult(acc, t, b); } + double i_ns = (now_ns() - t0) / it2; + sink += _ff_get_ui(acc); + static ff_t xs[BATCH], zs[BATCH]; + for (int i = 0; i < BATCH; i++) _ff_random(xs[i]); + long it3 = 20000; + t0 = now_ns(); + for (long k = 0; k < it3; k++) ff_parallel_invert(zs, xs, BATCH); + double pi_ns = (now_ns() - t0) / ((double)it3 * BATCH); + sink += _ff_get_ui(zs[0]); + printf("fieldop %-12s p=%-20lu M=%5.2f ns I=%6.2f ns batchedI=%5.2f ns/elem (I/M=%.1f)\n", + label, p, m_ns, i_ns, pi_ns, i_ns / m_ns); + (void)sink; +} +static void run_scalar(unsigned long p, const char *label) { + int n = build_pool(p); + if (n < 8) { printf("scalar %-12s p=%lu: only %d divisors\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]); - } + 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]); - } + 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("scalar %-12s p=%-20lu add=%7.1f ns double=%7.1f ns\n", label, p, add_ns, dbl_ns); + (void)sink; +} + +static void run_batched(unsigned long p, const char *label) { + int n = build_pool(p); + if (n < BATCH + 1) { printf("batched %-12s p=%lu: only %d divisors\n", label, p, n); return; } + unsigned long sink = 0; + long rounds = 2000; - 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); + double t0 = now_ns(); + for (long r = 0; r < rounds; r++) { + int m = 0; + for (int i = 0; i < BATCH; i++) { + ctx[i].state = 0; + int done = hecurve_g2_compose(RU[i], RV[i], U[i], V[i], U[i+1], V[i+1], f, &ctx[i]); + if (!done) { invs[m] = ctx[i].invert; idx[m] = i; m++; } + } + ff_parallel_invert(invs, invs, m); + for (int j = 0; j < m; j++) { + int i = idx[j]; + _ff_set(ctx[i].invert, invs[j]); + hecurve_g2_compose(RU[i], RV[i], U[i], V[i], U[i+1], V[i+1], f, &ctx[i]); + } + sink += _ff_get_ui(RU[0][0]); + } + double add_ns = (now_ns() - t0) / ((double)rounds * BATCH); + + t0 = now_ns(); + for (long r = 0; r < rounds; r++) { + int m = 0; + for (int i = 0; i < BATCH; i++) { + ctx[i].state = 0; + int done = hecurve_g2_square(RU[i], RV[i], U[i], V[i], f, &ctx[i]); + if (!done) { invs[m] = ctx[i].invert; idx[m] = i; m++; } + } + ff_parallel_invert(invs, invs, m); + for (int j = 0; j < m; j++) { + int i = idx[j]; + _ff_set(ctx[i].invert, invs[j]); + hecurve_g2_square(RU[i], RV[i], U[i], V[i], f, &ctx[i]); + } + sink += _ff_get_ui(RU[0][0]); + } + double dbl_ns = (now_ns() - t0) / ((double)rounds * BATCH); + printf("batched %-12s p=%-20lu add=%7.1f ns double=%7.1f ns (N=%d)\n", + label, p, add_ns, dbl_ns, BATCH); + (void)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"); + run_fieldops(65521UL, "16-bit"); + run_fieldops(72057594037927931UL, "56-bit Fp56"); + run_scalar(65521UL, "16-bit"); + run_batched(65521UL, "16-bit"); + run_scalar(72057594037927931UL, "56-bit Fp56"); + run_batched(72057594037927931UL, "56-bit Fp56"); return 0; } diff --git a/src/field.rs b/src/field.rs index d93e4ca..68f5a2e 100644 --- a/src/field.rs +++ b/src/field.rs @@ -80,6 +80,49 @@ pub trait Field: fn random(rng: &mut R) -> Self; } +/// Invert every element of `xs` in place using Montgomery's trick +/// (simultaneous inversion): one field inversion plus `~3(n−1)` multiplications +/// for the whole slice, instead of `n` inversions. +/// +/// Zero entries are left as zero (their inverse is undefined and simply skipped), +/// so callers can pass a slice that mixes invertible and zero values. This is the +/// same batching smalljac performs with `ff_parallel_invert` to amortize the one +/// expensive inversion across many independent group operations. +/// +/// ``` +/// use divisor_arithmetic::field::{batch_invert, Field, PrimeField}; +/// type F = PrimeField<65521>; +/// let mut xs = [F::new(2), F::new(3), F::new(0), F::new(7)]; +/// batch_invert(&mut xs); +/// assert_eq!(xs[0], F::new(2).inv()); +/// assert_eq!(xs[2], F::new(0)); // zero is left untouched +/// ``` +pub fn batch_invert(xs: &mut [F]) { + let n = xs.len(); + if n == 0 { + return; + } + // Forward pass: prefix[i] = product of the nonzero entries in xs[0..i]. + let mut prefix = Vec::with_capacity(n); + let mut acc = F::one(); + for &x in xs.iter() { + prefix.push(acc); + if !x.is_zero() { + acc *= x; + } + } + // One inversion for the whole batch: acc = 1 / (product of all nonzero xs). + let mut acc = acc.inv(); + // Backward pass: recover each individual inverse. + for i in (0..n).rev() { + if !xs[i].is_zero() { + let inv_i = acc * prefix[i]; + acc *= xs[i]; + xs[i] = inv_i; + } + } +} + /// A simple prime field implementation using u64 arithmetic. /// Suitable for small primes where p² fits in u128. #[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] diff --git a/src/g2/ramified/not_char2.rs b/src/g2/ramified/not_char2.rs index 1b52807..fcc1457 100644 --- a/src/g2/ramified/not_char2.rs +++ b/src/g2/ramified/not_char2.rs @@ -12,7 +12,7 @@ //! //! Based on: Sebastian Lindner, 2019 -use crate::field::Field; +use crate::field::{batch_invert, Field}; /// Curve constants for a not-char-2 ramified genus 2 curve. /// @@ -515,6 +515,289 @@ pub fn double(d: &DivisorCoords, cc: &CurveConstants) -> Divisor } } +// ============================================================================= +// Batched inversion (Montgomery's trick) +// ============================================================================= +// +// Each generic degree-2 group operation performs exactly one field inversion. +// When many *independent* operations are evaluated together, that inversion can +// be amortized: split each operation into a pre-inversion phase (everything up +// to the single `inv`) and a post-inversion phase (everything after), collect +// the to-be-inverted values from the whole batch, invert them all with one field +// inversion via [`batch_invert`], then finish each operation. This is exactly the +// strategy smalljac uses through its `hecurve_ctx_t` state machine + the +// `ff_parallel_invert` routine, and is what makes the affine group law fast in +// the generic-group algorithms (e.g. BSGS order computations) that use it. +// +// Only the fully generic degree-2 path is batched (the rare special branches — +// identity inputs, `d = 0`, degree-dropping `sp1 = 0` — fall back to the direct +// [`add`]/[`double`], which do their own single inversion). On random inputs the +// special cases occur with negligible probability, so the batch is essentially +// uniform. + +/// Captured state of a generic degree-2 [`deg2_add`], taken just before its one +/// inversion, so the result can be finished by [`deg2_add_post`] once the +/// inverse is known. Produced by [`deg2_add_pre`]. +#[derive(Clone, Copy, Debug)] +pub struct Deg2AddPre { + u1: F, + u0: F, + v1: F, + v0: F, + up1: F, + m1: F, + m3: F, + sp1: F, + sp0: F, + d: F, +} + +/// Phase 1 of a generic degree-2 + degree-2 addition: compute everything up to +/// the single inversion. Returns `Some((state, to_invert))` for the generic case +/// (`d ≠ 0` and `sp1 ≠ 0`), where `to_invert = d·sp1` is the value the caller +/// must invert before calling [`deg2_add_post`]; returns `None` when a special +/// branch applies, in which case the caller should fall back to [`deg2_add`]. +#[inline] +#[allow(clippy::too_many_arguments)] +pub fn deg2_add_pre( + u1: F, + u0: F, + v1: F, + v0: F, + up1: F, + up0: F, + vp1: F, + vp0: F, + _cc: &CurveConstants, +) -> Option<(Deg2AddPre, F)> { + let m3 = up1 - u1; + let m4 = u0 - up0; + let m1 = m4 + up1 * m3; + let m2 = -up0 * m3; + let d = m1 * m4 - m2 * m3; + if d.is_zero() { + return None; + } + let r0 = vp0 - v0; + let r1 = vp1 - v1; + let sp1 = r0 * m3 + r1 * m4; + let sp0 = r0 * m1 + r1 * m2; + if sp1.is_zero() { + return None; + } + Some(( + Deg2AddPre { + u1, + u0, + v1, + v0, + up1, + m1, + m3, + sp1, + sp0, + d, + }, + d * sp1, + )) +} + +/// Phase 2 of a generic degree-2 addition: finish using `inv = (d·sp1)⁻¹` +/// (the inverse of the value returned by [`deg2_add_pre`]). Identical arithmetic +/// to the `sp1 ≠ 0` tail of [`deg2_add`]. +#[inline] +pub fn deg2_add_post(pre: &Deg2AddPre, inv: F) -> DivisorCoords { + let Deg2AddPre { + u1, + u0, + v1, + v0, + up1, + m1, + m3, + sp1, + sp0, + d, + } = *pre; + + let w1 = inv; // (d·sp1)⁻¹ + let w2 = w1 * d; + let w3 = w2 * d; + let w4 = w3.square(); + let s1 = w1 * sp1.square(); + let spp0 = sp0 * w2; + + let t1 = spp0 - m3; + let t2 = t1 - w4; + let t3 = w3 * v1; + let upp1 = spp0 + t2; + let upp0 = spp0 * (t1 - m3) + m1 + t3 + t3 + w4 * (u1 + up1); + + let t0 = upp0 - u0; + let t1 = u1 - upp1; + let vpp1 = s1 * (t1 * t2 + t0) - v1; + let vpp0 = s1 * (spp0 * t0 + upp0 * t1) - v0; + + DivisorCoords::deg2(upp1, upp0, vpp1, vpp0) +} + +/// Captured state of a generic degree-2 [`deg2_dbl`], taken just before its one +/// inversion. Produced by [`deg2_dbl_pre`], consumed by [`deg2_dbl_post`]. +#[derive(Clone, Copy, Debug)] +pub struct Deg2DblPre { + u1: F, + u0: F, + v1: F, + v0: F, + sp1: F, + sp0: F, + d: F, +} + +/// Phase 1 of a generic degree-2 doubling. Returns `Some((state, to_invert))` +/// with `to_invert = d·sp1` for the generic case, else `None` (fall back to +/// [`deg2_dbl`]). +#[inline] +pub fn deg2_dbl_pre( + u1: F, + u0: F, + v1: F, + v0: F, + cc: &CurveConstants, +) -> Option<(Deg2DblPre, F)> { + let CurveConstants { f3, f2, .. } = *cc; + + let m3 = -v1 - v1; + let m4 = v0 + v0; + let m1 = m4 + m3 * u1; + let m2 = -m3 * u0; + let d = m4 * m1 - m2 * m3; + if d.is_zero() { + return None; + } + + let t0 = u1.square(); + let t1 = f3 + t0; + let t2 = u0 + u0; + let t3 = t1 - t2; + let r1 = t0 + t0 + t3; + let r0 = u1 * (t2 - t3) + f2 - v1.square(); + + let sp0 = r0 * m1 + r1 * m2; + let sp1 = r0 * m3 + r1 * m4; + if sp1.is_zero() { + return None; + } + Some(( + Deg2DblPre { + u1, + u0, + v1, + v0, + sp1, + sp0, + d, + }, + d * sp1, + )) +} + +/// Phase 2 of a generic degree-2 doubling: finish using `inv = (d·sp1)⁻¹`. +/// Identical arithmetic to the `sp1 ≠ 0` tail of [`deg2_dbl`]. +#[inline] +pub fn deg2_dbl_post(pre: &Deg2DblPre, inv: F) -> DivisorCoords { + let Deg2DblPre { + u1, + u0, + v1, + v0, + sp1, + sp0, + d, + } = *pre; + + let w1 = inv; // (d·sp1)⁻¹ + let w2 = w1 * d; + let w3 = w2 * d; + let w4 = w3.square(); + let s1 = w1 * sp1.square(); + let spp0 = sp0 * w2; + + let t2 = spp0 - w4; + let t3 = w3 * v1 + w4 * u1; + let upp1 = spp0 + t2; + let upp0 = spp0.square() + t3 + t3; + + let t0 = upp0 - u0; + let t1 = u1 - upp1; + let vpp1 = s1 * (t1 * t2 + t0) - v1; + let vpp0 = s1 * (spp0 * t0 + upp0 * t1) - v0; + + DivisorCoords::deg2(upp1, upp0, vpp1, vpp0) +} + +/// Add many divisor pairs with a single field inversion for the whole batch. +/// +/// Equivalent to mapping [`add`] over `pairs`, but the one inversion each generic +/// degree-2 addition would perform is amortized across the batch via +/// [`batch_invert`] (Montgomery's trick). The result is identical, element for +/// element, to calling [`add`] on each pair. +pub fn add_batch( + pairs: &[(DivisorCoords, DivisorCoords)], + cc: &CurveConstants, +) -> Vec> { + let mut out = vec![DivisorCoords::identity(); pairs.len()]; + let mut pres: Vec<(usize, Deg2AddPre)> = Vec::with_capacity(pairs.len()); + let mut to_invert: Vec = Vec::with_capacity(pairs.len()); + + for (i, (d1, d2)) in pairs.iter().enumerate() { + if d1.degree() == 2 && d2.degree() == 2 { + if let Some((pre, t)) = + deg2_add_pre(d1.u1, d1.u0, d1.v1, d1.v0, d2.u1, d2.u0, d2.v1, d2.v0, cc) + { + pres.push((i, pre)); + to_invert.push(t); + continue; + } + } + out[i] = add(d1, d2, cc); + } + + batch_invert(&mut to_invert); + for ((i, pre), &inv) in pres.iter().zip(to_invert.iter()) { + out[*i] = deg2_add_post(pre, inv); + } + out +} + +/// Double many divisors with a single field inversion for the whole batch. +/// Result is identical, element for element, to calling [`double`] on each. +pub fn double_batch( + ds: &[DivisorCoords], + cc: &CurveConstants, +) -> Vec> { + let mut out = vec![DivisorCoords::identity(); ds.len()]; + let mut pres: Vec<(usize, Deg2DblPre)> = Vec::with_capacity(ds.len()); + let mut to_invert: Vec = Vec::with_capacity(ds.len()); + + for (i, d) in ds.iter().enumerate() { + if d.degree() == 2 { + if let Some((pre, t)) = deg2_dbl_pre(d.u1, d.u0, d.v1, d.v0, cc) { + pres.push((i, pre)); + to_invert.push(t); + continue; + } + } + out[i] = double(d, cc); + } + + batch_invert(&mut to_invert); + for ((i, pre), &inv) in pres.iter().zip(to_invert.iter()) { + out[*i] = deg2_dbl_post(pre, inv); + } + out +} + #[cfg(test)] mod tests { use super::*; @@ -586,4 +869,49 @@ mod tests { let id = DivisorCoords::identity(); assert_eq!(double(&id, &cc), id); } + + #[test] + fn batched_matches_scalar() { + use crate::field::PrimeField; + use rand::rngs::StdRng; + use rand::SeedableRng; + type F = PrimeField<65521>; + + let mut rng = StdRng::seed_from_u64(7); + let cc = CurveConstants { + f3: F::random(&mut rng), + f2: F::random(&mut rng), + f1: F::random(&mut rng), + f0: F::random(&mut rng), + }; + let rand_d2 = |rng: &mut StdRng| { + DivisorCoords::deg2( + F::random(rng), + F::random(rng), + F::random(rng), + F::random(rng), + ) + }; + + // add_batch == add, element for element (covers generic + any special + // branches that fall through to the direct path). + let pairs: Vec<_> = (0..1000) + .map(|_| (rand_d2(&mut rng), rand_d2(&mut rng))) + .collect(); + let batched = add_batch(&pairs, &cc); + for (i, (d1, d2)) in pairs.iter().enumerate() { + assert_eq!(batched[i], add(d1, d2, &cc), "add_batch mismatch at {i}"); + } + + // double_batch == double, element for element. + let singles: Vec<_> = pairs.iter().map(|(a, _)| *a).collect(); + let batched_dbl = double_batch(&singles, &cc); + for (i, d) in singles.iter().enumerate() { + assert_eq!( + batched_dbl[i], + double(d, &cc), + "double_batch mismatch at {i}" + ); + } + } }