Pith. sign in

REVIEW 3 major objections 6 minor 24 references

Multiple Double Arithmetic on NVIDIA Tensor Cores

T0 review · 3 major / 6 minor · reviewed 2026-07-09 · glm-5.2

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Splitting doubles into 13-bit quarters unlocks tensor cores for extended precision

desk verdict The core idea is sound but the headline performance number is apples-to-oranges. read the letter →

arxiv 2607.06881 v1 pith:Y4ST3CVY submitted 2026-07-08 cs.MS cs.DCcs.NAmath.NA

Howard Chen , Jan Verschelde This is my paper · ORCID
classification cs.MScs.DCcs.NAmath.NA
keywords tensorcoresdoublemultiplearithmeticcorenvidiaozaki
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tackles a fundamental mismatch: multiple double arithmetic (representing numbers as unevaluated sums of doubles for extra precision) requires branching for renormalization after every operation, but tensor cores are pure matrix-multiply engines with no branching capability. The authors sidestep this by partitioning each 52-bit double fraction into four 13-bit quarters, each padded with 26 trailing zero bits. This padding guarantees that products and long inner products of quarters are computed exactly in standard double precision, because the results fit within the 52-bit fraction without overflow. A product of two double-double matrices is then rewritten as a single product of double matrices arranged in a block-convolution structure, which tensor cores can execute directly. The authors describe balancing algorithms that redistribute bits when input exponents deviate from the intended grid spacing, and report 13.75 TFLOPS on an A100 GPU for double-double matrix multiplication, exceeding the 9.7 TFLOPS theoretical peak of regular FP64 CUDA cores.

What carries the argument

The quartering decomposition (splitting a 52-bit double fraction into four 13-bit quarters with 26 trailing zero bits each), the block-convolution rewrite mapping a product of double-double matrices into a single product of double matrices (equations 2-4), and the balancing algorithms of Section 3 that redistribute bits to keep quarter exponents on the grid {0, -13, -26, -39}.

What would settle it

An input matrix containing doubles whose quarter exponents fall into worst-case gaps (e.g., -25, -38, -51 as described in Problem 3.4) and for which the balancing corrections cascade without converging, causing the trailing-zero guarantee to break and introducing rounding errors into the tensor-core computation.

Watch

Extended reading notes

Core claim

The central mechanism is the quartering decomposition: splitting each double's 52-bit fraction into four 13-bit pieces with 26 trailing zero bits each. This padding makes products of quarters exact in double precision and allows the renormalization step (which normally requires branching) to be deferred entirely, because no information is lost during the tensor-core matrix multiplication itself. The renormalization is then performed only once, after the matrix product is complete, as a massively parallel summation on regular cores. The product of two multiple-double matrices is rewritten as a single large double-precision matrix product using a block structure where one factor stacks quarter

Load-bearing premise

The balancing algorithms are assumed to reliably force all quarter exponents onto the grid {0, -13, -26, -39} for arbitrary inputs, ensuring exact inner products. The paper acknowledges worst-case exponent gaps and proposes corrective bit-shuffling steps, but does not formally prove these corrections always terminate or guarantee exactness, and tests only random inputs rather than adversarial ones.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 6 minor

Summary. The paper proposes a method for performing multiple double arithmetic on NVIDIA FP64 tensor cores, which are specialized for matrix multiplication and unsuited to the branching required by renormalization. The core idea is to partition each double into 13-bit quarters with trailing zero bits, rewrite the product of two multiple-double matrices as a single product of enlarged double matrices, and execute it on tensor cores. The approach is inspired by the Ozaki scheme. The paper presents both a bottom-up method (custom WMMA kernel) and a top-down method (adapted dmmaTensorCoreGemm), reporting 13.75 TFLOPS on an A100 for the latter.

Significance. The central technical idea—using 13-bit quartering with 26 trailing zero bits to enable exact inner products on tensor cores (Section 2, Eq. 5)—is clean and well-motivated. The overflow bound N = 67,108,865 is a direct and useful calculation. The software is released under GPL on GitHub, which supports reproducibility. The problem is timely given the growing interest in mixed-precision and extended-precision arithmetic on accelerators.

