engineering/lumbda
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
lumbda/quantum/mod-karatsuba.lsp
russell@unturf.com 1665893321
factory + quantum + sweep-doctrine: AGPLv3 share-back from foxhop ecdsa 
29 new files publish factory infra (V2 autoscaler with live VRAM
sampling + EWMA peak tracking, HUGE solo-dispatch, two-tier DLQ/rDLQ
classifier + retry), general quantum circuit primitives (Cuccaro
ripple-carry adder, Clifford gate library, Clifford tableau simulator,
mod-arith family, dialog GCD reversible inverse, Karatsuba multiplier,
Solinas fast reduction), and a TCRAUDT reducer harness. Originally
developed in ~/git/www.foxhop.net/ecdsa/ for secp256k1 attack-surface
research; published upstream as obligated by AGPLv3.

Parametrization contract at factory/CONTRACT.md. Consumers export
LUMBDA_REPO_DIR + LUMBDA_QUEUE_DIR + LUMBDA_BACKEND_CMD + LUMBDA_EMITTER_CMD
then exec factory scripts. No fork-and-modify; single source of truth
upstream.

Integration tests gate 7 V2 defect classes that wedged a live factory
on 2026-06-12 (skewed-demand starve, zero-floor reservation,
multi-tier greedy, +-25%% damping, cold-start ramp, DLQ surge halve,
post-damp CPU ceiling) + 28 DLQ classifier cases (auto-retry vs
escalate partition) + bash -n syntax lint across every script.

GPU backend stays in consumer trees; rationale in
factory/GPU-BACKEND-NOTE.md. Bend wire protocol + gpu-worker.lsp
already upstream at examples/cuda-fanout/.

