;;; 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)))))