major comments (3)
  1. §4.3: The headline performance claim of 13.75 TFLOPS is reported for the top-down method using dmmaTensorCoreGemm with shared memory, but this figure appears to measure only the raw matrix multiply kernel. The paper does not report full-pipeline timing for the top-down method that includes the costs of splitting doubles into quarters, constructing the enlarged A and B matrices (8192×4096 per Table 2), applying balancing algorithms, and reassembling results. These overhead costs are precisely what caused the bottom-up method to be 3–4x slower than regular CUDA cores (Table 1: t_TC vs t_CUDA). Since the top-down method uses the same quartering and matrix expansion, there is no reason to expect its full-pipeline timing to differ qualitatively. The comparison '13.75 > 9.7' is therefore not apples-to-apples: the 9.7 TFLOPS figure for regular CUDA cores includes all overhead, while the 13.75 T
  2. §4.3: No accuracy validation (ε_max) is reported for the top-down method. The bottom-up method reports ε_max in Table 1 (e.g., 3.5E-28 for n=1024 double-double), but the top-down method, which produces the headline performance number, has no corresponding accuracy data. Without this, the reader cannot confirm that the top-down method actually achieves double-double accuracy.
  3. §3, Problems 3.1–3.4: The balancing algorithms are described informally through problem-solution pairs but lack formal proofs of correctness or termination. Problem 3.4 acknowledges worst-case exponent gaps (e.g., -25, -38, -51) and proposes corrective steps (Eq. 8), but it is not shown that these corrections always converge to grid-aligned exponents for arbitrary inputs. Additionally, Section 4.1 states that only random numbers were used in experiments. Adversarial inputs (e.g., numbers with maximal exponent gaps, cancellation scenarios) would stress the balancing algorithms and provide evidence for or against their robustness.
minor comments (6)
  1. §2: The notation in Eqs. (1)–(3) uses subscripts like a_{i,k,0} but the matrix A_{i,k} in Eq. (2) is written with a transpose that may confuse readers; clarifying whether the transpose applies to the column vector or the indexing would help.
  2. Table 1: The units for ε_max are not specified. Are these absolute errors or relative errors? Adding units or a reference to the double-double precision range would contextualize the values.
  3. §4.2: The bottom-up method reports f_WMMA and f_CUDA in teraflops, but the basis for the FLOP count (e.g., 2n³ for n×n matrix multiply) is not stated. Clarifying the formula used would make the performance figures verifiable.
  4. §3, Problem 3.3: The corrective step in Eq. (7) assumes x_h > 0. The case x_h < 0 is not discussed; presumably the sign is handled symmetrically, but stating this explicitly would be clearer.
  5. The abstract states the software is available 'under the GPU GPL license'; this should likely be 'GPL license' or specify the version (e.g., GPL v3).
  6. §1: The phrase 'multiple the accuracy' should be 'multiply the accuracy'.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity found: the derivation is self-contained, parameter-free, and validated against an independent implementation.

full rationale

The paper's derivation chain is self-contained and does not exhibit circularity. The 13-bit quartering scheme follows directly from the IEEE 754 double format (52-bit fraction split into 4×13 bits), and the overflow bound (Eq. 5) is a direct arithmetic calculation yielding N ≤ 67,108,865. The Ozaki scheme is cited as inspiration from external work (Ozaki et al., Numerical Algorithms, 2012 — authors are Chen and Verschelde, not Ozaki). No fitted parameters are introduced anywhere in the derivation. The accuracy claim is validated against an independent CUDA core implementation (Eq. 9 defines ε_max as the max element-wise difference between the tensor core result C_TC and the CUDA core result C_CUDA). The balancing algorithms in Section 3 are constructive bit-manipulation procedures, not fitted models. While the skeptic correctly identifies that the 13.75 TFLOPS headline from the top-down method excludes pipeline overhead that the bottom-up method reveals as dominant (Table 1: t_TC is 3-4× slower than t_CUDA), this is a correctness/completeness concern about performance reporting, not a circularity issue — the performance claim is not defined in terms of itself, and the accuracy validation uses an independently computed reference. The self-citations (refs [16]-[20] by Verschelde) are to software packages (PHCpack) and prior applications, not to load-bearing theorems or uniqueness results that would force the present conclusion. No step in the derivation reduces to its inputs by construction.

