ML-DSA on x86-64: NTT and NTT⁻¹ in SSE2 (sign −48%, verify −35%)

README.md

@@ -795,7 +795,7 @@ yours to keep:
 
 <td>✅</td>
 
-<td>✅</td>
+<td>✅ SSE2 NTT</td>
 
 <td>✅ SHA extensions</td>
 
@@ -811,7 +811,7 @@ yours to keep:
 
 <td>✅</td>
 
-<td>✅</td>
+<td>✅ SSE2 NTT</td>
 
 <td>✅ SHA extensions</td>
 
@@ -827,7 +827,7 @@ yours to keep:
 
 <td>✅</td>
 
-<td>✅</td>
+<td>✅ SSE2 NTT</td>
 
 <td>✅ SHA extensions</td>
 

docs/algorithms/ml-dsa-44.toml

@@ -3,3 +3,4 @@ family = "Signatures"
 specs = ["MlDsa"]
 modules = ["src/mldsa44.rs"]
 asm = ["mldsa44", "mldsa"]
+optimized = { x86_64 = "SSE2 NTT" }

docs/algorithms/ml-dsa-65.toml

@@ -3,3 +3,4 @@ family = "Signatures"
 specs = ["MlDsa"]
 modules = ["src/mldsa65.rs"]
 asm = ["mldsa65", "mldsa"]
+optimized = { x86_64 = "SSE2 NTT" }

docs/algorithms/ml-dsa-87.toml

@@ -3,3 +3,4 @@ family = "Signatures"
 specs = ["MlDsa"]
 modules = ["src/mldsa87.rs"]
 asm = ["mldsa87", "mldsa"]
+optimized = { x86_64 = "SSE2 NTT" }

lean/VerifiedGarbage/Artifacts/MlDsaArith/X86_64.lean

