mpoz2 — Multiple Precision Ozaki scheme II

mpoz2 is a C library that computes multiple-precision GEMM and LU decomposition by mapping the arithmetic onto the low-precision matrix engines of modern hardware (Intel AMX-INT8, Arm SVE2 i8mm, NVIDIA tensor cores), using Ozaki scheme II: exact integer modular products followed by an explicit CRT reconstruction, with a single rounding at the end.

Because the low-precision products are exact and only the final result is rounded, the computed matrix product is more accurate than a naive multiple-precision GEMM, which rounds at each of the K multiply-adds — while also being faster.

Features

• Arbitrary precision through MPFR (mpfr_t in, mpfr_t out).
• Multi-component precision (DD/TD/QD) with a direct conversion between the non-overlapping expansion format and the internal fixed-point representation — no MPFR round trip, no per-element heap allocation.
• Blocked LU decomposition with partial pivoting (LAPACK getrf style), where only the Schur complement update goes through Ozaki scheme II.
• Selectable back-ends: INT8 (default), binary64, and, on the GPU, FP16, FP8 (E4M3) and FP4 (E2M1). The INT8 back-end switches to binary64 automatically when the target precision exceeds its modulus-table capacity.
• CPU kernels for x86-64 (AMX-INT8, TDPBSSD) and aarch64 (SVE2 i8mm, SMMLA), selected at build time from uname -m.
• GPU support as header-only CUDA (cuda/*.cuh), operating directly on gdd_real/gtd_real/gqd_real of gdtq and on cu_freal<PB> of MPC_CUDA.

Requirements

mpoz2 is split into a CPU library, which is self-contained, and an optional CUDA part under cuda/. Nothing in the CPU library depends on the GPU packages, so the CPU side can be built and used on a machine with no GPU and no CUDA toolkit at all.

CPU library (required)

Package	Purpose	Make variable	Default
GMP	multiple-precision integers behind MPFR and the CRT	GMP_PREFIX	/usr/local
MPFR	arbitrary-precision input/output type (mpfr_t)	GMP_PREFIX	/usr/local
OpenBLAS	cblas_dgemm for the binary64 back-end	BLAS_PREFIX	/opt/OpenBLAS (aarch64), /usr/local (x86-64)

A C11 compiler with OpenMP is also required. For the INT8 back-end the CPU needs AMX-INT8 (x86-64, Sapphire Rapids or later) or SVE2 with the i8mm extension (aarch64); the kernel is chosen at build time from uname -m. If the matrix engine turns out to be unavailable at run time, the library reports this once and falls back to the binary64 back-end.

make GMP_PREFIX=$HOME/local BLAS_PREFIX=/opt/OpenBLAS

dtq interoperability (optional)

Package	Purpose	Make variable	Default
dtq	dd_real / td_real / qd_real	DTQ_PREFIX	/usr/local

mpoz2 never includes a dtq header. dtq stores the components of its types in a contiguous double x[nc], so an array of qd_real is an array of double[4] and can be passed to oz2_gemm_qd() by a reinterpret cast, with no conversion and no copy. examples/example_dtq.cpp demonstrates and checks this:

make dtq-example DTQ_PREFIX=$HOME/local
./examples/example_dtq 256

CUDA part (optional, cuda/)

Package	Purpose	Make variable	Default
CUDA 12+	nvcc, cuBLASLt	CUDA_ARCH	sm_90
MPC_CUDA	cu_freal<PB>, the GPU fixed-precision type	MPC_CUDA_INC	/usr/local/include
gdtq	gdd_real / gtd_real / gqd_real	GDTQ_INC	/usr/local/include/gdtq
MPFR / GMP	host-side conversion and reference	GMP_PREFIX	/usr/local

make -C cuda CUDA_ARCH=sm_90 MPC_CUDA_INC=$HOME/local/include \
             GDTQ_INC=$HOME/local/include/gdtq
sudo make install-cuda            # installs cuda/include/*.cuh

cuda/ is header-only; make -C cuda builds the example drivers only.

Build and install

make                               # builds libmpoz2.a (CPU only)
make check                         # runs the CPU test suite
sudo make install PREFIX=/usr/local

The CUDA part is built separately and is never required:

make -C cuda                       # optional

Quick start

#include <mpoz2.h>

mpfr_t *A, *B, *C;              /* M*K, K*N, M*N arrays, precision p */
oz2_opts opts;
oz2_opts_default(&opts);        /* INT8 back-end, 16 guard bits,
                                   automatic slice count */

oz2_gemm(M, N, K, A, K, B, N, C, N, &opts);

Slice splitting

The significand is cut into opts.nslice chunks of w bits and each digit group is reconstructed separately. The conversion cost (residue computation and CRT) falls as 1/nslice because fewer moduli are needed per pass, while the number of low-precision GEMMs grows with nslice, so there is an interior optimum. opts.nslice = 0 (the default) picks it from a measured cost model; 1 reproduces the original single-pass construction; a larger value forces that many slices. OZ2_NSLICE overrides the heuristic at run time.

Slice splitting also removes the capacity ceiling of the INT8 back-end (previously t <~ 11607 bits), because the chunk width is set by the slice count rather than by the target precision.

Two notes on accuracy. With nslice = 1 the product is reconstructed exactly by the CRT, so the result is correct regardless of cancellation in the inner products. With nslice > 1 the low-order digit groups are truncated, and the truncation is bounded relative to K * 2^(2t) rather than to the computed value; problems with severe cancellation therefore need opts.guard_bits (or the retained group count, OZ2_EXTRA_GROUPS) raised accordingly. One extra guard bit buys one bit of accuracy, one extra retained group buys about w bits.

For DD/TD/QD, pass the expansion arrays directly — the components of dd_real/td_real/qd_real are contiguous doubles, so an array of them is an array of double[nc]:

oz2_gemm_qd(M, N, K, (const double (*)[4])A, K,
                     (const double (*)[4])B, N,
                     (double (*)[4])C, N, &opts);

Blocked LU with partial pivoting:

oz2_lu(n, A, lda, ipiv, /*panel width*/ 256, OZ2_LU_SCHUR_OZAKI, &opts);
oz2_lu_solve(n, A, lda, ipiv, x);

See examples/example.c for a complete program.

Layout

Path	Contents	Needs
include/	Public headers (oz2.h, mpoz2.h)	---
src/	CPU implementation (context, split, CRT, GEMM, LU, I/O, kernels)	GMP, MPFR, CBLAS
test/	Correctness tests (make check)	as above
examples/	example.c (MPFR), example_dtq.cpp (dtq interop)	dtq for the latter
cuda/include/	Header-only CUDA implementation	CUDA, MPC_CUDA, gdtq
cuda/examples/	GPU LU drivers	as above

Choosing a back-end

The right back-end is the one that maximises bits per modulus × engine throughput, and it changes with the machine and the target precision:

• DD/TD/QD: INT8 is normally fastest; on GPUs with strong FP64 the binary64 back-end can win at QD.
• Arbitrary precision on the CPU: binary64 is faster than INT8 at every precision we measured (p ≥ 1024), because INT8's two-digit moduli need three GEMM planes each.
• Arbitrary precision on the GPU: INT8.
• FP4: not recommended; its base-13 two-digit decomposition costs four sub-GEMMs per modulus.

oz2_auto_select() provides a first-cut choice from N and the target precision.

Citing

If you use mpoz2 in your research, please cite:

@misc{mpoz2,
  author = {Tomonori Kouya},
  title  = {mpoz2: Multiple Precision Ozaki scheme II},
  version= {0.0.1},
  year   = {2026},
  url    = {https://github.com/tkouya/mpoz2}
}

and the algorithm it implements:

K. Ozaki, Y. Uchino and T. Imamura, "Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique", arXiv:2504.08009, 2025.

License

Modified BSD (3-clause). See LICENSE. Note that GMP and MPFR are LGPL v3, so redistributing binaries requires complying with their terms too.

Acknowledgements

Supported by JSPS KAKENHI Grant Number JP26K14846.