Assumptions & free parameters 2 free parameters · 4 assumptions · 0 invented entities

The paper introduces no new physical entities or mathematical objects beyond the 'quarters' terminology, which is a renaming of standard bit-slicing. The axioms are mostly standard or domain assumptions, with one ad-hoc assumption about the balancing algorithms that lacks formal proof.

free parameters (2)
  • Quarter bit width (13 bits) = 13
    Chosen so that 4 quarters fill the 52-bit fraction of a double, and products of two quarters have at least 26 trailing zero bits. This is determined by the IEEE 754 format, not fitted to data.
  • Overflow bound N = 67108865
    Derived from equation (5) as the maximum number of terms that can be summed exactly. Not a free parameter but a consequence of the 13-bit choice.
assumptions (4)
  • standard math IEEE 754 double-precision floating-point format has a 52-bit fraction and 11-bit exponent
    Used throughout Section 2 to justify the 13-bit quartering and trailing zero bits.
  • domain assumption FP64 tensor cores on A100 perform IEEE-compliant double-precision matrix multiply-accumulate
    Stated in Section 1 citing NVIDIA A100 whitepaper [12]. The exactness of the computation depends on the tensor core performing true IEEE-compliant additions without intermediate rounding.
  • domain assumption The Ozaki scheme can be adapted from mixed-precision to multiple-double arithmetic
    Invoked in Section 1 and Section 2 as the conceptual basis for the rewriting approach.
  • ad hoc to paper Balancing corrections (Problems 3.1-3.4) preserve the value of the multiple double and can be applied iteratively to achieve grid-aligned exponents
    Section 3 describes corrective steps for sign mismatches and exponent gaps but does not prove they always converge or preserve exactness.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Multiple Double Arithmetic on NVIDIA Tensor Cores." pith.science (2026). https://pith.science/paper/Y4ST3CVY

@misc{pith2026260706881,
  author       = {Pith},
  title        = {Pith review of: Multiple Double Arithmetic on NVIDIA Tensor Cores},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/Y4ST3CVY}},
  note         = {Machine review of arXiv:2607.06881}
}
read the original abstract

A multiple double is an unevaluated sum of doubles. An NVIDIA tensor core is a specialized high performance compute core for matrix multiplication. The Ampere A100, released in 2020, introduced tensor cores capable of 64-bit floating-point arithmetic. Every multiple double arithmetical operation requires renormalization, which involves branching, for which tensor cores are unsuited. To solve this problem caused by renormalization, we apply a solution similar to the Ozaki scheme [Ozaki et al, Numerical Algorithms, 2012]. Our software is available under the GPU GPL license on github.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

24 extracted references · 24 canonical work pages

  1. [1]

    Abdelfattah, H

    A. Abdelfattah, H. Anzt, E. G. Boman, E. Carson, T. Cojean, J. Dongarra, A. Fox, M. Gates, N. J. Higham, X. S. Li, J. Loe, P. Luszczek, S. Pranesh, S. Rajamanickam, T. Ribizel, B. F. Smith, K. Swirydowicz, S. Thomas, S. Tomov, Y. M. Tsai, and U. M. Yang. A survey of numerical linear algebra methods utilizing mixed-precision arithmetic. International Journ...

    work page 2021

  2. [2]

    T. J. Dekker. A floating-point technique for extending the available precision.Numerische Mathematik, 18(3):224–242, 1971

    work page 1971

  3. [3]

    Y. Hida, X. S. Li, and D. H. Bailey. Algorithms for quad-double precision floating point arithmetic. In15th IEEE Symposium on Computer Arithmetic (Arith-15 2001), pages 155–

    work page 2001

  4. [4]

    IEEE Computer Society, 2001

    work page 2001

  5. [5]

    N. J. Higham and T. Mary. Mixed precision algorithms in numerical linear algebra.Acta Numerica, pages 347–414, 2022

    work page 2022

  6. [6]

    Joldes, J.-M

    M. Joldes, J.-M. Muller, V. Popescu, and Tucker. W. CAMPARY: Cuda Multiple preci- sion arithmetic library and applications. InMathematical Software – ICMS 2016, the 5th International Conference on Mathematical Software, pages 232–240. Springer-Verlag, 2016

    work page 2016

  7. [7]

    Kashi, H

    A. Kashi, H. Lu, W. Brewer, D. Rogers, M. Matheson, M. Shankar, and F. Wang. Mixed- precision numerics in scientific applications: survey and perspectives.The Journal of Su- percomputing, 82(287):1–67, 2026

    work page 2026

  8. [8]

    T. Kouya. Acceleration of multicomponent multiple-precision arithmetic with branch-free algorithms and SIMD vectorization.arXiv:2603.14926v2 [cs.MS] 7 May 2026

    work page Pith review arXiv 2026