make factory-lint                bash -n on every factory/*.sh
make test-integration            V2 reducer + DLQ classifier + syntax gate
make sweep-doctrine              TCRAUDT reducer gate (serial)
make sweep-doctrine-parallel     xargs -P fan-out

Verified on neoblanka: factory-lint 12 scripts PASS; test-integration
14 V2 cases + 28 DLQ classifier cases + 12 syntax cases all PASS.
2026-06-14 10:37:35 -04:00

838 lines
36 KiB
Text
Raw Permalink Blame History

;;; mod-karatsuba.lsp — Phase B step 8: Karatsuba reversible wide-multiply.
;;;
;;; Replaces the wide-product stage of mod-mul-solinas! (the Litinski
;;; schoolbook-mul-into-addsub!) with a Karatsuba recursion. Below the
;;; threshold *karatsuba-threshold*, falls through to the existing
;;; schoolbook primitive — so at p=11 (n=4) and p=251 (n=8) the karatsuba
;;; path produces byte-identical output to mod-mul-solinas!. The win
;;; shows up at n+1 >= 18.
;;;
;;; Algorithm (single-level expansion):
;;;
;;; karatsuba-wide!(a, b, dst, n):
;;; if n <= THRESHOLD: schoolbook-wide!(a, b, dst, n) ; fallback
;;; else:
;;; m_lo = ceil(n/2), m_hi = n - m_lo ; m_hi <= m_lo
;;; a_lo = a[0..m_lo), a_hi = a[m_lo..n)
;;; b_lo = b[0..m_lo), b_hi = b[m_lo..n)
;;;
;;; ;; Compute three sub-products into fresh 2*m wide accumulators.
;;; alloc z0 (2*m_lo), z2 (2*m_lo), sum_a (m_lo+1), sum_b (m_lo+1),
;;; z1 (2*(m_lo+1))
;;; karatsuba-wide!(a_lo, b_lo, z0, m_lo) ; z0 = a_lo*b_lo
;;; karatsuba-wide!(a_hi, b_hi, z2-low(2*m_hi),m_hi) ; z2 high bits |0>
;;; sum_a = a_lo (copy) + a_hi (add)
;;; sum_b = b_lo (copy) + b_hi (add)
;;; karatsuba-wide!(sum_a, sum_b, z1, m_lo+1)
;;;
;;; ;; Combine into dst.
;;; ;; dst += z0 at offset 0 (2*m_lo bits)
;;; ;; dst += z2 at offset 2*m_lo (2*m_hi bits)
;;; ;; dst += z1 at offset m_lo (2*(m_lo+1) bits)
;;; ;; dst -= z0 at offset m_lo (2*m_lo bits)
;;; ;; dst -= z2 at offset m_lo (2*m_hi bits)
;;;
;;; ;; Uncompute every sub-product & sum register.
;;; karatsuba-wide!-INVERSE(sum_a, sum_b, z1, m_lo+1)
;;; uncompute sum_b (CX b_hi out, CX b_lo out — reverse order)
;;; uncompute sum_a
;;; karatsuba-wide!-INVERSE(a_hi, b_hi, z2-low, m_hi)
;;; karatsuba-wide!-INVERSE(a_lo, b_lo, z0, m_lo)
;;; free buffers
;;;
;;; Op-count scaling: at the leaves we still pay schoolbook (n_leaf^2
;;; controlled-add-subtracts). With k = log_2(n/THRESHOLD) levels we get
;;; 3^k leaf calls instead of 4^k. At THRESHOLD=8, n=32 -> k=2 -> 9 leaves
;;; vs 16 schoolbook (44% Toffoli reduction in the multiply stage); n=128
;;; -> k=4 -> 81 leaves vs 256 schoolbook (68% reduction).
;;;
;;; Ancilla overhead: each recursion allocates 4 width-2m_lo buffers + 2
;;; width-(m_lo+1) buffers = ~6n+4 ancilla qubits per level. Total
;;; cumulative ~6n*log2(n) ancilla — pays for itself well above n=16.
;;;
;;; Cross-tier portal validation: the schoolbook fallback below threshold
;;; means at p=11 / p=251 mod-mul-karatsuba! emits THE SAME gate sequence
;;; as mod-mul-solinas! — sweep-001/002 oracles still match.
(load "quantum/gates.lsp")
(load "quantum/adder.lsp")
(load "quantum/mod-arith.lsp")
(load "quantum/mod-solinas.lsp")
;;; ── threshold knob (caller may set! before building) ───────────
(define *karatsuba-threshold* 8)
;; n <= *karatsuba-threshold* falls through to schoolbook. Default 8
;; chosen so n+1 in {5, 9} (the sweep-001/002 widths) skip Karatsuba
;; entirely → byte-identity preserved by construction.
;;; ── ancilla-name generator ────────────────────────────────────
;;;
;;; Recursion needs fresh ancilla names per level. We append the depth
;;; suffix to each base name; depth threads through every call.
(define (kara-name base depth)
(string->symbol (string-append (symbol->string base) "-d"
(number->string depth))))
(define (kara-name-w base depth width)
"Like kara-name but ALSO tags the width. Used by primitives whose
ancilla width varies between sibling calls at the same depth — the
simulator's widths-hash retains the LAST-set width per name, so
re-allocating the same name at different widths within one circuit
lifetime crashes the reverse-pass register re-creation. Width-tagging
sidesteps the collision."
(string->symbol (string-append (symbol->string base) "-d"
(number->string depth) "-w"
(number->string width))))
;;; ── public entry: karatsuba-wide! ─────────────────────────────
;;;
;;; dst (2n bits |0> in) receives a*b. a and b are n-wide registers
;;; (passed by name + offset for slicing). Caller allocates dst, a, b,
;;; and the schoolbook-shared scratch (sb-low, sb-xext, sb-const-tmp,
;;; cin). All ancilla return to |0>.
(define (karatsuba-wide! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth)
"Reversible wide multiply: dst[dst-off..dst-off+2n) := a*b. a, b
preserved; dst |0> in. Falls through to kara-leaf-mul! below threshold
or when n <= 3 (recursion-floor guard: at n=3, m_sum=3=n so the sum
sub-multiply would re-enter karatsuba-step at the same n → infinite
recursion). Threshold default 8 keeps us well above this guard."
(cond
((= n 0) #t)
((or (<= n *karatsuba-threshold*) (<= n 3))
(kara-leaf-mul! c a-reg a-off b-reg b-off n
dst-reg dst-off
cin-reg cin-idx
depth))
(else
(karatsuba-step! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth))))
(define (karatsuba-wide-inverse! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth)
"Inverse: zero-out dst that currently holds a*b, restoring |0>."
(cond
((= n 0) #t)
((or (<= n *karatsuba-threshold*) (<= n 3))
(kara-leaf-mul-inverse! c a-reg a-off b-reg b-off n
dst-reg dst-off
cin-reg cin-idx
depth))
(else
(karatsuba-step-inverse! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth))))
;;; ── leaf multiplier: kara-leaf-mul! ──
;;;
;;; Reversible "schoolbook" multiplier with unique-per-depth ancilla
;;; names. Why not use the existing Litinski schoolbook (mod-arith.lsp
;;; schoolbook-mul-into-addsub!)?
;;;
;;; The Litinski primitive allocates `sb-xfull` internally with width
;;; 2n+1 = function of leaf n. Different Karatsuba leaves run schoolbook
;;; at different n values within the SAME circuit lifetime (e.g. m_lo
;;; sub-mult at n=m_lo, sum_a*sum_b sub-mult at n=m_lo+1). When the
;;; simulator's reverse pass replays alloc ops, it uses the widths-hash
;;; value at that point, but the widths hash retains the LAST forward
;;; set value — so for an op-pair (alloc sb-xfull w5, ..., free sb-xfull)
;;; followed later by (alloc sb-xfull w7, ..., free), the reverse pass
;;; re-creates the w5 register at width 7, mis-sizing the bit vector
;;; and crashing on the next vector-ref.
;;;
;;; Cleanest workaround: leaf multiplier has its OWN ancilla names that
;;; are depth-tagged, so the width per name stays consistent.
;;;
;;; Algorithm (O(n²) controlled-add):
;;; for i in 0..n:
;;; for k in 0..n: tmp[k] ^= b[b-off+i] AND a[a-off+k] (CCX)
;;; dst[dst-off+i..dst-off+i+n) += tmp (cuccaro-add-offset width=n)
;;; uncompute tmp (same CCX self-inverse)
;;;
;;; Reversible, dst |0> in / dst = a*b out. tmp is depth-tagged.
;;; Width-cost analysis: tmp width n.
(define (kara-leaf-mul! c a-reg a-off b-reg b-off n
dst-reg dst-off
cin-reg cin-idx
depth)
"dst[dst-off..dst-off+2n) += a*b. a, b each n wide (preserved).
dst-reg must be >= dst-off + 2n wide. Uses kara-leaf-tmp-d<depth>
ancilla of width n+1 (top bit |0> carry-extension to propagate the
add's carry-out into dst[dst-off+i+n)).
Algorithm:
for i in 0..n:
for k in 0..n: tmp[k] ^= b[i] AND a[k] (CCX)
dst[dst-off+i..dst-off+i+n+1) += tmp[0..n+1) (cuccaro-add width n+1)
uncompute tmp
The width-n+1 add propagates the carry from bit (i+n-1) into bit
(i+n) cleanly. tmp[n] stays |0> across each row (just the carry slot)."
(let ((tmpn (kara-name-w 'kara-leaf-tmp depth n)))
(alloc! c tmpn (+ n 1))
(let row-loop ((i 0))
(when (< i n)
;; tmp ^= b[i] AND a[k]
(let loop ((k 0))
(when (< k n)
(gate-ccx! c b-reg (+ b-off i) a-reg (+ a-off k) tmpn k)
(loop (+ k 1))))
;; dst[dst-off+i..+i+n+1) += tmp[0..n+1) — width n+1 lets carry propagate
(cuccaro-add-offset! c tmpn 0 dst-reg (+ dst-off i)
cin-reg cin-idx (+ n 1))
;; uncompute tmp
(let loop ((k 0))
(when (< k n)
(gate-ccx! c b-reg (+ b-off i) a-reg (+ a-off k) tmpn k)
(loop (+ k 1))))
(row-loop (+ i 1))))
(free! c tmpn)))
(define (kara-leaf-mul-inverse! c a-reg a-off b-reg b-off n
dst-reg dst-off
cin-reg cin-idx
depth)
"Inverse: dst[dst-off..dst-off+2n) -= a*b. Same convention as
kara-leaf-mul! with cuccaro-add replaced by cuccaro-sub. Walks rows
in REVERSE so the inverse mirrors the forward operation pair-by-pair."
(let ((tmpn (kara-name-w 'kara-leaf-tmp depth n)))
(alloc! c tmpn (+ n 1))
(let row-loop ((i (- n 1)))
(when (>= i 0)
(let loop ((k 0))
(when (< k n)
(gate-ccx! c b-reg (+ b-off i) a-reg (+ a-off k) tmpn k)
(loop (+ k 1))))
(cuccaro-sub-offset! c tmpn 0 dst-reg (+ dst-off i)
cin-reg cin-idx (+ n 1))
(let loop ((k 0))
(when (< k n)
(gate-ccx! c b-reg (+ b-off i) a-reg (+ a-off k) tmpn k)
(loop (+ k 1))))
(row-loop (- i 1))))
(free! c tmpn)))
;;; ── karatsuba-step! — one recursion level ─────────────────────
;;;
;;; Splits a, b at m_lo := ceil(n/2). m_hi := n - m_lo (m_hi <= m_lo).
;;; Allocates buffers, recurses three times, combines into dst, uncomputes.
(define (karatsuba-step! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth)
(let* ((m-lo (quotient (+ n 1) 2)) ; ceil(n/2)
(m-hi (- n m-lo)) ; floor(n/2) <= m-lo
(m-sum (+ m-lo 1)) ; sum register width
(z0-w (* 2 m-lo))
;; z2 has true width 2*m_hi but we allocate 2*m_lo for uniform
;; offset arithmetic — high (2*m_lo - 2*m_hi) bits stay |0>.
(z2-w (* 2 m-hi))
(z1-w (* 2 m-sum))
;; Width-tag every ancilla name. Sibling recursions at the same
;; depth but different n (e.g. m_lo vs m_sum sub-mults) all hit
;; this same step!; without width tags they'd collide on the
;; simulator widths-hash.
(z0n (kara-name-w 'kara-z0 depth n))
(z2n (kara-name-w 'kara-z2 depth n))
(san (kara-name-w 'kara-sa depth n))
(sbn (kara-name-w 'kara-sb depth n))
(z1n (kara-name-w 'kara-z1 depth n))
(next-depth (+ depth 1)))
(alloc! c z0n z0-w)
(alloc! c z2n z2-w)
(alloc! c san m-sum)
(alloc! c sbn m-sum)
(alloc! c z1n z1-w)
;; ── compute z0 := a_lo * b_lo (m_lo bits each → 2*m_lo product) ──
(karatsuba-wide! c a-reg a-off b-reg b-off m-lo
z0n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
;; ── compute z2 := a_hi * b_hi (m_hi bits each → 2*m_hi product) ──
(karatsuba-wide! c a-reg (+ a-off m-lo) b-reg (+ b-off m-lo) m-hi
z2n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
;; ── sum_a := a_lo + a_hi via a width-m_sum cuccaro-add ──
;;
;; Naive approach (cuccaro-add-offset width=m_hi of a_hi into san)
;; TRUNCATES the carry into san[m_lo] = san[m_sum-1]. Instead, pad
;; a_hi into a fresh m_sum-wide register (high bits |0>), copy a_lo
;; into san low m_lo bits, then do a full width-m_sum cuccaro-add.
;; Carry out cleanly into san[m_sum-1].
(let ((ahpn (kara-name-w 'kara-ahp depth n))
(bhpn (kara-name-w 'kara-bhp depth n)))
(alloc! c ahpn m-sum)
(alloc! c bhpn m-sum)
;; sum_a = a_lo (copy)
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c a-reg (+ a-off k) san k)
(loop (+ k 1))))
;; pad a_hi into ahpn low m_hi bits
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
;; sum_a += a_hi_padded (width m_sum, carries cleanly)
(cuccaro-add! c ahpn san cin-reg cin-idx m-sum)
;; sum_b = b_lo (copy)
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c b-reg (+ b-off k) sbn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(cuccaro-add! c bhpn sbn cin-reg cin-idx m-sum)
;; uncompute padded copies (CX self-inverse)
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(free! c bhpn)
(free! c ahpn))
;; ── z1 := sum_a * sum_b ──
(karatsuba-wide! c san 0 sbn 0 m-sum
z1n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
;; ── combine into dst ──
;; Karatsuba: a*b = z0 + (z1 - z0 - z2)*2^m_lo + z2*2^(2*m_lo)
;;
;; Naive direct adds/subs into dst at NARROW widths LOSE BORROW into
;; higher dst bits — subtracting z0 (width 2*m_lo) at offset m_lo
;; truncates inside the slice [m_lo .. m_lo+2*m_lo) instead of letting
;; the borrow propagate up. Standard fix: compute the middle term in
;; a SEPARATE accumulator (here z1n itself) at full width 2*m_sum,
;; where 2*m_sum >= 2*m_lo >= 2*m_hi → borrows stay inside z1n. Then
;; add z1n into dst at offset m_lo.
;;
;; Sequence (z1n holds z1 = sum_a * sum_b on entry):
;; 1. z1n -= z0 (sub at width 2*m_sum — z0 pads with zero high bits)
;; 2. z1n -= z2
;; 3. dst[0..2*m_lo) += z0
;; 4. dst[2*m_lo..2*m_lo+2*m_hi) += z2
;; 5. dst[m_lo..m_lo+2*m_sum) += z1n (the middle term in full)
;; 6. z1n += z2 (restore z1n = z1 for the inverse step)
;; 7. z1n += z0
;;
;; Steps 1-2 and 6-7 use the SAME width-2*m_sum sub via padded copies
;; (z0 and z2 are narrower; we copy into a 2*m_sum-wide temp).
(let ((z0extn (kara-name-w 'kara-z0ext depth n))
(z2extn (kara-name-w 'kara-z2ext depth n)))
(alloc! c z0extn z1-w)
(alloc! c z2extn z1-w)
;; copy z0 into z0extn (low z0-w bits)
(let loop ((k 0))
(when (< k z0-w)
(gate-cx! c z0n k z0extn k)
(loop (+ k 1))))
;; copy z2 into z2extn (low z2-w bits)
(let loop ((k 0))
(when (< k z2-w)
(gate-cx! c z2n k z2extn k)
(loop (+ k 1))))
;; z1n -= z0extn (full width 2*m_sum)
(cuccaro-sub! c z0extn z1n cin-reg cin-idx z1-w)
;; z1n -= z2extn
(cuccaro-sub! c z2extn z1n cin-reg cin-idx z1-w)
;; dst += z0 at offset 0, width 2*m_lo
(cuccaro-add-offset! c z0n 0 dst-reg dst-off
cin-reg cin-idx z0-w)
;; dst += z2 at offset 2*m_lo, width 2*m_hi
(cuccaro-add-offset! c z2n 0 dst-reg (+ dst-off (* 2 m-lo))
cin-reg cin-idx z2-w)
;; dst += z1n at offset m_lo. The middle term value = z1 - z0 - z2 =
;; a_lo*b_hi + a_hi*b_lo < 2^(n+1). So a width-(n+1) add captures
;; the full value; the extra bits of z1n (width 2*m_sum) up to bit
;; 2*m_sum-1 are guaranteed |0> and don't need gate emission. Using
;; the narrower width also keeps the write inside the dst's 2*n
;; slice (avoids writes past dst-reg[dst-off + 2n - 1]).
(cuccaro-add-offset! c z1n 0 dst-reg (+ dst-off m-lo)
cin-reg cin-idx (+ n 1))
;; restore z1n = z1 by adding z2extn then z0extn back
(cuccaro-add! c z2extn z1n cin-reg cin-idx z1-w)
(cuccaro-add! c z0extn z1n cin-reg cin-idx z1-w)
;; uncopy z2extn, z0extn
(let loop ((k 0))
(when (< k z2-w)
(gate-cx! c z2n k z2extn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k z0-w)
(gate-cx! c z0n k z0extn k)
(loop (+ k 1))))
(free! c z2extn)
(free! c z0extn))
;; ── uncompute z1, sum_b, sum_a, z2, z0 ──
(karatsuba-wide-inverse! c san 0 sbn 0 m-sum
z1n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
;; uncompute sum_a, sum_b — mirror the padded build with sub
(let ((ahpn (kara-name-w 'kara-ahp depth n))
(bhpn (kara-name-w 'kara-bhp depth n)))
(alloc! c ahpn m-sum)
(alloc! c bhpn m-sum)
;; reapply the pads
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
;; sub padded operands (inverse of add)
(cuccaro-sub! c bhpn sbn cin-reg cin-idx m-sum)
(cuccaro-sub! c ahpn san cin-reg cin-idx m-sum)
;; uncopy a_lo from san, b_lo from sbn
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c b-reg (+ b-off k) sbn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c a-reg (+ a-off k) san k)
(loop (+ k 1))))
;; uncompute the pads
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(free! c bhpn)
(free! c ahpn))
(karatsuba-wide-inverse! c a-reg (+ a-off m-lo) b-reg (+ b-off m-lo) m-hi
z2n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(karatsuba-wide-inverse! c a-reg a-off b-reg b-off m-lo
z0n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(free! c z1n)
(free! c sbn)
(free! c san)
(free! c z2n)
(free! c z0n)))
(define (karatsuba-step-inverse! c a-reg a-off b-reg b-off n
dst-reg dst-off
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
depth)
"Inverse of karatsuba-step!: subtract a*b from dst, restoring |0>.
Walks the forward gate-sequence in REVERSE order with each gate's
inverse — the standard gate-level inverse pattern.
Why this works (and why naive anti-combine fails): the forward
combine adds z0 → z2 → z1n_middle into dst at narrow widths.
Each add carries cleanly because dst grows monotonically. The
FORWARD sequence visits intermediate dst states:
0 → z0 → z0 + z2*2^(2*m_lo) → z0 + z2*2^(2*m_lo) + z1m*2^m_lo
The REVERSE-with-inversion runs subs in REVERSE order, walking
dst back through the SAME intermediate states. Each sub mirrors a
forward add → it cannot underflow because the corresponding slice
already holds (the post-add value), and subtracting recovers the
pre-add value cleanly.
In contrast, anti-combine subs in FORWARD order (z0 first) tries
to subtract z0 from dst-low-slice that holds (a*b mod 2^(2*m_lo))
which is NOT z0 — so the slice may underflow, dropping a borrow."
(let* ((m-lo (quotient (+ n 1) 2))
(m-hi (- n m-lo))
(m-sum (+ m-lo 1))
(z0-w (* 2 m-lo))
(z2-w (* 2 m-hi))
(z1-w (* 2 m-sum))
(z0n (kara-name-w 'kara-z0 depth n))
(z2n (kara-name-w 'kara-z2 depth n))
(san (kara-name-w 'kara-sa depth n))
(sbn (kara-name-w 'kara-sb depth n))
(z1n (kara-name-w 'kara-z1 depth n))
(next-depth (+ depth 1)))
(alloc! c z0n z0-w)
(alloc! c z2n z2-w)
(alloc! c san m-sum)
(alloc! c sbn m-sum)
(alloc! c z1n z1-w)
;; FORWARD: compute z0, z2, sums, z1 — same as kara-step! up to combine
(karatsuba-wide! c a-reg a-off b-reg b-off m-lo
z0n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(karatsuba-wide! c a-reg (+ a-off m-lo) b-reg (+ b-off m-lo) m-hi
z2n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(let ((ahpn (kara-name-w 'kara-ahp depth n))
(bhpn (kara-name-w 'kara-bhp depth n)))
(alloc! c ahpn m-sum)
(alloc! c bhpn m-sum)
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c a-reg (+ a-off k) san k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(cuccaro-add! c ahpn san cin-reg cin-idx m-sum)
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c b-reg (+ b-off k) sbn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(cuccaro-add! c bhpn sbn cin-reg cin-idx m-sum)
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(free! c bhpn)
(free! c ahpn))
(karatsuba-wide! c san 0 sbn 0 m-sum
z1n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
;; ── REVERSE-ORDER ANTI-COMBINE — walk forward-combine gates in
;; reverse with each gate inverted. dst transit through forward
;; intermediate states in reverse, so no underflow occurs at any
;; step.
;;
;; Forward combine order:
;; alloc z0extn, z2extn
;; copy z0 → z0extn, z2 → z2extn
;; z1n -= z0extn, z1n -= z2extn
;; dst += z0 (offset 0, width z0-w)
;; dst += z2 (offset 2*m_lo, width z2-w)
;; dst += z1n (offset m_lo, width n+1)
;; z1n += z2extn, z1n += z0extn
;; uncopy z2 ← z2extn, z0 ← z0extn
;; free z0extn, z2extn
;;
;; Reverse-with-inversion:
;; alloc z0extn, z2extn ← reverse of free
;; copy z0 → z0extn, z2 → z2extn ← reverse of uncopy
;; z1n -= z0extn, z1n -= z2extn ← reverse of += pair
;; dst -= z1n (offset m_lo width n+1)
;; dst -= z2 (offset 2*m_lo width z2-w)
;; dst -= z0 (offset 0 width z0-w)
;; z1n += z2extn, z1n += z0extn
;; uncopy z2extn, z0extn
;; free z0extn, z2extn
(let ((z0extn (kara-name-w 'kara-z0ext depth n))
(z2extn (kara-name-w 'kara-z2ext depth n)))
(alloc! c z0extn z1-w)
(alloc! c z2extn z1-w)
(let loop ((k 0))
(when (< k z0-w)
(gate-cx! c z0n k z0extn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k z2-w)
(gate-cx! c z2n k z2extn k)
(loop (+ k 1))))
(cuccaro-sub! c z0extn z1n cin-reg cin-idx z1-w)
(cuccaro-sub! c z2extn z1n cin-reg cin-idx z1-w)
;; ── REVERSED dst subs (z1n_middle first, then z2, then z0) ──
(cuccaro-sub-offset! c z1n 0 dst-reg (+ dst-off m-lo)
cin-reg cin-idx (+ n 1))
(cuccaro-sub-offset! c z2n 0 dst-reg (+ dst-off (* 2 m-lo))
cin-reg cin-idx z2-w)
(cuccaro-sub-offset! c z0n 0 dst-reg dst-off
cin-reg cin-idx z0-w)
(cuccaro-add! c z2extn z1n cin-reg cin-idx z1-w)
(cuccaro-add! c z0extn z1n cin-reg cin-idx z1-w)
(let loop ((k 0))
(when (< k z2-w)
(gate-cx! c z2n k z2extn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k z0-w)
(gate-cx! c z0n k z0extn k)
(loop (+ k 1))))
(free! c z2extn)
(free! c z0extn))
;; Uncompute z1, sums, z2, z0 — same as kara-step! tail.
(karatsuba-wide-inverse! c san 0 sbn 0 m-sum
z1n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(let ((ahpn (kara-name-w 'kara-ahp depth n))
(bhpn (kara-name-w 'kara-bhp depth n)))
(alloc! c ahpn m-sum)
(alloc! c bhpn m-sum)
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(cuccaro-sub! c bhpn sbn cin-reg cin-idx m-sum)
(cuccaro-sub! c ahpn san cin-reg cin-idx m-sum)
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c b-reg (+ b-off k) sbn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-lo)
(gate-cx! c a-reg (+ a-off k) san k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c b-reg (+ b-off m-lo k) bhpn k)
(loop (+ k 1))))
(let loop ((k 0))
(when (< k m-hi)
(gate-cx! c a-reg (+ a-off m-lo k) ahpn k)
(loop (+ k 1))))
(free! c bhpn)
(free! c ahpn))
(karatsuba-wide-inverse! c a-reg (+ a-off m-lo) b-reg (+ b-off m-lo) m-hi
z2n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(karatsuba-wide-inverse! c a-reg a-off b-reg b-off m-lo
z0n 0
sb-low-reg sb-low-idx sb-xext-reg
cin-reg cin-idx sb-const-tmp-reg
next-depth)
(free! c z1n)
(free! c sbn)
(free! c san)
(free! c z2n)
(free! c z0n)))
;;; ── mod-mul-karatsuba! — wrap Karatsuba wide-mul + Solinas reduce ─
;;;
;;; Same calling convention as mod-mul! and mod-mul-solinas!. When
;;; n <= *karatsuba-threshold*, dispatches verbatim to mod-mul-solinas!
;;; (byte-identical output). Otherwise runs Karatsuba for Stage 1, the
;;; Solinas reducer for Stage 2, Karatsuba-inverse for Stage 3.
(define (mod-mul-karatsuba! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)
"out := (a * b) mod p via Karatsuba wide-multiply + Solinas reduce.
Below threshold, byte-identical to mod-mul-solinas!."
(let ((n (- n+1 1)))
(cond
((= n 0) #t)
((<= n *karatsuba-threshold*)
;; Byte-identical schoolbook+Solinas path.
(mod-mul-solinas! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg))
(else
;; ── Stage 0: alloc wide-product buffer ──
(alloc! c (quote sb-tmp-ext) (* 2 n))
;; ── Stage 1: Karatsuba wide product sb-tmp-ext := a * b ──
(karatsuba-wide! c a-reg 0 b-reg 0 n
(quote sb-tmp-ext) 0
#f #f #f
cin-reg cin-idx #f
0)
;; ── Stage 2: Solinas reduce sb-tmp-ext mod p into out ──
(alloc! c (quote sol-lo-ext) (+ n 1))
(cx-copy-slice! c (quote sb-tmp-ext) 0 (quote sol-lo-ext) n)
(mod-add! c (quote sol-lo-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx)
(cx-copy-slice! c (quote sb-tmp-ext) 0 (quote sol-lo-ext) n)
(free! c (quote sol-lo-ext))
(alloc! c (quote sol-hi-ext) (+ n 1))
(cx-copy-slice! c (quote sb-tmp-ext) n (quote sol-hi-ext) n)
(let loop ((rest c-expansion) (current-shift 0))
(cond
((null? rest)
(mod-halve-n! c (quote sol-hi-ext) n+1 p current-shift
cin-reg cin-idx tmp-reg flag-reg flag-idx))
(else
(let* ((pair (car rest))
(sign (car pair))
(shift (cdr pair))
(delta (- shift current-shift)))
(when (< delta 0)
(error "mod-mul-karatsuba: c-expansion must be ascending"
(list current-shift shift)))
(mod-double-n! c (quote sol-hi-ext) n+1 p delta
cin-reg cin-idx tmp-reg flag-reg flag-idx)
(cond
((= sign 1)
(mod-add! c (quote sol-hi-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx))
((= sign -1)
(mod-sub! c (quote sol-hi-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx))
(else
(error "mod-mul-karatsuba: sign must be ±1" sign)))
(loop (cdr rest) shift)))))
(cx-copy-slice! c (quote sb-tmp-ext) n (quote sol-hi-ext) n)
(free! c (quote sol-hi-ext))
;; ── Stage 3: uncompute Karatsuba wide product ──
(karatsuba-wide-inverse! c a-reg 0 b-reg 0 n
(quote sb-tmp-ext) 0
#f #f #f
cin-reg cin-idx #f
0)
(free! c (quote sb-tmp-ext))))))
;;; ── mod-mul-karatsuba-sub! — out -= a*b mod p via Karatsuba ─────
(define (mod-mul-karatsuba-sub! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)
"out := (out - a * b) mod p. Below threshold, byte-identical to
mod-mul-solinas-sub!."
(let ((n (- n+1 1)))
(cond
((= n 0) #t)
((<= n *karatsuba-threshold*)
(mod-mul-solinas-sub! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg))
(else
(alloc! c (quote sb-tmp-ext) (* 2 n))
(karatsuba-wide! c a-reg 0 b-reg 0 n
(quote sb-tmp-ext) 0
#f #f #f
cin-reg cin-idx #f
0)
;; Stage 2 with SIGN INVERTED (mod-sub instead of mod-add).
(alloc! c (quote sol-lo-ext) (+ n 1))
(cx-copy-slice! c (quote sb-tmp-ext) 0 (quote sol-lo-ext) n)
(mod-sub! c (quote sol-lo-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx)
(cx-copy-slice! c (quote sb-tmp-ext) 0 (quote sol-lo-ext) n)
(free! c (quote sol-lo-ext))
(alloc! c (quote sol-hi-ext) (+ n 1))
(cx-copy-slice! c (quote sb-tmp-ext) n (quote sol-hi-ext) n)
(let loop ((rest c-expansion) (current-shift 0))
(cond
((null? rest)
(mod-halve-n! c (quote sol-hi-ext) n+1 p current-shift
cin-reg cin-idx tmp-reg flag-reg flag-idx))
(else
(let* ((pair (car rest))
(sign (car pair))
(shift (cdr pair))
(delta (- shift current-shift)))
(when (< delta 0)
(error "mod-mul-karatsuba-sub: c-expansion must be ascending"
(list current-shift shift)))
(mod-double-n! c (quote sol-hi-ext) n+1 p delta
cin-reg cin-idx tmp-reg flag-reg flag-idx)
(cond
((= sign 1)
(mod-sub! c (quote sol-hi-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx))
((= sign -1)
(mod-add! c (quote sol-hi-ext) out-reg n+1 p
cin-reg cin-idx tmp-reg flag-reg flag-idx))
(else
(error "mod-mul-karatsuba-sub: sign must be ±1" sign)))
(loop (cdr rest) shift)))))
(cx-copy-slice! c (quote sb-tmp-ext) n (quote sol-hi-ext) n)
(free! c (quote sol-hi-ext))
(karatsuba-wide-inverse! c a-reg 0 b-reg 0 n
(quote sb-tmp-ext) 0
#f #f #f
cin-reg cin-idx #f
0)
(free! c (quote sb-tmp-ext))))))
;;; ── dispatch flag ─────────────────────────────────────────────
;;;
;;; mod-arith.lsp's mod-mul! dispatcher routes to mod-mul-solinas! when
;;; *mod-mul-use-solinas* is true. To slot Karatsuba in WITHOUT editing
;;; mod-arith.lsp (byte-identity discipline), we REDEFINE mod-mul-solinas!
;;; itself: when *mod-mul-use-karatsuba* is true, the redefined version
;;; routes to mod-mul-karatsuba!; otherwise it calls the original
;;; mod-mul-solinas-impl! (the original body, captured via name change).
;;;
;;; Effect: a caller flipping (*mod-mul-use-solinas* = #t,
;;; *mod-mul-use-karatsuba* = #t) gets Karatsuba; flipping just
;;; *mod-mul-use-solinas* = #t alone gets the original Solinas.
(define *mod-mul-use-karatsuba* #f)
;; Capture the original mod-mul-solinas! into mod-mul-solinas-impl!
;; (alias). Only happens once at load time. Subsequent set! to
;; mod-mul-solinas! redirects callers through our dispatcher.
(define mod-mul-solinas-impl! mod-mul-solinas!)
(define mod-mul-solinas-sub-impl! mod-mul-solinas-sub!)
(set! mod-mul-solinas!
(lambda (c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)
(cond
(*mod-mul-use-karatsuba*
(mod-mul-karatsuba! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg))
(else
(mod-mul-solinas-impl! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)))))
(set! mod-mul-solinas-sub!
(lambda (c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)
(cond
(*mod-mul-use-karatsuba*
(mod-mul-karatsuba-sub! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg))
(else
(mod-mul-solinas-sub-impl! c a-reg b-reg out-reg n+1 p c-expansion
cin-reg cin-idx tmp-reg flag-reg flag-idx
red-tmp-reg)))))
Reference in a new issue View git blame Copy permalink