# thirteen: double-precision math from 4-bit numbers, about 8× faster than the GPU's own FP64 --- # Part 1: For you ## The idea in one line A consumer GPU is terrible at double precision (FP64) but monstrous at 4-bit floats (FP4). Ozaki Scheme II turns one FP64 matrix multiply into 75 FP4 matrix multiplies whose integer results are exact, then glues them back together with the Chinese Remainder Theorem (CRT). The answer is more accurate than the GPU's own FP64 matmul (cuBLAS DGEMM) and about 8× faster. The trick that makes 4 bits enough is the number 13. You write it by hand in CUDA and inline PTX on a rented RTX 5090. Then a second version runs in the visitor's browser, which has no FP64 type at all, by using WebGPU's 8-bit integer dot product. ## Why it's exciting - **The hook is strange and true.** An August 2026 paper, [DGEMM with Ozaki Scheme I/II on FP4 Tensor Cores](https://arxiv.org/abs/2608.06812), emulates FP64 matrix multiplication on an RTX PRO 6000 Blackwell. At N = 16384, its FP4 Ozaki II reaches 15.26 equivalent TFLOPS end to end, against 1.85 TFLOPS for cuBLAS DGEMM. That's 8.2× faster, and it's also 4-7 bits more accurate at every dynamic range it tests. - **The number 13 does the work.** FP4 (E2M1) holds only 15 distinct values. Doubled, the set {0, ±1, ±2, ±3, ±4, ±6, ±8, ±12} hits all 13 residue classes modulo 13, so every integer has a base-13 expansion whose digits each fit in one FP4 number. - **It's exact on hardware that doesn't add like IEEE.** The [MMA-Sim paper](https://arxiv.org/abs/2511.10909) shows one dot product giving six answers across GPUs, because tensor cores align and truncate partial sums. Ozaki stays exact because every partial sum is a small integer that survives the truncation. You probe your own 5090. - **The browser has no double type, and it multiplies doubles anyway.** WGSL has only f32 and an optional f16, so the in-tab version works in 8-bit integers and only JavaScript sees FP64. - **Nobody has built this.** The only public code is the authors' Triton reference and GEMMul8. Searches found no hand-written PTX version of FP4 Ozaki II, no RTX 5090 numbers, and no browser demo of Ozaki emulation. - **It's low-level.** You write inline PTX `mma.sync` with block scales, `ldmatrix`, a `cp.async` pipeline in 99 KB of shared memory, and 160-bit software integers for the CRT. - **It has an honest twist.** On this GPU generation, INT8 emulation is still about 2.4× faster than FP4. The authors' case for FP4 is future GPUs that cut back INT8. The page says so. ## What the demo looks like - **A hook panel:** "Your browser's GPU has no double type. Watch it multiply doubles anyway." A live 1024 × 1024 (or device-appropriate) DGEMM emulation runs with timing, next to a naive f32 WGSL matmul, with effective-bits readouts on a sampled block. - **A limb explorer.** Type any double and see it scaled to an integer, its base-13 digits, its 19 residues, and each residue's two FP4 limbs, with the FP4 bit patterns. - **A CRT reassembler** that animates 19 residues combining into one 123-bit integer. - **A six-answers panel:** the MMA-Sim probe vector, that paper's answers per architecture, and what your RTX 5090 returned, with the alignment and truncation shown step by step. - **A 5090 results dashboard:** TFLOPS against N, effective bits against dynamic range, a time breakdown, and your hand-PTX kernel against Triton, cuBLAS DGEMM, and GEMMul8. - **A dynamic-range slider (φ)** that shows accuracy degrading, the authors' stated limit. - **A footer** with credits and "not affiliated with the authors." ## How it works 1. **Measure the card.** On a rented RTX 5090, measure cuBLAS DGEMM, cuBLAS SGEMM, and a hand-written PTX FP4 `mma.sync` peak. Run the six-answers probe through the tensor cores. 2. **Write the math in Python.** On your Mac, build an exact Ozaki II reference with Python integers: base-13 limbs, 19 residues, 75 exact integer matmuls, and a big-integer CRT. Check it against exact fractions, and reproduce the paper's accuracy curve. 3. **Run the authors' code as the oracle.** Run their Triton kernels and GEMMul8 on the 5090 and record the first RTX 5090 numbers for the paper's Table I. 4. **Write the kernels by hand.** In CUDA C++ with inline PTX, write the scale, encode, fused residue GEMM, and CRT kernels, bit-identical to the Python reference and the Triton oracle. 5. **Go beyond the paper.** Add rectangular sizes, masking, and moduli tiling so N = 16384 fits in 32 GB. Measure a Karatsuba variant, and sweep accuracy over φ. 6. **Port the idea to the browser.** Write an INT8 Ozaki II in WGSL with `dot4I8Packed`, check it bit for bit against a JavaScript BigInt reference, and ship the page on `vm.ifkash.dev`. ## Weekend plan | When | What | Done when | |---|---|---| | Saturday morning | Python Ozaki II reference and exact checks on the Mac; start the pod; microbenchmarks (DGEMM, FP4 `mma` peak, and the six-answers probe) | The FP4 peak and DGEMM are measured, and the reference is bit-exact on small matrices | | Saturday afternoon | Run the Triton oracle and GEMMul8 on the 5090 (Table I at 4096 and 8192); start the hand-written encode and FP4 `mma` kernels | The 5090 Table I rows are recorded, and one modulus's product bit-matches the reference | | Saturday evening | Fused residue GEMM and the CRT kernel | The full DGEMM is bit-identical to the reference at N = 4096 | | Sunday morning | Tuning, rectangular sizes and masking, moduli tiling for 16384, and the accuracy sweep; the WGSL INT8 half on the Mac | The kernel beats cuBLAS DGEMM by 5× or more at 8192 (stretch: match the Triton oracle), and the browser result matches the BigInt reference | | Sunday afternoon | Build the demo page, deploy it, record the video, and write the post | The link works | ## Cost - **GPU:** 1× RTX 5090 on RunPod Community cloud, $0.69 per hour ($0.99 on Secure cloud). The fallback is an RTX PRO 6000 96 GB at $1.69 per hour, the paper's own card. - **Time on the pod:** about 15-25 hours, by estimate. You edit on the Mac and run on the pod, which stops between jobs. - **Total:** about $10-17, by estimate, with a hard stop at $30 (about 43 Community 5090 hours). - **Hosting:** free. It's a static page on your VM; visitors' GPUs run the in-tab version. ## What you have at the end - A live page where a browser with no FP64 multiplies doubles more accurately than f32, and explains why 13. - Your own RTX 5090 Table I (the first published 5090 numbers) and an accuracy-against-φ curve. - A hand-written PTX FP4 Ozaki II DGEMM, bit-identical to the reference. - The six-answers probe of your 5090's tensor cores, and a video of the page. - A blog post titled something like *"My gaming GPU does double precision 8× faster with 4-bit numbers."* Keep it honest, and use the factor you measured. ## What might go wrong | Problem | What to do | |---|---| | The GeForce FP4 rate is lower than the whitepaper's | The microbenchmarks measure it first. Report the real speedup; the accuracy story still holds. For a direct paper comparison, fall back to an RTX PRO 6000 pod. | | The PTX fragment or scale layout is wrong | Build a one-tile `mma` test with known integer inputs first. Use the NVIDIA forum thread and CUTLASS example 79 as references, and compare with the PTX that Triton emits. | | The hand kernel is slower than Triton | That's acceptable. The deliverable is a correct, readable PTX kernel plus the measurement; keep Triton as the "paper" row. | | INT8 still beats FP4 on this card | The paper already shows it: 36.39 against 15.26 TFLOPS at 16384. Say so on the page; FP4's case is future GPUs. | | 32 GB is too small at N = 16384 | Tile over moduli. If that fails, cap at 8192 and say why. | | Wide dynamic range kills accuracy | It's the authors' stated limitation. The φ slider shows it. | | The browser lacks the packed 8-bit dot product | Fall back to i32 multiply-add in WGSL. It gives the same bits, only slower. | | CRT bugs | Use the Python big-integer reference and property tests in which random residues round-trip. | | No 5090 is available | Use Secure cloud ($0.99 per hour) or an RTX PRO 6000. | | The paper isn't CC-licensed | Link it and don't copy its figures or text. Redraw every chart from your own data. | ## Other ideas the research turned up - **[sixanswers](https://arxiv.org/abs/2511.10909).** One dot product gives six different answers depending on the GPU, in a bit-accurate web simulator of ten tensor-core generations. It lost because it likely only confirms the paper on sm_120, and the paper is CC BY-NC-SA, so it became thirteen's warm-up panel instead. - **[fouroversix](https://arxiv.org/abs/2512.02010).** Scaling an NVFP4 block to FP4's largest value (6) makes it worse: [10, 20, 30, 40] has an MSE of 4.33 scaled to 6 and 0 scaled to 4. It lost because its MIT code already runs on sm_120, so it's a reproduction. - **[gausssum](https://arxiv.org/abs/2609.04910).** FlashAttention is secretly a fast Gaussian-kernel sum solver, 2-21× faster than the PyKeOps forward pass for D > 8 on an RTX 5090. It lost because stock SDPA already does it, so a custom kernel might not win. - **[samebits](https://arxiv.org/abs/2609.25624).** Deterministic LLM inference gets faster by removing tensor cores; unmitigated BF16 diverges across GPUs on 30.81-100% of runs. It lost on no code, no Blackwell data, and too many unknowns in Mac-to-5090 bit equality. - **popcorn-120.** Port the B200-only GPU MODE reference-kernels NVFP4 GEMV and GEMM problems to consumer Blackwell with raw `mma.sync`. It lost because it has no hook for people who don't write kernels, and the reference repository's custom license needs checking. - **[dispatch](https://arxiv.org/abs/2608.08730).** In browser LLM inference at batch 1, the dispatch count matters more than kernel quality: 24-36 µs per dispatch on Vulkan and 32-71 µs on Metal. It lost because it's an optimization story rather than a surprise. ## Reading, if you want it - [DGEMM with Ozaki Scheme I/II on FP4 Tensor Cores](https://arxiv.org/abs/2608.06812): the paper this project builds on. - [Ozaki Scheme II](https://arxiv.org/abs/2504.08009) (Ozaki, Uchino, and Imamura): the CRT-based emulation scheme itself, on INT8. - [GEMMul8](https://github.com/RIKEN-RCCS/GEMMul8) and [Oz-FP4](https://github.com/FP4-is-All-you-Need/Oz-FP4): the INT8 and FP8 Ozaki II library that the paper compares against, and the authors' Triton reference. - [CUTLASS tutorial: NVFP4 block-scaled GEMM on RTX PRO Blackwell GPUs (SM12x)](https://research.colfax-intl.com/cutlass-tutorial-nvfp4-blockscaled-gemm-on-nvidia-rtx-pro-blackwell-gpus-sm12x/) (Colfax): the FP4 `mma` shape and performance on this architecture. - [PTX ISA](https://docs.nvidia.com/cuda/parallel-thread-execution/): the `mma` section, with the block-scaled variants. - [Bit-Accurate Modeling of GPU Matrix Multiply-Accumulate Units](https://arxiv.org/abs/2511.10909) (Xie et al.): MMA-Sim and the six answers. - [NVIDIA RTX Blackwell GPU architecture whitepaper](https://images.nvidia.com/aem-dam/Solutions/geforce/blackwell/nvidia-rtx-blackwell-gpu-architecture.pdf): the RTX 5090's peak rates. - [WGSL specification](https://www.w3.org/TR/WGSL/) (W3C) and [What's new in WebGPU (Chrome 123)](https://developer.chrome.com/blog/new-in-webgpu-123): the shader language and the packed 8-bit dot product. --- # Part 2: For the coding agent ## Mission Build `thirteen`, an FP64 matrix multiply emulated with FP4 tensor cores on an RTX 5090, and a browser demo of the same idea in INT8 that runs entirely on the client: - **Reference:** An exact Python Ozaki Scheme II on FP4 limbs (arXiv 2608.06812), with Python integers and a big-integer CRT, checked against exact fractions. Ozaki I is explainer-only. - **Measurements:** RTX 5090 microbenchmarks (cuBLAS DGEMM and SGEMM, the FP4 `mma.sync` peak, and the MMA-Sim probe), plus the authors' Triton kernels and GEMMul8 as the oracle. - **Kernels:** Hand-written CUDA C++ with inline PTX for the scale, encode, fused residue GEMM, and CRT steps, bit-identical to the reference, with no Triton or CUTLASS in the hot path. Beyond the paper: rectangular sizes, masking, moduli tiling for N = 16384 in 32 GB, a Karatsuba experiment, and an accuracy sweep over φ. - **Browser:** INT8 Ozaki II in WGSL with `dot4I8Packed` and a u32-word CRT, with JavaScript FP64 glue, bit-exact against a BigInt reference, inside a seven-section demo page. The final artifact is a static website with no inference server, plus the CUDA repository and its results. The authors released a Triton reference, but you write the kernels from the paper's text. Record every choice that the paper leaves open in `NOTES.md`. Work through milestones M0-M8 in order. Each milestone has acceptance criteria. Don't start a milestone until the previous one passes, except where a milestone says it runs in parallel. After each milestone, commit your work and write a short entry in `NOTES.md` with the results and numbers. ## Hard constraints - **Secrets:** Read `RUNPOD_API_KEY` from the environment only. Never write it into any file, log, commit, or echoed command. Commit a `.env.example` that has placeholder values only, and add `.env` to `.gitignore`. The project needs no Hugging Face or Weights & Biases token. - **Budget:** The hard cap is $30 of RunPod spend; the estimate is about 15-25 pod-hours, or $10-17. Every pod runs a watchdog that stops it after `MAX_POD_HOURS` hours (default `6`). If no RTX 5090 is available, use an RTX PRO 6000 ($1.69 per hour), the paper's own card. - **Pod cleanup:** Stop the pod whenever it isn't running a job. Kernel development means long idle stretches, so edit on the Mac, `rsync`, and compile and run on the pod. At the end, terminate every pod that you created and report the total spend. - **Compute split:** The Mac runs the Python reference, the goldens, the WGSL half, the web build, and the browser tests. The pod runs CUDA builds, microbenchmarks, the Triton oracle, GEMMul8, and the sweeps. - **Correctness rule:** Compare every emulated DGEMM result against an exact reference (Python integers, `fractions`, or `mpmath`) on sampled blocks, and make the hand kernel bit-identical to the Python reference. Never report speed for a configuration that failed this check. - **Benchmark hygiene:** Use CUDA events, warm-ups, and the median of 12 or more runs. Record clocks and power from `nvidia-smi` with each result, which the paper couldn't (no NVML). Lock only what a rented pod lets you lock, and note the rest in `NOTES.md`. - **Outward actions:** Ask the user before you do any of the following: create the GitHub repository, deploy to a VM, change DNS, or post anything publicly. Use the `kashifulhaque` GitHub account (`gh auth switch -u kashifulhaque`). - **Shared VM:** `vm.ifkash.dev` runs other production apps behind one shared Caddy, which owns ports 80 and 443. Its config lives at `~/docs/caddy`. - The thirteen container lives in `~/docs/thirteen` and must not publish any ports. It joins the external Docker network `edge` with a stable alias. Check existing aliases with `docker network inspect edge` first; use `thirteen` if it's free, otherwise `thirteen-demo`. - Add the vhost only by appending to `~/docs/caddy/Caddyfile` with `>>`. Never rewrite, rename, or replace that file. It's a single-file bind mount, and a rewrite orphans the inode, so `caddy reload` then reports "config is unchanged" while serving the old config. Reload from stdin, as M8 shows. - Docker on the VM has no BuildKit for plain `docker build`. Don't use `COPY --chmod`, Dockerfile heredocs, or `RUN --mount`. Multi-stage builds and `COPY --from` work. - Don't stop, restart, reconfigure, or remove any other container, network, or vhost. - **Licenses:** The paper has the arXiv perpetual non-exclusive license (nonexclusive-distrib/1.0), not Creative Commons: cite and link it, but don't copy its figures or text. Oz-FP4 and GEMMul8 are MIT; keep their notices in anything that you adapt. The MMA-Sim paper is CC BY-NC-SA 4.0: cite its Table 8 numbers as facts with attribution, and don't copy its text or figures; don't port its MIT code logic without the MIT notice. CUTLASS is BSD-3 if you copy any. State that thirteen isn't affiliated with the authors. - **Credit:** The page's first paragraph credits Ozaki, Uchino, and Imamura for Ozaki Scheme II, and Hayashi et al. for the FP4 base-13 limb representation. ## Science background These facts come from "DGEMM with Ozaki Scheme I/II on FP4 Tensor Cores: A Base-13 E2M1 Limb Representation" by Shun-ichiro Hayashi, Daichi Mukunoki, Tetsuya Hoshino, and Takahiro Katagiri (all Nagoya University), on arXiv on August 7, 2026 (v1 only, cs.DC): . Funding: JSPS KAKENHI JP25K24387 and JHPCN jh260065. Implement the scheme to match them: - **Setup (§VI-A):** An RTX PRO 6000 Blackwell Workstation Edition (sm_120, 188 SMs, 96 GB), with nominal dense peaks of 2000 TFLOPS FP4, 1000 TFLOPS FP8, and 1000 TOPS INT8. FP64 isn't on the datasheet; cuBLAS measured 1.85 TFLOPS. Bandwidth is 1792 GB/s nominal and 1490 GB/s measured. Software: driver 595.71.05, CUDA 12.8, PyTorch 2.11.0, and Triton 3.6.0, with GEMMul8 v3.1.0 built with CUDA 13.2 `nvcc` for `compute_120`/`sm_120`. Timing: CUDA events, 3 warm-ups, the median of 12 runs, and no host-device transfers. NVML was unavailable, so there's no power, clock, or temperature data. The paper never mentions the RTX 5090. - **Digit set (§III):** E2M1 values are {0, ±0.5, ±1, ±1.5, ±2, ±3, ±4, ±6}. Doubled, they give S = {0, ±1, ±2, ±3, ±4, ±6, ±8, ±12}. Each limb c ∈ S is stored as the FP4 value c/2, and the MMA result is multiplied by 4 at the end. Lemma 1: S hits all 13 residue classes mod 13, so every integer N = c + 13n with c ∈ S, and N = Σ 13^i c_i: base 13 with an unusual digit set, at log₂ 13 ≈ 3.70 bits per limb. - **Greedy decomposition (§III):** The reference `_limb_step` (`fused_decompose.py`) computes r = ((x mod 13) + 13) mod 13, takes c = r if r ≤ 6 and c = r − 13 otherwise, replaces c = 5 with −8 and c = −5 with 8 (±5 isn't in S), and continues with (x − c)/13. Generated limbs have |c| ≤ 8. The gap-free range for p limbs is |N| ≤ (13^p − 1)/3: 4, 56, and 732 for p = 1, 2, and 3. - **Integer conversion (§IV-B), shared by both schemes:** For each row of A, find the exponent e_i with 2^{e_i} ≤ max_j |A_ij| < 2^{e_i+1}, and shift so the row becomes integers with |x| < 2^53 (in code, shift = 52 − e). An all-zero row gets shift 0. Do the same for each column of B. The output is C_ij = C̃_ij · 2^{−(s_i + t_j)}, rounded to FP64 once (round-half-even). Here, small elements in a row with a big element get truncated, which the authors name as their main limitation. - **Ozaki I on FP4 (§IV), explainer only:** ⟨A, B⟩ = Σ_{i P/2, use x − P. Round to FP64 once. - **The Triton kernel:** `tl.dot_scaled(a, sa, "e2m1", b, sb, "e2m1", acc)` with every E8M0 scale set to 127 (2^0 = 1) at block size 32, which gives MXFP4 semantics with unit scales, and an FP32 accumulator. FP4 packing is `nib = sign << 3 | idx`, where idx 0-7 maps to 0, 0.5, 1, 1.5, 2, 3, 4, and 6, with two nibbles per byte along K and even k in the low nibble. The layout is TN, with both operands K-contiguous. The fused residue kernel uses 128 × 128 tiles, BK = 128, 8 warps, and 3 stages, the largest configuration that fits the 99 KB shared-memory limit. The paper names no PTX instruction, so the exact lowering is unverified. The paper's Table I reports equivalent DGEMM TFLOPS (2N³/t) for N × N × N, as compute-only / end-to-end: | Method | 4096 | 8192 | 16384 | |---|---|---|---| | cuBLAS DGEMM | 1.85 / 1.85 | 1.85 / 1.85 | 1.85 / 1.85 | | OzI-FP4 | 5.60 / 4.96 | 5.41 / 5.09 | 5.58 / 5.41 | | OzII-FP4 | 16.28 / 10.81 | 17.11 / 13.52 | 17.31 / 15.26 | | GEMMul8-FP8 | 13.68 / 11.64 | 15.41 / 13.88 | 15.68 / 14.65 | | GEMMul8-INT8 | 32.41 / 25.45 | 39.78 / 33.49 | 40.94 / 36.39 | - **Speed:** OzII-FP4 beats cuBLAS by 8.8×, 9.2×, and 9.4× compute-only and 5.8×, 7.3×, and 8.2× end to end, and GEMMul8-FP8 by 1.10-1.19× in compute. The model (Eq. 4) is t ≈ 2QN³/R + D/β, with R = 1408 TFLOPS and β = 1490 GB/s; measured compute is within 2-5% of it, at 61-65% of the FP4 peak. Take the definitions of Q and D from the paper. - **At N = 8192:** scale 9.26 ms, encode 7.84 ms, fused residue GEMM 55.63 ms, CRT 8.63 ms, total 81.35 ms. In the ablation, Ozaki II goes from 249.96 ms (Garner CRT) to 183.32 ms (direct CRT) to 81.35 ms (fused), and Ozaki I from 625.83 ms to 216.20 ms. - **Accuracy (§VII-A, Fig. 6):** Inputs are ±(1 + u) · 2^e with u ~ U[0, 1) and e ~ U{0..φ}, for φ ∈ {0, 1, 2, 4, 8, 16, 32}, with K = 8192. All elements of a 128 × 128 output are compared against 220-bit reference arithmetic. The metric is ε = max |err| / (|A||B|), elementwise, and effective bits = −log₂ ε. Ozaki I and II are bit-identical. Speed benchmarks use φ = 4. The paper's effective bits by φ are as follows: | φ | FP4 Ozaki | GEMMul8-FP8 | cuBLAS DGEMM | |---|---|---|---| | 0 | 58.17 | 57.23 | 51.44 | | 4 | 56.49 | 56.41 | 51.16 | | 8 | 55.32 | 55.58 | 50.05 | | 16 | 54.08 | 54.92 | 49.40 | | 32 | 52.26 | 53.77 | 48.00 | FP4 Ozaki beats cuBLAS DGEMM by 4-7 bits at every φ. The caveat that the page must carry: on this GPU generation, INT8 emulation (GEMMul8-INT8) is still about 2.4× faster than FP4 (36.39 against 15.26 TFLOPS at 16384). The authors' case for FP4 is future GPUs, and they name the B300 and Rubin, where INT8 throughput is cut back. The authors' repository, [Oz-FP4](https://github.com/FP4-is-All-you-Need/Oz-FP4) (MIT, "Copyright (c) 2026 the authors", 2 commits, last on 2026-08-11), needs sm_120, CUDA 13.2 `nvcc` for the GEMMul8 bridge, PyTorch 2.11 with cu128, Triton 3.6, and Python 3.12, with no pinned requirements file. Its kernels have no masking, so sizes must be multiples of 128 (fused) or 256. `dgemm_ii_fused_tn` and `dgemm_ii_tnk7` pass `n` as K, so only square matrices work, and there are no transposes, alpha, beta, or leading dimensions. Inf, NaN, subnormals, and exponent overflow are unhandled, and the INT8-emulation kernel directories that the code references are missing. By estimate (not the paper's), Ozaki II at N = 16384 needs roughly 45 GB, of which the int32 residue buffer, L · N², is about 20.4 GB. That doesn't fit a 32 GB RTX 5090 without tiling over N or over moduli; N = 8192 fits. The related work and hardware facts are as follows: - **Ozaki Scheme II:** Ozaki, Uchino, and Imamura, "Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique", arXiv 2504.08009 (v4, September 16, 2026). It's CRT-based and reports 7.4-9.8 TFLOPS on an RTX 4090 and 56.6-80.2 TFLOPS on a GH200, both with INT8. - **GEMMul8:** [RIKEN-RCCS/GEMMul8](https://github.com/RIKEN-RCCS/GEMMul8) (MIT), release v3.4.4 (2026-09-18). It implements Ozaki II emulation (`gemm`, `symm`, `syrk`, `trsm`, and more) on INT8 or FP8, for CUDA and HIP, with a cuBLAS hook mode, and requires CUDA 12.9 or later and C++20. The paper's harness used `NM_FP8=13`, `NM_INT8=15`, and `FASTMODE=0` ("accuracy-oriented mode"). - **MMA-Sim:** "Bit-Accurate Modeling of GPU Matrix Multiply-Accumulate Units" by Xie et al., arXiv 2511.10909 (v2, April 14, 2026, CC BY-NC-SA 4.0), with code at [microsoft/MMA-Sim](https://github.com/microsoft/MMA-Sim) (MIT). It covers ten architectures, including "RTX Blackwell (sm120)" on an RTX PRO 6000. Its Table 4 lists sm120 FP8, FP6, and FP4 MMA as FP32 accumulate with 25 fractional alignment bits and round-toward-zero. Its Table 8 vector, a = (−2¹³, −0.5, −0.25, −0.125, 0, …), b = (2¹⁰, 1, 1, 1, 0, …), c = 2²³, gives six outputs across GPUs and instructions: 0.0, −0.375, −0.5, −0.75, −0.875, and −1.0; exact is −0.875, and the B200 and RTX Blackwell return −0.75. Ozaki's exactness relies on every partial sum being a small integer that survives this truncation. Use the vector as a probe; reimplement from the paper and cite it rather than port its code logic. - **RTX 5090 peaks (RTX Blackwell whitepaper), dense:** FP4 with FP32 accumulate 1676 TFLOPS, FP8 419 TFLOPS, BF16 209.5 TFLOPS, and FP32 non-tensor 104.8 TFLOPS; 1792 GB/s; 170 SMs. FP64 isn't listed. The usual GeForce 1/64 ratio gives about 1.6 TFLOPS, an estimate; measure cuBLAS DGEMM on the pod. - **sm_120:** Warp-level `mma.sync` only, with no `tcgen05`, no TMEM, and no `wgmma` (`ptxas` rejects it). The GeForce cluster shape is effectively 1 × 1 × 1, with no multicast. Shared memory is 100 KB per SM and 99 KB per block. Block-scaled `mma` needs the architecture-specific target `sm_120a`; confirm this against the PTX ISA in M0. Triton has a native FP4 `dot_scaled` path on sm_120. - **FP4 `mma` on SM12x:** The Colfax tutorial states that "the MMA shape for NVFP4 blockscaling on SM12x is fixed as 16x8x64" and reaches about 60% of the 2000 TFLOPS peak on an RTX PRO 6000. An [NVIDIA forum thread](https://forums.developer.nvidia.com/t/fragment-layout-for-mma-sync-aligned-m16n8k64-with-fp4-e2m1-and-block-scaling-on-sm120/364020) (2026-06-23) covers the m16n8k64 FP4 fragment and scale-factor lane layout. The [CUTLASS Blackwell docs](https://docs.nvidia.com/cutlass/latest/media/docs/cpp/blackwell_functionality.html) list the sm_120 kinds `kind::f8f6f4`, `kind::mxf8f6f4.block_scale`, `kind::mxf4.block_scale`, and `kind::mxf4nvf4.block_scale` with `scale_vec::[2X|4X]`, TN only. CUTLASS example 79 (BSD-3) has an sm_120 NVFP4 GEMM; CUTLASS 4.8.0 is dated 2026-09-17. - **WebGPU:** WGSL has no f64 type, only f32 and an optional f16. It has `fn dot4I8Packed(a: u32, b: u32) -> i32` and `dot4U8Packed` under the language feature `packed_4x8_integer_dot_product` (Chrome 123); check it with `navigator.gpu.wgslLanguageFeatures.has(...)`. Without it, fall back to i32 multiply-add, which is slower and exact. Safari and Firefox support is unverified; test it in M7. ## Design decisions - **Ozaki II only in kernels:** Ozaki I is explainer-only: 225 GEMMs and about 3× slower. - **PTX path:** Use `mma.sync.aligned.m16n8k64.row.col.kind::mxf4.block_scale.scale_vec::2X` with e2m1 operands, f32 accumulate, and a ue8m0 scale of 127 (unit), to match the reference's MXFP4 semantics. Confirm the exact mnemonic and the operand and scale fragment layouts against the PTX ISA and the forum thread in M0 and M1, and record them. Measure two alternatives in M1: `kind::mxf4nvf4.block_scale.scale_vec::4X` with ue4m3 unit scales, and `kind::f8f6f4` with e2m1 operands, which probably runs at the FP8 rate. - **Equality tests:** Exactness holds by construction, so tests check equality. Use tolerances only for timing. - **Fast exact reference:** Each limb-plane product has entries of at most 144K in magnitude, which is about 2.4 × 10⁶ at K = 16,384, far below 2^53. So `ref/ozaki2.py` can compute the 75 products as ordinary `float64` NumPy (BLAS) matmuls and still be exact, then convert to Python integers for the mod-m and CRT steps. This makes full goldens at N = 4096 practical on the Mac. Assert the bound in code, and cross-check against pure-integer matmuls on small sizes. - **K limit:** Support K ≤ 16,384 exactly, as the paper does. Document that larger K needs split-K with integer summation, a stretch goal. - **Memory:** Tile over groups of moduli, for example by processing moduli in chunks and accumulating CRT partial sums, so that N = 16384 fits in 32 GB. Measure the cost. - **Digit tables and rounding:** Use the reference's rule, the first (a, b) in S order, so outputs match the Triton oracle, and record it. Round half to even at the final FP64 conversion. - **Special values:** Reject Inf and NaN inputs with an error. Don't flush subnormal inputs; handle them with the same exponent rule, test them, and record the behavior. - **Browser half in INT8:** WebGPU has no FP4 matrix hardware, so the in-tab engine uses INT8, like GEMMul8. The page says that the 5090 numbers come from the FP4 kernels and that the in-tab run is their INT8 cousin. JavaScript handles scaling and the final rounding with `Float64Array`, because JavaScript numbers are FP64; the GPU never sees an FP64 value. - **Browser moduli:** Read GEMMul8's source for its INT8 moduli, and cite it if you use them. Otherwise, pick pairwise-coprime moduli up to 256 greedily until 2K · 2^106 < P for the page's largest K. Record the set, and derive the i32 exactness bound in `NOTES.md`. - **Speed (estimates to verify):** The 5090 has 170 SMs against the PRO 6000's 188, so expect roughly 90% of the paper's Ozaki II numbers: about 12-14 end-to-end TFLOPS at N = 8192, against about 1.6 TFLOPS for cuBLAS DGEMM, or about 8×. Measure it. The hook depends on this factor, and if the 5090 lands lower, the page reports the real number. ## Tech stack Pin every version, and record it in `NOTES.md`. The stack is as follows: - **Python (Mac):** Python 3.12, `numpy`, `mpmath`, `matplotlib`, and `pytest`, managed with `uv`. Log runs to JSONL under `results/`, with no hosted tracker. - **CUDA (pod):** CUDA C++ with inline PTX, built with `nvcc -arch=sm_120a` from a `Makefile`, with cuBLAS only as a baseline. Record the toolkit and driver versions from the pod. - **Oracle (pod):** Oz-FP4 (PyTorch 2.11 with cu128, Triton 3.6, Python 3.12, and CUDA 13.2 `nvcc` for the GEMMul8 bridge) and GEMMul8 (CUDA 12.9 or later, C++20), in `third_party/`. - **Browser (Mac):** TypeScript, bundled with `vite`, and raw WebGPU with WGSL, tested with Playwright on Chromium with WebGPU enabled. Don't use a compute library. - **Pods:** The `runpod` Python SDK, with `runpodctl` inside pods. `infra/pod.py` creates, stops, and terminates pods named `thirteen-*` and prints their spend. ## Repository layout Create the following layout: ``` thirteen/ pyproject.toml .env.example .gitignore README.md NOTES.md ref/ limbs.py # S, base-13 decomposition, two-limb digit tables per modulus ozaki2.py # exact OzII with Python ints; direct CRT; FP64 rounding ozaki1.py # explainer only exact.py # fractions/mpmath reference, effective-bits metric, input generator (phi) golden.py # dump inputs, residues, packed limbs, outputs cuda/ common.cuh scale.cu encode.cu mma_fp4.cuh residue_gemm.cu crt.cu dgemm.cu bench/ # dgemm_cublas.cu fp4_peak.cu probe_table8.cu bench_ozaki.cu tests/ # parity vs goldens Makefile # -arch=sm_120a oracle/ # scripts that run Oz-FP4 Triton and GEMMul8 (third_party/, not vendored) infra/pod.py infra/watchdog.sh infra/bootstrap.sh infra/sync.sh web/ src/gpu/ # dp4a_gemm.wgsl crt.wgsl f32_gemm.wgsl engine.ts src/ozaki/ # moduli.ts bigint_ref.ts split.ts src/panels/ # hook, limbs, crt, sixanswers, results, phi src/main.ts index.html public/ (results JSON, video) tests/ deploy/Dockerfile deploy/compose.yml deploy/nginx.conf results/ post/draft.md ``` ## M0: Setup Do this milestone on the Mac, except for the pod smoke test. **Tasks:** 1. Scaffold the repository, the `uv` project, and the `vite` app. Add `.gitignore` and `.env.example`. 2. Write `infra/pod.py` (create, stop, and terminate `thirteen-*` pods, and print spend), `infra/watchdog.sh`, and `infra/bootstrap.sh` and `infra/sync.sh`, which `rsync` the tree. 3. Clone Oz-FP4 and GEMMul8 into `third_party/`. Don't commit their code: either ignore `third_party/` or add them as submodules. Pin both to a commit, and record the commits. 4. Read §III-§VII of the paper and the PTX ISA's `mma` block-scale section. Confirm that block-scaled `mma` needs `sm_120a`. 5. Write the "Not stated" list into `NOTES.md`, with your choice for each item: the digit table rule, rounding, special values, large K, the square-only bug, memory at 16384, and the exact PTX kind. **Acceptance criteria:** - `infra/pod.py` creates, stops, and terminates a test pod and prints its spend, and the watchdog stops a test pod when `MAX_POD_HOURS` is set to a small value. - `git ls-files` shows no `.env` and no third-party source. - `NOTES.md` records the two pinned commits, the `sm_120a` finding, and the "Not stated" table. ## M1: Microbenchmarks and the probe Run this milestone on the pod. Start it on Saturday morning, in parallel with M2. **Tasks:** - Write `cuda/bench/dgemm_cublas.cu`: cuBLAS DGEMM and SGEMM at N = 4096 and 8192. - Write `cuda/bench/fp4_peak.cu`: a register-resident loop of FP4 block-scaled `mma.sync.aligned.m16n8k64` with unit scales, for each candidate kind in the design decisions. Report TFLOPS against the whitepaper's 1676. - Confirm the operand and scale fragment layouts with a one-tile `mma` test on known integers. - Write `cuda/bench/probe_table8.cu`: run the MMA-Sim Table 8 vector through `mma.sync` for the FP16, BF16, TF32, FP8, and FP4 kinds, and record what the 5090 returns. MMA-Sim's sm120 model predicts −0.75 for FP8 and FP4. Some probe values fall outside the FP8 and FP4 ranges, so work out how to encode the vector for each kind, and record the encoding. **Acceptance criteria:** - `NOTES.md` has the median DGEMM, SGEMM, and FP4 peak numbers, with `nvidia-smi` clocks and power. - `NOTES.md` has a probe table (kind, the 5090's output, and MMA-Sim's reported outputs), and records the chosen PTX kind, its confirmed mnemonic and fragment layouts, and why. ## M2: Python reference Do this milestone on the Mac. **Tasks:** - Write `ref/limbs.py`: S, the greedy base-13 decomposition, and the two-limb digit tables for all 19 moduli, using the first (a, b) in S order. - Write `ref/ozaki2.py`: per-row and per-column scaling, 19 residues, the 75 products as exact integer matmuls, the direct CRT with Python big integers, the symmetric range, and one round-half-even conversion to FP64. - Write `ref/ozaki1.py` (explainer only); `ref/exact.py`, with the `fractions` and `mpmath` reference, the effective-bits metric, and the §VII-A φ input generator; and `ref/golden.py`, which dumps inputs, residues, packed limbs, and outputs for the kernel tests. **Acceptance criteria:** - On 200 random small cases (M, N ≤ 64, K ≤ 8192, φ ∈ {0, 4, 32}), the result matches the correctly rounded exact result (`fractions`) within the paper's accuracy. Report the bits. - Tests confirm Lemma 1, two-limb coverage for every residue of all 19 moduli, pairwise coprimality, and log₂ P = 123.22, and a property test round-trips random residues through the CRT. - The accuracy curve at K = 8192 with a 128 × 128 output is within 1 bit of Fig. 6 at φ = 0, 4, and 8. - Golden dumps exist for the kernel tests at N = 4096 and the M6 rectangular sizes, with full outputs from the `float64` BLAS path in the design decisions, and sampled blocks at 8192. - The `float64` BLAS path equals the pure-integer path bit for bit on sizes up to 256. ## M3: Oracle on the pod **Tasks:** - Write the `oracle/` scripts that run the Oz-FP4 Triton kernels and GEMMul8 INT8 and FP8 at N = 4096 and 8192, and at 16384 only if memory allows. Use the paper's GEMMul8 settings: `NM_FP8=13`, `NM_INT8=15`, and `FASTMODE=0`. - Dump the PTX and SASS that Triton emits for `dot_scaled` with e2m1, and record which `mma` kind it uses. **Acceptance criteria:** - The Triton output is bit-identical to the Python reference on the golden inputs. - `NOTES.md` has the 5090 Table I rows (cuBLAS DGEMM, OzII-FP4, GEMMul8-FP8, and GEMMul8-INT8, compute-only and end to end) next to the paper's, and the `mma` kind that Triton uses. ## M4: Hand kernels, part 1 **Tasks:** - Write `cuda/scale.cu` (per-row and per-column exponents, shift = 52 − e, and shift 0 for zero rows) and `cuda/encode.cu` (FP64 to integer, to 19 residues, to two FP4 limbs each, packed as nibbles in the `mma` operand layout). - Write `cuda/mma_fp4.cuh` and a single-modulus FP4 `mma` GEMM: 4 (or 3, for m = 169) products accumulated in FP32 registers and combined mod m into int32 residues. Use `ldmatrix` and the chosen `mma.sync` kind with unit scales. - Write `cuda/tests/` parity tests against the goldens. **Acceptance criteria:** - The residues and packed limbs are byte-identical to the reference, and the single-modulus int32 residues are identical to it at N = 4096. ## M5: Hand kernels, part 2 **Tasks:** - Write `cuda/residue_gemm.cu`: the fused 19-modulus residue GEMM, with a `cp.async` (or TMA) pipeline inside 99 KB of shared memory. Start from the paper's configuration (128 × 128 tiles, BK = 128, 8 warps, 3 stages), and record what you tune. - Write `cuda/crt.cu`: 160-bit software-integer accumulation over 19 residues, the FP64 quotient estimate with the ±1 fix, the symmetric range, and one rounding to FP64. - Write `cuda/dgemm.cu`, which chains the kernels, and `cuda/bench/bench_ozaki.cu`, which times compute-only and end-to-end runs with the per-step breakdown. **Acceptance criteria:** - The full DGEMM FP64 output is bit-identical to the Triton oracle at N = 4096 and 8192, and running the same input twice gives identical bits. - The kernel is 5× or more faster than cuBLAS DGEMM end to end at N = 8192. Report the result against Triton, whatever it is. - `NOTES.md` has the time breakdown (scale, encode, residue GEMM, and CRT) next to the paper's. ## M6: Beyond the paper **Tasks:** - Add rectangular M × N × K support, which fixes the square-only bug, and masking for sizes that aren't multiples of 128. - Add tiling over moduli so that N = 16384 fits in 32 GB, and measure its cost. - Implement the 68-GEMM Karatsuba variant as a measured experiment. The paper says it lost under fusion at this shared-memory size; measure it on the 5090. - Run the accuracy sweep over φ ∈ {0, 1, 2, 4, 8, 16, 32}, and build a roofline with the Eq. 4 model, with R re-measured on the 5090. **Acceptance criteria:** - Rectangular sizes, such as 1000 × 3000 × 777, are bit-identical to the reference. - N = 16384 runs in under 32 GB, or `NOTES.md` documents the cap and the reason. - `NOTES.md` has tables for the Karatsuba experiment, the φ sweep next to Fig. 6, and the roofline. ## M7: Browser half and demo page Do this milestone on the Mac. The WGSL engine can start in parallel with M4, once M2 passes. **Tasks:** - Choose the INT8 moduli set as the design decisions describe. Write `src/ozaki/moduli.ts`, `src/ozaki/split.ts`, and `src/ozaki/bigint_ref.ts`, a JavaScript BigInt reference. - Write `src/gpu/dp4a_gemm.wgsl`: int8 residues, `dot4I8Packed` GEMMs accumulating in i32, and the i32 multiply-add fallback. - Write `src/gpu/crt.wgsl`: the CRT with software big integers on u32 words. - Write `src/gpu/f32_gemm.wgsl`, the naive f32 comparison, and `engine.ts`, with the FP64 scaling and final rounding in `Float64Array`. - Build the seven page sections that "What the demo looks like" in Part 1 describes. The hook panel's effective bits come from an exact reference on a sampled block, the six-answers panel shows your M1 result, and the dashboard reads M3-M6 results from `public/` JSON. Without WebGPU, show a recorded video and the static charts. - Write Playwright tests. **Acceptance criteria:** - The WGSL result is bit-identical to the JavaScript BigInt reference for 20 random cases up to 256 × 256 × 1024. - A 1024³ run completes in Chrome on an Apple M-series Mac. Record the time. - The page works in Chrome and Safari and on one phone (record which support the packed 8-bit dot product), and the fallback path works with WebGPU disabled. - `NOTES.md` records the moduli set, its source, and the i32 exactness bound. ## M8: Deploy and write-up **Tasks:** - Write `deploy/Dockerfile` (`nginx:alpine` serving `web/dist`, buildable with the legacy builder) and `deploy/compose.yml`, which publishes no ports and joins the external network `edge` with the alias chosen in the hard constraints. - Ask the user before deploying. The user confirms a domain such as `thirteen.ifkash.dev` and adds a non-proxied Cloudflare A record. - After the user approves, copy the build and `deploy/` to `~/docs/thirteen`, and run `docker compose up -d` there. Then compare the live and host Caddyfiles with `docker exec caddy cat /etc/caddy/Caddyfile | diff - ~/docs/caddy/Caddyfile`. If they differ, stop and ask the user. Otherwise, back up, append, validate, and reload: ``` cd ~/docs/caddy cp Caddyfile "Caddyfile.bak-thirteen-$(date +%Y%m%d)" printf '\nthirteen.ifkash.dev {\n\treverse_proxy thirteen:80\n}\n' >> Caddyfile docker exec -i caddy caddy validate --config - --adapter caddyfile < Caddyfile docker exec -i caddy caddy reload --config - --adapter caddyfile < Caddyfile ``` Replace `thirteen` in `reverse_proxy` with the alias you chose. Omit any `tls` block: for a non-proxied A record, the default ACME HTTP-01 challenge works. - Write `README.md`: what the project is, one command to reproduce each milestone, the results tables, every open choice and deviation, the licenses, and the credits. - Write `post/draft.md`, a blog post of 1,200-1,800 words. Structure it as follows: - The hook: *"My gaming GPU does double precision 8× faster with 4-bit numbers,"* with the factor you measured, and why FP64 is slow on GeForce. - The math: the 13 trick, then the CRT with 19 moduli, two limbs each, and 75 GEMMs. - Exactness and the six answers, and the PTX kernel. - Results against cuBLAS DGEMM, Triton, and GEMMul8, with the INT8 caveat. - Limitations: dynamic range, the K limit, and memory at 16384. - Ask the user before you push anything; after approval, push on the `kashifulhaque` account. The owner lists the project on projects.dotslasha.me; that isn't your job. **Acceptance criteria:** - The deployed URL loads over HTTPS and runs the in-tab DGEMM in Chrome. - Every other site on the VM still responds as it did before the deploy. - No pods are left running, and `NOTES.md` records the total spend, which is less than $30. ## Final report to the user When you finish, report the following: - The 5090 microbenchmarks (cuBLAS DGEMM, the FP4 peak, and the chosen PTX kind) and the six-answers probe results, from M1. - The 5090 Table I (Triton oracle, GEMMul8, and your hand kernel) next to the paper's, from M3 and M5, and accuracy against φ next to Fig. 6, from M2 and M6. - The parity results for the hand kernel and the browser engine, from M5 and M7. - The beyond-paper results (rectangular sizes, N = 16384, and Karatsuba), from M6. - The browser timing and which browsers support the packed 8-bit dot product, from M7. - Every open choice from the "Not stated" list and what you picked. - The URL, if the site is deployed, and the total spend. - Anything that you skipped or that failed.