Show all 24 references
  1. [9]

    N. Maho. MPLAPACK version 2.0.1. user manual.arXiv:2109.13406v2 [cs.MS] 12 Sep 2022. 8

    2022 arXiv

  2. [10]

    Mukunoki, K

    D. Mukunoki, K. Ozaki, T. Ogita, and T. Imamura. DGEMM using tensor cores, and its accurate and reproducible versions. In P. Adayappan, B. L. Chamberlain, G. Juckeland, and H. Ltaief, editors,High Performance Computing. ISC High Performance 2020, volume 12151 ofLecture Notes i...

    2020

  3. [11]

    Muller, N

    J.-M. Muller, N. Brunie, F. de Dinechin, C.-P. Jeannerod, M. Joldes, V. Lefèvre, G. Melquiond, N. Revol, and S. Torres.Handbook of Floating-Point Arithmetic. Springer- Verlag, second edition, 2018

    2018

  4. [12]

    CUDA Samples

    NVIDIA. CUDA Samples. Athttps://github.com/NVIDIA/cuda-samples

  5. [13]

    NVIDIA A100 Tensor Core GPU Architecture, 2020

    NVIDIA. NVIDIA A100 Tensor Core GPU Architecture, 2020. Whitepaper available via https://www.nvidia.com

    2020

  6. [14]

    Ozaki, T

    K. Ozaki, T. Ogita, S. Oishi, and S. M. Rump. Error-free transformations of matrix mul- tiplication by using fast routines of matrix multiplication and its applications.Numerical Algorithms, 59:95–118, 2012

    2012

  7. [15]

    S. M. Rump. Verification methods: Rigorous results using floating-point arithmetic.Acta Numerica, 19:287–449, 2010

    2010

  8. [16]

    J. R. Shewchuk. Adaptive precision floating-point arithmetic and fast robust geometric predicates.Discrete Comput. Geom., 18(3):305–363, 1997

    1997

  9. [17]

    Verschelde

    J. Verschelde. Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation.ACM Trans. Math. Softw., 25(2):251–276, 1999. Available at https://github.com/janverschelde/PHCpack

    1999

  10. [18]

    Verschelde

    J. Verschelde. Accelerated polynomial evaluation and differentiation at power series in mul- tiple double precision. InThe 2021 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 740–749. IEEE, 2021

    2021

  11. [19]

    Verschelde

    J. Verschelde. Least squares on GPUs in multiple double precision. InThe 2022 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 828–837. IEEE, 2022

    2022

  12. [20]

    Verschelde

    J. Verschelde. GPU accelerated Newton for Taylor series solutions of polynomial homotopies in multiple double precision. In F. Boulier, C. Mou, T. M. Sadykov, and E. V. Vorozhtsov, editors,Proceedings of the 26th International Workshop on Computer Algebra in Scientific Computi...

    2024

  13. [21]

    Springer-Verlag, 2024

    2024

  14. [22]

    Verschelde

    J. Verschelde. Multiword arithmetic and parallel computing.ACM SIGAda Ada Letters, 45(2):67–68, 2025

    2025

  15. [23]

    D. K. Zhang. MultiFloats.jl.https://github.com/dzhang314/MultiFloats.jl

  16. [24]

    D. K. Zhang and A. Aiken. High-performance branch-free algorithms for extended-precision floating-point arithmetic. InProceedings of the International Conference for High Perfor- mance Computing, Networking, Storage and Analysis, pages 695–710. ACM, 2025. 9

    2025

Pith tools

Reviewed July 9, 2026 · model on record in the stance chip above.