torch-dd

Double-double (emulated quad precision) arithmetic for PyTorch.

torch-dd represents each number as an unevaluated sum of two float64 values (hi + lo), giving roughly 106 significant bits — about twice the precision of float64 — while running on ordinary CPU and CUDA tensors. It is built on error-free transformations (two_sum, two_prod, Dekker split) and the QD-library algorithms, so results carry a proven relative-error bound instead of drifting.

Use it like a built-in numeric type through the DDTensor tensor subclass: operators, comparisons, torch functions, autograd and most of torch.nn.functional route through __torch_dispatch__ and run at DD precision. Build one with tensor(...) or the torch-style factories (zeros, ones, full, eye, arange, linspace, rand, randn, …); its .hi and .lo are the two float64 words, and .item() reads the value back as a rounded Python float.

import torch, torch_dd as dd

# float64 cannot even represent 1 + 2**-60, so the difference collapses to 0.
one = torch.ones(1, dtype=torch.float64)
tiny = torch.full((1,), 2.0**-60, dtype=torch.float64)
print(((one + tiny) - one).item())                       # 0.0   (wrong)

one_dd = dd.ones(1)
tiny_dd = dd.full((1,), 2.0**-60)
r = (one_dd + tiny_dd) - one_dd
print(r.item())                                          # 8.67e-19  (= 2**-60)

Why

float64 silently loses low-order bits in cancellation, long accumulations, and ill-conditioned formulas. The classic example is the small root of x² + 10⁸·x + 1 = 0 via the textbook formula: the subtraction amplifies the rounding error of sqrt by ~10¹⁶, so float64 keeps no correct digits while DD still keeps ~15.

import torch, torch_dd as dd

x = dd.full((1,), 1e8)
root = (-x + torch.sqrt(x * x - 4.0)) / 2.0                  # like builtin math
# relative error vs. the exact root:  float64 ~0.25   |   torch-dd ~4e-18

Install

uv sync           # or: pip install -e .

Requirements: Python ≥ 3.12, PyTorch ≥ 2.13. The core runs in pure PyTorch. On first use torch-dd JIT-compiles fused C++/CUDA kernels via torch.utils.cpp_extension for speed; if no compiler/CUDA toolkit is present it transparently falls back to the pure-PyTorch implementation, which produces bit-identical results.

Features

• Operators & reductions: + - * / @, comparisons, sum, and matmul on a DDTensor, all with proven O(2**-106) error bounds.
• Transcendentals: torch.exp, torch.log, torch.sin/cos/tan, torch.sinh/cosh/tanh, torch.rsqrt, torch.pow, torch.erf/erfc on a DDTensor evaluate in DD from mpmath-verified constants.
• Autograd: reverse-mode differentiation runs at DD precision; gradients come back as DDTensors.
• torch.nn.functional: forward and backward for the whole functional surface — activations, softmax, losses, normalizations, pooling, convolution, attention, interpolation, grid_sample, embeddings, ctc_loss, and more.
• Fast matmul: @ uses an Ozaki-scheme product that turns one DD matmul into error-free lower-precision GEMMs — INT8 tensor cores on CUDA, blocked DGEMM on CPU — above a size threshold, with a dynamic-range guard and FP64 fallback.

import torch, torch_dd as dd
from torch_dd import DDTensor
import torch.nn.functional as F

w = dd.randn(4, 3, requires_grad=True)
x = dd.randn(2, 3)
F.mse_loss(F.linear(x, w), dd.zeros(2, 4)).backward()
assert isinstance(w.grad, DDTensor)      # a double-double gradient

Examples

examples/ trains an MNIST classifier entirely in double-double precision — weights, activations, gradients and optimizer state are all DDTensors — using a from-scratch DD AdamW optimizer:

python examples/train_mnist.py            # quick subset, a few seconds on GPU
python examples/train_mnist.py --full     # all 60k images (~97% test accuracy)

See examples/README.md for details.

Accuracy

Every operation stays within the DD contract: element-wise ops are correct to a few units of 2**-106 relative, and reductions/mm to O(k · 2**-105 · Σ|terms|). IEEE inf/nan propagate as in float64. The representable range matches float64 (the lo word degrades below ~2**-970, the documented double-double limit).

Development

uv run pytest             # the full test suite (CPU + CUDA)
uv run black src tests examples
uv run mypy               # file set and config live in pyproject.toml

References verify DD results against 256-bit mpmath (and exact Fraction arithmetic for the error-free transformations), well beyond float64 resolution.

Layout

src/torch_dd/
  eft.py        error-free transformations (two_sum, two_prod, split)
  core.py       DD type and arithmetic (add/mul/div/sqrt/sum/mm, ...)
  mathfn.py     transcendentals (exp/log/trig/hyperbolic/erf/pow)
  dd_const.py   exact hex-pair constants (mpmath-generated)
  tensor.py     DDTensor subclass + __torch_dispatch__ handlers
  factory.py    tensor-creation factories (zeros/ones/arange/linspace/...)
  ozaki.py      Ozaki-scheme mm engines (CUDA INT8 / CPU DGEMM)
  native.py     JIT loader for the fused C++/CUDA kernels
  csrc/         C++/CUDA kernel sources
examples/       DD AdamW + MNIST training
tests/          pytest suite (mpmath / Fraction references)