I’ve released a C implementation for computing G = A Aᵀ when consecutive rows differ from a one-position shift at only a few coordinates.
The part I find most interesting is the arithmetic order. The update replaces all left factors first, then all right factors. Every intermediate accumulator remains a sum of at most n products of bounded original input values.
For inputs A = 2^e X, with integer coefficients bounded by B and n B² ≤ 2^53, this gives exact Gram entries in binary64 under the paper’s exponent and IEEE arithmetic assumptions. The implementation validates the input rather than quantizing it.
On an Apple M4, using single-thread Accelerate and warm repeated inputs, one full-output 768 × 8192 case takes 9.05 ms versus 19.63 ms for SYRK. Validation, allocation and output materialization are included. Smaller cases lose: every tested 192 × 768 case was slower.
The repository includes the paper, C code, Lean formalization, integer-oracle checks and all benchmark results. The formalization covers the mathematical model; the C implementation is tested separately.
https://github.com/LeonardSEO/exact-sparse-gram
I’d be interested in workloads where sparse changes between shifted rows occur naturally. Would this input structure be useful in your work, and which structured baseline would you want compared next?