@@ -22,7 +22,10 @@ def artifacts : List Artifact := [
   { Spec.MlDsa.nttApi with
     target := X86_64.target
     doc := Spec.MlDsa.nttApi.doc
-      (notes := ["The function stores a table of the 256 zetas in `scratch`."])
+      (notes := ["The function computes on four coefficients at a time in SSE2 registers, with a table of \
+        the 256 zetas that it stores in `scratch`. It sets MXCSR to `0x1FBF` around its multiplications \
+        (Intel's mitigation of MXCSR-configuration-dependent timing) and loads the caller's MXCSR back \
+        before returning."])
     code := Impl.MlDsa.X86_64.Arith.ntt
     contract := Spec.MlDsa.nttContract X86_64.abi
     verified := Proof.MlDsa.X86_64.Arith.ntt_verified
@@ -31,7 +34,10 @@ def artifacts : List Artifact := [
   { Spec.MlDsa.nttInvApi with
     target := X86_64.target
     doc := Spec.MlDsa.nttInvApi.doc
-      (notes := ["The function stores a table of the 256 negated zetas in `scratch`."])
+      (notes := ["The function computes on four coefficients at a time in SSE2 registers, with a table of \
+        the 256 zetas that it stores in `scratch`. It sets MXCSR to `0x1FBF` around its multiplications \
+        (Intel's mitigation of MXCSR-configuration-dependent timing) and loads the caller's MXCSR back \
+        before returning."])
     code := Impl.MlDsa.X86_64.Arith.nttInv
     contract := Spec.MlDsa.nttInvContract X86_64.abi
     verified := Proof.MlDsa.X86_64.Arith.nttInv_verified

lean/VerifiedGarbage/Impl/MlDsa/X86_64/Arith/Common.lean

@@ -16,10 +16,7 @@ Pieces of code that the ML-DSA arithmetic functions share, for
   `⌊rax / q⌋ - 1`, so `rax` less that quotient times `q` (a second `mul`) is
   less than `2q`, and `csubQ` reduces it. `mul` is the only multiplication of
   the model, and its timing does not depend on its operands (it is on
-  Intel's DOIT list). It uses `rax`, `rdx` and `r11`;
-* `storeTab t n`: the table `t 0, …, t (n - 1)` of constants stored as
-  `u32`s at `r9` (in the working space: the code has no other memory), with
-  immediates. It uses `rax`.
+  Intel's DOIT list). It uses `rax`, `rdx` and `r11`.
 -/
 
 namespace VG.Impl.MlDsa.X86_64.Arith
@@ -47,11 +44,4 @@ def reduce : List Instr :=
   [.mov .r10 (.reg .rax), .movImm64 .r11 barrettImm, .mul .r11, .mov .rax (.reg .rdx),
     .mov .r11 (.imm qImm), .mul .r11, .alu .sub .r10 (.reg .rax)] ++ csubQ .r10 .r11
 
-/-- `t i` to `[r9 + 4i]`. -/
-def tabStep (t : Nat → Nat) (i : Nat) : List Instr :=
-  [.mov32 .rax (.imm (BitVec.ofNat 32 (t i))), .store32 (at_ .r9 (4 * i)) .rax]
-
-/-- The table `t 0, …, t (n - 1)` at `r9`. -/
-def storeTab (t : Nat → Nat) (n : Nat) : List Instr := (List.range n).flatMap (tabStep t)
-
 end VG.Impl.MlDsa.X86_64.Arith

lean/VerifiedGarbage/Impl/MlDsa/X86_64/Arith/Ntt.lean

@@ -1,97 +1,104 @@
-import VerifiedGarbage.Impl.MlDsa.X86_64.Arith.Common
+import VerifiedGarbage.Impl.MlDsa.X86_64.Arith.Vec
 import VerifiedGarbage.Spec.MlDsa
 
 /-!
 # ML-DSA on x86-64: `vg_mldsa_ntt` and `vg_mldsa_inv_ntt`
 
-`ntt(f = rdi, scratch = rsi)` and `nttInv(f = rdi, scratch = rsi)`: the
-prologue stores a table of 256 zetas to `scratch` as `u32`s (`storeTab`,
-with `r9` = `scratch`): `ζ^BitRev8(m) mod q` for `NTT`, and its negation
-`-ζ^BitRev8(m) mod q` for `NTT⁻¹` (the `z` of Algorithm 42). Then each of
-the eight layers runs its blocks, with `rsi` pointing at coefficient `j` of
-`f`, `r8` at the zeta of the block, `rdi` counting the blocks down and `rcx`
-the butterflies of a block; the zeta of the block is in `r9`.
-
-* `NTT` (Algorithm 41): the layers with `len` = 128, 64, …, 1, whose zetas
-  are consecutive, from `m = 1` up. A butterfly computes
-  `t = ζ · f[j + len] mod q` (with `reduce`, a Barrett reduction with
-  `mul`), and stores `f[j] - t` (`f[j] + q - t`, reduced with `csubQ`) to
-  `f[j + len]` and `f[j] + t` (reduced) to `f[j]`.
-* `NTT⁻¹` (Algorithm 42): the layers with `len` = 1, 2, …, 128, whose zetas
-  are consecutive from `m = 255` down. A butterfly stores `f[j] + f[j + len]`
-  (reduced) to `f[j]` and `z · (f[j] - f[j + len]) mod q` to `f[j + len]`.
-  Then every coefficient is multiplied by `8347681 = 256⁻¹ mod q` and
-  reduced.
-
-A block ends with `rsi` advanced past its upper half, so a layer ends with
-`rsi` at `f + 1024`, and moves it back. Every address and branch depends
-only on the pointers.
+`ntt(f = rdi, scratch = rsi)` and `nttInv(f = rdi, scratch = rsi)` compute
+on four coefficients of `f` at a time, in place, as doublewords of SSE
+registers (see `Vec.lean`). `scratch` holds the 256 `u32`s
+`ζ^BitRev8(m) · 2³² mod q` (`zmTab`, stored with immediates: the code has
+no other memory), and MXCSR's at bytes 768 to 775 before and after
+(`withMxcsr`, which saves the caller's MXCSR in `r11` in between).
+
+Within `withMxcsr`, the prologue stores the table and the constants. Then
+the layers, each a pass over `f` with `rdx` pointing at the coefficients it
+loads and `r8` at the zetas of the table:
+
+* `NTT` (Algorithm 41): the layers with `len` = 128, 64, 32, 16, 8 and 4
+  (`vlay`) run their blocks (counted in `rax`), each its zeta in the
+  doublewords of `xmm13` (`vzeta`), and `len / 4` times the butterflies of
+  four coefficients `w[j]` and of the four `w[j + len]` (`vbfly`, counted
+  in `rcx`). The layer with `len = 2` (`vlay2`) loads the 8 coefficients of
+  two blocks, gathers their lower and upper halves into `xmm0` and `xmm1`
+  (`punpcklqdq`, `punpckhqdq`), with the two zetas in the halves of
+  `xmm13`; the layer with `len = 1` (`vlay1`) those of four blocks, their
+  pairs gathered (with `pshufd` first), with the four zetas in the
+  doublewords of `xmm13`.
+* `NTT⁻¹` (Algorithm 42): the same layers in the opposite order, with the
+  zetas from `m = 255` down and the inverse butterflies (`vibfly`); then
+  every coefficient is multiplied by `8347681 = 256⁻¹ mod q` (`vscale`,
+  with `vmont` by `8347681 · 2³² mod q = 16382`).
+
+Every address and branch depends only on the pointers.
 -/
 
 namespace VG.Impl.MlDsa.X86_64.Arith
 
 open VG.X86_64
-
-/-- `ζ^BitRev8(m) mod q`. -/
-def zetaTab (m : Nat) : Nat := 1753 ^ Spec.MlDsa.bitRev8 m % 8380417
-
-/-- `-ζ^BitRev8(m) mod q`. -/
-def negZetaTab (m : Nat) : Nat := (8380417 - zetaTab m) % 8380417
-
-/-- A butterfly of `NTT` on `[rsi]` and `[rsi + 4len]`, with the zeta in `r9`. -/
-def bfly (len : Nat) : List Instr :=
-  [.mov32 .rax (.mem (at_ .rsi (4 * len))), .mul .r9] ++ reduce ++
-    [.mov32 .rax (.mem (at_ .rsi 0)), .mov32 .rdx (.reg .rax), .alu32 .add .rdx (.imm qImm),
-      .alu32 .sub .rdx (.reg .r10)] ++ csubQ .rdx .r11 ++
-    [.store32 (at_ .rsi (4 * len)) .rdx, .alu32 .add .rax (.reg .r10)] ++ csubQ .rax .r11 ++
-    [.store32 (at_ .rsi 0) .rax, .alu .add .rsi (.imm 4), .alu .sub .rcx (.imm 1)]
-
-/-- A butterfly of `NTT⁻¹` on `[rsi]` and `[rsi + 4len]`, with the zeta in `r9`. -/
-def bflyInv (len : Nat) : List Instr :=
-  [.mov32 .rax (.mem (at_ .rsi 0)), .mov32 .r10 (.mem (at_ .rsi (4 * len))), .mov32 .rdx (.reg .rax),
-    .alu32 .add .rdx (.reg .r10)] ++ csubQ .rdx .r11 ++
-    [.store32 (at_ .rsi 0) .rdx, .alu32 .add .rax (.imm qImm), .alu32 .sub .rax (.reg .r10)] ++
-    csubQ .rax .r11 ++ [.mul .r9] ++ reduce ++
-    [.store32 (at_ .rsi (4 * len)) .r10, .alu .add .rsi (.imm 4), .alu .sub .rcx (.imm 1)]
-
-/-- A block of `len` butterflies `b`, with the zeta at `r8`, which then moves
-by `dz` bytes (4 or -4). -/
-def nttBlk (b : List Instr) (len : Nat) (dz : BitVec 32) : Prog isa :=
-  .seq (.block [.mov32 .r9 (.mem (at_ .r8 0)), .alu .add .r8 (.imm dz), .mov32 .rcx (.imm (BitVec.ofNat 32 len))])
-    (.seq (.loop (.block b) .ne)
-      (.block [.alu .add .rsi (.imm (BitVec.ofNat 32 (4 * len))), .alu .sub .rdi (.imm 1)]))
-
-/-- A layer: its `128 / len` blocks, then `rsi` back to `f`. -/
-def nttLay (b : List Instr) (len : Nat) (dz : BitVec 32) : Prog isa :=
-  .seq (.block [.mov32 .rdi (.imm (BitVec.ofNat 32 (128 / len)))])
-    (.seq (.loop (nttBlk b len dz) .ne) (.block [.alu .sub .rsi (.imm 1024)]))
-
-/-- The layers of `NTT` with `len` in `lens`. -/
-def nttLays : List Nat → Prog isa
-  | [] => .block []
-  | len :: lens => .seq (nttLay (bfly len) len 4) (nttLays lens)
-
-/-- The layers of `NTT⁻¹` with `len` in `lens`. -/
-def nttInvLays : List Nat → Prog isa
-  | [] => .block []
-  | len :: lens => .seq (nttLay (bflyInv len) len (-4)) (nttInvLays lens)
-
-/-- The table `t` to `scratch`, and `rsi` = `f`. -/
-def nttPro (t : Nat → Nat) : List Instr :=
-  [.mov .r9 (.reg .rsi)] ++ storeTab t 256 ++ [.mov .rsi (.reg .rdi), .mov .r8 (.reg .r9)]
-
-def ntt : Prog isa :=
-  .seq (.block (nttPro zetaTab ++ [.alu .add .r8 (.imm 4)])) (nttLays [128, 64, 32, 16, 8, 4, 2, 1])
-
-/-- A coefficient times `8347681`, reduced. -/
-def scaleBody : List Instr :=
-  [.mov32 .rax (.mem (at_ .rsi 0)), .mul .r9] ++ reduce ++
-    [.store32 (at_ .rsi 0) .r10, .alu .add .rsi (.imm 4), .alu .sub .rcx (.imm 1)]
-
-def nttInv : Prog isa :=
-  .seq (.block (nttPro negZetaTab ++ [.alu .add .r8 (.imm (4 * 255))]))
-    (.seq (nttInvLays [1, 2, 4, 8, 16, 32, 64, 128])
-      (.seq (.block [.mov32 .r9 (.imm 8347681)])
-        (.seq (.block [.mov32 .rcx (.imm 256)]) (.loop (.block scaleBody) .ne))))
+open VG.Impl.MlKem.X86_64 (xb xmov withMxcsr rcxLoop)
+
+/-- `ζ^BitRev8(m) · 2³² mod q`. -/
+def zmTab (m : Nat) : Nat := 1753 ^ Spec.MlDsa.bitRev8 m * 2 ^ 32 % 8380417
+
+/-- `d ← r + off`. -/
+def leaR (d r : Reg) (off : Nat) : List Instr := [.mov d (.reg r), .alu .add d (.imm (BitVec.ofNat 32 off))]
+
+/-- The zetas at `[r8]` in the doublewords of `xmm13`, arranged by `pshufd`
+with `o`, and its odd doublewords in the even ones of `xmm12`. -/
+def vzeta (o : BitVec 8) : List Instr :=
+  [.movdquLoad .xmm13 (at_ .r8 0), .xop (.pshufd .xmm13 .xmm13 o), .xop (.pshufd .xmm12 .xmm13 0xF5)]
+
+/-- A layer with `len ≥ 4` and butterflies `bf`: its `128 / len` blocks, the
+first with the zeta `k`, the zeta pointer moving by `dz` bytes. -/
+def vlay (bf : List Instr) (len k : Nat) (dz : BitVec 32) : Prog isa :=
+  .seq (.block ([.mov .rdx (.reg .rdi)] ++ leaR .r8 .rsi (4 * k) ++
+    [.mov32 .rax (.imm (BitVec.ofNat 32 (128 / len)))])) <|
+  .loop (.seq (.block (vzeta 0 ++ [.alu .add .r8 (.imm dz)]))
+    (.seq (rcxLoop (len / 4) ([.movdquLoad .xmm0 (at_ .rdx 0), .movdquLoad .xmm1 (at_ .rdx (4 * len))] ++
+        bf ++ [.movdquStore (at_ .rdx 0) .xmm0, .movdquStore (at_ .rdx (4 * len)) .xmm3,
+          .alu .add .rdx (.imm 16)]))
+      (.block [.alu .add .rdx (.imm (BitVec.ofNat 32 (4 * len))), .alu .sub .rax (.imm 1)]))) .ne
+
+/-- The layer with `len = 2`, two blocks at a time: the zetas at `[r8]`
+arranged by `pshufd` with `o`, the zeta pointer moving by `dz` bytes. -/
+def vlay2 (bf : List Instr) (k : Nat) (o : BitVec 8) (dz : BitVec 32) : Prog isa :=
+  .seq (.block ([.mov .rdx (.reg .rdi)] ++ leaR .r8 .rsi (4 * k))) <|
+  rcxLoop 32 ([.movdquLoad .xmm0 (at_ .rdx 0), .movdquLoad .xmm1 (at_ .rdx 16)] ++ vzeta o ++
+    [.alu .add .r8 (.imm dz), xmov .xmm2 .xmm0, xb .punpcklqdq .xmm0 .xmm1, xb .punpckhqdq .xmm2 .xmm1,
+      xmov .xmm1 .xmm2] ++ bf ++
+    [xmov .xmm1 .xmm0, xb .punpcklqdq .xmm0 .xmm3, xb .punpckhqdq .xmm1 .xmm3,
+      .movdquStore (at_ .rdx 0) .xmm0, .movdquStore (at_ .rdx 16) .xmm1, .alu .add .rdx (.imm 32)])
+
+/-- The layer with `len = 1`, four blocks at a time: the zetas at `[r8]`
+arranged by `pshufd` with `o`, the zeta pointer moving by `dz` bytes. -/
+def vlay1 (bf : List Instr) (k : Nat) (o : BitVec 8) (dz : BitVec 32) : Prog isa :=
+  .seq (.block ([.mov .rdx (.reg .rdi)] ++ leaR .r8 .rsi (4 * k))) <|
+  rcxLoop 32 ([.movdquLoad .xmm0 (at_ .rdx 0), .movdquLoad .xmm2 (at_ .rdx 16)] ++ vzeta o ++
+    [.alu .add .r8 (.imm dz), .xop (.pshufd .xmm0 .xmm0 0xD8), .xop (.pshufd .xmm2 .xmm2 0xD8),
+      xmov .xmm1 .xmm0, xb .punpcklqdq .xmm0 .xmm2, xb .punpckhqdq .xmm1 .xmm2] ++ bf ++
+    [xmov .xmm1 .xmm0, xb .punpckldq .xmm0 .xmm3, xb .punpckhdq .xmm1 .xmm3,
+      .movdquStore (at_ .rdx 0) .xmm0, .movdquStore (at_ .rdx 16) .xmm1, .alu .add .rdx (.imm 32)])
+
+/-- Every coefficient times `8347681 = 256⁻¹ mod q`, reduced. -/
+def vscale : Prog isa :=
+  .seq (.block [.mov .rdx (.reg .rdi), .mov32 .rax (.imm 16382), .xop (.movq .xmm13 .rax),
+      .xop (.pshufd .xmm13 .xmm13 0), xmov .xmm12 .xmm13])
+    (rcxLoop 64 ([.movdquLoad .xmm3 (at_ .rdx 0)] ++ vmont .xmm3 .xmm13 .xmm12 .xmm2 .xmm4 ++
+      vcsub .xmm3 .xmm2 ++ [.movdquStore (at_ .rdx 0) .xmm3, .alu .add .rdx (.imm 16)]))
+
+/-- The table and the constants. -/
+def vpro : List Instr := dwordTab zmTab 256 .rsi ++ vconsts
+
+def ntt : Prog isa := withMxcsr .rsi 768 <|
+  .seq (.block vpro) (.seq (vlay vbfly 128 1 4) (.seq (vlay vbfly 64 2 4) (.seq (vlay vbfly 32 4 4)
+    (.seq (vlay vbfly 16 8 4) (.seq (vlay vbfly 8 16 4) (.seq (vlay vbfly 4 32 4)
+      (.seq (vlay2 vbfly 64 0x50 8) (vlay1 vbfly 128 0xE4 16))))))))
+
+def nttInv : Prog isa := withMxcsr .rsi 768 <|
+  .seq (.block vpro) (.seq (vlay1 vibfly 252 0x1B (-16)) (.seq (vlay2 vibfly 126 0x05 (-8))
+    (.seq (vlay vibfly 4 63 (-4)) (.seq (vlay vibfly 8 31 (-4)) (.seq (vlay vibfly 16 15 (-4))
+      (.seq (vlay vibfly 32 7 (-4)) (.seq (vlay vibfly 64 3 (-4)) (.seq (vlay vibfly 128 1 (-4))
+        vscale))))))))
 
 end VG.Impl.MlDsa.X86_64.Arith

lean/VerifiedGarbage/Impl/MlDsa/X86_64/Arith/Vec.lean

@@ -0,0 +1,88 @@
+import VerifiedGarbage.Impl.MlDsa.X86_64.Arith.Common
+import VerifiedGarbage.Impl.MlKem.X86_64.Vec
+
+/-!
+# ML-DSA on x86-64: arithmetic modulo `q` in the doublewords of SSE registers
+
+The NTT and its inverse compute on four coefficients at a time, as the
+doublewords of SSE2 registers, with `q` in the doublewords of `xmm15` and
+`-q⁻¹ mod 2³² = 4236238847` in those of `xmm14` (`vconsts`).
+
+* `vmont d z zo t u`: `d ← d · z · 2⁻³² mod q`, in `[0, 2q)`, for any
+  doublewords `d` and `z < q` (a Montgomery reduction). `pmuludq` multiplies
+  the even doublewords of its operands into quadwords, so the even
+  doublewords of `d` are multiplied by those of `z`, and the odd ones,
+  moved to the even places of `u` by `pshufd`, by the even doublewords of
+  `zo`, which hold the odd doublewords of `z` (`pshufd` with `0xF5`). For
+  each product `P < 2³² · q`, `m = (P mod 2³²) · (-q⁻¹) mod 2³²` makes
+  `P + m · q` a multiple of `2³²` (with `pmuludq` by `xmm14`, then by
+  `xmm15`, which use the low doubleword of the product), less than `2⁶⁴`,
+  and its high doubleword is `(P + m · q) / 2³² < 2q`, congruent to
+  `P · 2⁻³²` modulo `q`. The quotients of the even doublewords are moved
+  down to their places by `psrlq`; those of the odd ones are in place, and
+  the low doublewords of their quadwords are 0, so `por` merges them.
+* `vcadd d t`: `d ← d + q` for the doublewords of `d` that are negative
+  (with `psrad` by 31, a mask), from `(-q, q)` to `[0, q)`; `vcsub d t`:
+  `d ← d - q`, then `vcadd`, from `[0, 2q)` to `[0, q)`.
+
+A coefficient `x` is multiplied by `ζ` as `vmont` with `ζ · 2³² mod q`,
+which the tables hold.
+
+`pmuludq` has data-dependent timing on processors with MCDT unless MXCSR is
+`0x1FBF` (see `TCB/X86_64/Isa.lean`): the functions run inside ML-KEM's
+`withMxcsr`.
+-/
+
+namespace VG.Impl.MlDsa.X86_64.Arith
+
+open VG.X86_64
+open VG.Impl.MlKem.X86_64 (xb xmov)
+
+/-- `q` in the doublewords of `xmm15` and `-q⁻¹ mod 2³² = 4236238847` in
+those of `xmm14`, through `rax`. -/
+def vconsts : List Instr :=
+  [.mov32 .rax (.imm 8380417), .xop (.movq .xmm15 .rax), .xop (.pshufd .xmm15 .xmm15 0),
+    .mov32 .rax (.imm 4236238847), .xop (.movq .xmm14 .rax), .xop (.pshufd .xmm14 .xmm14 0)]
+
+/-- The Montgomery reductions of the quadword products in `d`, with a
+temporary `t`: each quadword becomes `P + m · q`. -/
+def vredc (d t : XReg) : List Instr :=
+  [xmov t d, xb .pmuludq t .xmm14, xb .pmuludq t .xmm15, xb .paddq d t]
+
+/-- `d ← d · z · 2⁻³² mod q`, in `[0, 2q)`, with the odd doublewords of `z`
+in the even doublewords of `zo`, and temporaries `t` and `u`. -/
+def vmont (d z zo t u : XReg) : List Instr :=
+  [.xop (.pshufd u d 0xF5), xb .pmuludq d z, xb .pmuludq u zo] ++ vredc d t ++
+    [.xop (.shift .psrlq d 32)] ++ vredc u t ++ [xb .por d u]
+
+/-- `d ← d + q` for the negative doublewords of `d`, with a temporary `t`. -/
+def vcadd (d t : XReg) : List Instr :=
+  [xmov t d, .xop (.shift .psrad t 31), xb .pand t .xmm15, xb .paddd d t]
+
+/-- `d ← d mod q` for doublewords in `[0, 2q)`, with a temporary `t`. -/
+def vcsub (d t : XReg) : List Instr := xb .psubd d .xmm15 :: vcadd d t
+
+/-- The butterflies of Algorithm 41 on the doublewords of `xmm0` (`w[j]`)
+and `xmm1` (`w[j + len]`) with the zetas `ζ · 2³² mod q` in `xmm13` (and
+its odd doublewords in the even ones of `xmm12`): `xmm0 ← xmm0 + ζ · xmm1`
+and `xmm3 ← xmm0 - ζ · xmm1`. -/
+def vbfly : List Instr :=
+  vmont .xmm1 .xmm13 .xmm12 .xmm2 .xmm4 ++ vcsub .xmm1 .xmm2 ++
+    (xmov .xmm3 .xmm0 :: xb .paddd .xmm0 .xmm1 :: vcsub .xmm0 .xmm2) ++
+    (xb .psubd .xmm3 .xmm1 :: vcadd .xmm3 .xmm2)
+
+/-- The butterflies of Algorithm 42 on the doublewords of `xmm0` (`w[j]`)
+and `xmm1` (`w[j + len]`) with the zetas `ζ · 2³² mod q` in `xmm13` (and
+`xmm12`): `xmm0 ← xmm0 + xmm1` and `xmm3 ← ζ · (xmm1 - xmm0)` (Algorithm
+42 multiplies `w[j] - w[j + len]` by `-ζ`), from `xmm1 - xmm0 + q`. -/
+def vibfly : List Instr :=
+  (xmov .xmm3 .xmm1 :: xb .psubd .xmm3 .xmm0 :: xb .paddd .xmm3 .xmm15 :: xb .paddd .xmm0 .xmm1 ::
+    vcsub .xmm0 .xmm2) ++ vmont .xmm3 .xmm13 .xmm12 .xmm2 .xmm4 ++ vcsub .xmm3 .xmm2
+
+/-- The `u64`s `t (2i) + 2³² · t (2i + 1)` for `i < n / 2` at `[r + 8i]`,
+through `r9`. -/
+def dwordTab (t : Nat → Nat) (n : Nat) (r : Reg) : List Instr :=
+  (List.range (n / 2)).flatMap fun i =>
+    [.movImm64 .r9 (BitVec.ofNat 64 (t (2 * i) + 2 ^ 32 * t (2 * i + 1))), .store (at_ r (8 * i)) .r9]
+
+end VG.Impl.MlDsa.X86_64.Arith