thirteen

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, 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 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

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

Other ideas the research turned up

• sixanswers. 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. 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. 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. 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. 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: the paper this project builds on.
• Ozaki Scheme II (Ozaki, Uchino, and Imamura): the CRT-based emulation scheme itself, on INT8.
• GEMMul8 and 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) (Colfax): the FP4 mma shape and performance on this architecture.
• PTX ISA: the mma section, with the block-scaled variants.
• Bit-Accurate Modeling of GPU Matrix Multiply-Accumulate Units (Xie et al.): MMA-Sim and the six answers.
• NVIDIA RTX Blackwell GPU architecture whitepaper: the RTX 5090's peak rates.
• WGSL specification (W3C) and What's new in WebGPU (Chrome 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): https://arxiv.org/abs/2608.06812. 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, j<q} 13^{i+j} ⟨a_i, b_j⟩ (Eq. 1). FP64 needs 15 limbs, so 225 GEMMs, grouped by g = i + j and summed with weight 13^g in software 128-bit integers (four 32-bit words). Lemma 2 needs K ≤ 2^18 / n_g, or K ≤ 17,476 for the deepest group (n_g = 15). It's about 3× slower than Ozaki II (Table I).
• Ozaki II moduli (§V): Work modulo L = 19 pairwise-coprime moduli: {169, 115, 113, 112, 111, 109, 107, 103, 101, 97, 89, 83, 79, 73, 71, 67, 61, 59, 53}. Their product P ≈ 2^123.2 (log₂ P = 123.22, checked, with pairwise coprimality; the repository's docstring says 2^124, which is wrong). Eq. 3 requires 2K · 2^106 < P, which holds for K ≤ 16,384, where the left side is at most 2^121. No 18-modulus set exceeds 2^117.80, so 19 is minimal.
• Two limbs per residue (§V): Each residue is r = a0 + 13·a1 with a0, a1 ∈ S. Every modulus up to 113 works, plus {115, 117, 143, 169}, and all 19 moduli have a two-limb representative for every residue (checked). 169 = 13² excludes other multiples of 13. The code's digit tables take the first (a, b) in S order per residue, so limbs can reach |12|.
• 75 GEMMs (§V): Per modulus, ab ≡ a0b0 + 13(a0b1 + a1b0) + 169·a1b1 (mod m), or 4 GEMMs. For m = 169 the a1b1 term vanishes, so the total is 4L − 1 = 75. Karatsuba could make 68, but under fusion it was slower (§VI-D).
• Exactness bounds (Lemma 3): Each product needs 144K ≤ 2^24 (K ≤ 116,508). The cross-term accumulator needs 288K ≤ 2^24 (K ≤ 58,254), which is the binding limit. INT32 residue composition needs 28,224K ≤ 2^31 (K ≤ 76,087), and the CRT needs K ≤ 76,549. There's no K-chunking: the whole K accumulates in one FP32 accumulator.
• Direct CRT (§V-B): Precompute 128-bit weights w_i = M_i (M_i^{−1} mod m_i), where M_i = P/m_i. Sum S = Σ r_i w_i in 160-bit software integers. Estimate ⌊S/P⌋ in FP64, and fix it with a ±1 correction (the quotient is at most 1753; 754 was observed). If x > 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:

• 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 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 (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 (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 (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 (2026-06-23) covers the m16n8k64 FP4 fragment and scale-factor lane layout. The CUTLASS Blackwell docs 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.