arXiv CSopen access

Improved Scaling for Fast Mode of Ozaki Scheme II

Journal Title

Improved Scaling for Fast Mode of Ozaki Scheme II

XX(X):1–13 ©The Author(s) 2016 Reprints and permission: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/ToBeAssigned www.sagepub.com/

SAGE 1

2

Shota Kawakami and Daisuke Takahashi

Abstract Ozaki scheme II emulates high-precision matrix multiplication using low-precision integer matrix operations based on the Chinese remainder theorem (CRT). It first scales the high-precision matrices to convert them into integer matrices.

arXiv:2606.29129v1 [cs.MS] 28 Jun 2026

For this scaling step, Ozaki scheme II provides two modes: accurate mode, which uses INT8 matrix multiplication to estimate scaling factors, and fast mode, which applies the Cauchy–Schwarz inequality at lower computational cost. We show that the existing formula lacks scale invariance; multiplying the input matrices by a constant changes the effective bit width of the integer matrices in the scaling step, causing accuracy degradation or CRT recovery failure. To address this, we propose a revised scaling formula derived from the CRT uniqueness condition via the Cauchy– Schwarz inequality. The proposed formula is scale-invariant by construction, guarantees that the CRT uniqueness condition is always satisfied, and introduces no additional overhead over the original fast mode. Experiments on an NVIDIA GH200 GPU show that the proposed method achieves accuracy comparable to that of accurate mode while maintaining throughput comparable to that of fast mode. In the accuracy–throughput trade-off, the proposed method overcomes the accuracy limitation of fast mode and the throughput constraint of accurate mode, offering a superior accuracy and performance. Keywords matrix multiplication, floating-point arithmetic, matrix engine, Ozaki scheme

1

Introduction

Modern accelerators such as graphics processing units (GPUs) have been designed with increasingly powerful low-precision arithmetic units, driven by growing demand from the artificial intelligence community. For example, on the NVIDIA GH200 GPU, the FP64 Tensor Core achieves 67 TFLOPS (tera floating-point operations per second) and the INT8/FP8 Tensor Cores achieve 1979 TOPS (tera operations per second)/TFLOPS (NVIDIA 2025). On the more recent NVIDIA B200 GPU, the performance gap is wider: the FP64 Tensor Core achieves 37 TFLOPS and the INT8/FP8 Tensor Cores achieve 4500 TOPS/TFLOPS (NVIDIA 2024).

The Ozaki scheme emulates high-precision matrix multiplication using low-precision matrix operations. When there is a large performance gap between low- and highprecision units, exploiting the high-throughput low-precision hardware enables fast and energy-efficient high-precision matrix multiplication. Mukunoki et al. (2020) were the first to implement single- and double-precision general matrix multiply (SGEMM and DGEMM, respectively) using the Ozaki scheme with FP16 Tensor Cores, which are lowprecision matrix computation units. Ootomo et al. (2024) implemented DGEMM using even lower-precision INT8 Tensor Cores. More recently, Mukunoki (2026) proposed an implementation of DGEMM using FP8 Tensor Cores.

Scientific computing requires high-precision arithmetic such as FP64 or FP32. Methods that emulate high-precision computation using low-precision hardware, as well as mixedprecision algorithms that exploit low-precision operations for part of the computation, have thus attracted considerable

1 Graduate School of Science and Technology, University of Tsukuba, Japan 2 Center for Computational Sciences, University of Tsukuba, Japan

scheme (Ozaki et al. 2012, 2025) has emerged as a promising

Corresponding author: Shota Kawakami, Graduate School of Science and Technology, University of Tsukuba Tennodai 1-1-1, Tsukuba, Ibaraki, 305-8573, Japan

approach.

Email: [email protected]

attention (Dongarra et al. 2024). Among them, the Ozaki

Prepared using sagej.cls [Version: 2017/01/17 v1.20]

2

Journal Title XX(X)

These implementations achieve performance that exceeds

• We show that the existing scaling formula for fast

that of cuBLAS on NVIDIA GPUs. Since the core of the

mode lacks scale invariance, which can cause accuracy

implementations relies on BLAS (Basic Linear Algebra

degradation or even CRT recovery failure.

Subprograms), no highly specialized implementation or • We propose an improved scaling formula for fast mode

hand-tuned optimization is required, making the approach

that addresses this limitation.

highly portable. The Ozaki scheme has also been applied to matrix multiplication in scientific computing applications,

• We compare the accuracy, throughput, and their trade-

including quantum circuit simulation (Ootomo et al. 2024),

off of the proposed method against accurate mode,

quantum chemistry calculations (Dawson et al. 2024), and ab

the original fast mode, and cuBLAS on an NVIDIA

initio tensor network state methods (Brower et al. 2026).

GH200 GPU.

Recently, Ozaki scheme II (Ozaki et al. 2025) was Throughout this paper, we focus on DGEMM and

proposed as an improvement over the conventional Ozaki scheme (Ozaki et al. 2012). Ozaki scheme II converts highprecision matrices into integer matrices and decomposes the computation into a series of low-precision integer matrix multiplications via the Chinese remainder theorem (CRT). Compared to the conventional Ozaki scheme, Ozaki

SGEMM using INT8 matrix engines. However, since an analogous scaling formula is used in the fast mode of DGEMM and SGEMM with FP8 matrix engines and of ZGEMM and CGEMM with INT8 matrix engines, the proposed improvement is applicable to these cases as well. The rest of this paper is organized as follows. Section 2

scheme II improves performance by reducing the number of low-precision matrix multiplications required. On the NVIDIA GH200 GPU, cuBLAS DGEMM achieves a peak performance of 60.9 TFLOPS, whereas DGEMM based on Ozaki scheme II reaches up to 80.2 TFLOPS (Ozaki et al. 2025). Uchino et al. proposed methods for implementing DGEMM and SGEMM via Ozaki scheme II using INT8 matrix engines (Uchino et al. 2025) and FP8 matrix engines (Uchino et al. 2026b), and extended the approach to single- and double-precision complex general matrix

introduces notations and provides an overview of Ozaki scheme II and its two scaling modes. Section 3 analyzes the fast mode scaling formula and identifies its limitation through theoretical analysis and numerical demonstration. Section 4 proposes an improved scaling formula that addresses this limitation. Section 5 evaluates the proposed method in terms of accuracy, throughput, and their trade-off, comparing it against accurate mode, the original fast mode, and cuBLAS. Section 6 concludes the paper.

multiply (CGEMM and ZGEMM, respectively) (Uchino et al. 2026a). These methods are available in the open-source

2

Preliminaries

2.1

Notations

library GEMMul8 (Uchino 2026). In Ozaki scheme II, the high-precision matrices are first scaled to convert them into integer matrices. During this scaling process, the scaling factors must be chosen to satisfy the conditions for correct CRT recovery. Ozaki scheme II provides two modes for computing the scaling factors: accurate mode, which uses INT8 matrix multiplication to estimate scaling factors, and fast mode, which applies the Cauchy–Schwarz inequality at lower computational cost. The trade-off between accuracy and computation time is different between these two modes. A performance evaluation of the two modes is provided by Ozaki et al. (2025) and the implementation details are described by Uchino et al. (2025).

We follow the notations of Uchino et al. (2025). F denotes a floating-point number system, such as F32 for FP32 and F64 for FP64. The operations mod and rmod (·, ·) are respectively defined as x mod p := x − p · ⌊x/p⌋ , rmod (x, p) := x − p · round (x/p) , where round(·) denotes rounding to the nearest integer with ties broken by round-to-nearest-even. trunc(·) denotes truncation toward zero. diag(·) constructs a diagonal matrix from a vector. single(·) evaluates its argument in single precision with round-to-nearest-even rounding, while

In this paper, we propose an improved fast mode scaling formula for DGEMM and SGEMM via Ozaki scheme II with

single△ (·) and single▽ (·) denote single-precision evaluation

INT8 matrix engines that achieves better accuracy without

of these functions is applied to a matrix, it acts element-wise.

sacrificing throughput. The main contributions of this paper

ai,: and b:,j denote the i-th row of A and the j-th column of

are as follows.

B, respectively. Throughout this paper, we assume that an

Prepared using sagej.cls

with round-up and rounding-down, respectively. When any

Kawakami and Takahashi

3

INT8 matrix engine accepts two INT8 matrices as inputs and

Step 4. Rescale C ′′ using µ and ν to obtain C ≈ AB as follows:

accumulates the inner products in INT32.

C := diag (µ)

2.2

−1

· C ′′ · diag (ν)

−1

.

Ozaki Scheme II

Ozaki scheme II (Ozaki et al. 2025) is a high-precision

By choosing pi as positive integers representable in INT8,

matrix multiplication emulation algorithm based on the

the matrix products rmod (A′ , pi ) · rmod (B ′ , pi ) in (2) can

following CRT, using low-precision matrix multiplication.

be computed using an INT8 matrix engine, provided that

Theorem 1. Chinese remainder theorem. Let p1 , . . . , pN ∈ QN N≥2 be pairwise coprime integers, and let P := i=1 pi . For x ∈ Z, we define

k ≤ 217 so that the INT32 accumulator does not overflow; for k > 217 , block matrix multiplication can be applied to avoid overflow. By choosing µ and ν to satisfy (1), A′ B ′ can be exactly

xi := x mod pi

recovered in (3). The bit width of A′ and B ′ affects the

(i = 1, . . . , N ).

accuracy of C; increasing it requires increasing the number Let Pi := P/pi and let qi be the modular inverse of Pi modulo pi . Then, the following holds: x≡

N X

of moduli N , and hence enlarging P . The dominant cost is Step 3, which performs INT8 general matrix multiply (GEMM) on A′ and B ′ for each modulus

Pi qi xi mod P.

in (2). Since the O(mnk) INT8 GEMM is performed N

i=1

times (once per modulus), the total computational cost is

Ozaki scheme II computes the product C ≈ AB ∈ Fm×n

O(N mnk).

of matrices A = (aij ) ∈ Fm×k and B = (bij ) ∈ Fk×n via the CRT in the following steps.

2.3 Step 1. Choose pairwise coprime integers p1 , . . . , pN ∈ QN N≥2 . Compute P := i=1 pi , Pi := P/pi , and qi := Pi−1 mod pi .

Scaling in Ozaki Scheme II

Ozaki scheme II provides two modes for determining µ and ν: accurate mode and fast mode. The two modes differ in how they evaluate condition (1).

Step 2. Using µ = (µi ) ∈ Zm and ν = (νj ) ∈ Zn , convert

In accurate mode, an auxiliary INT8 GEMM is used to

and B ′ = (b′ij ) ∈ Zk×n , respectively, as follows:

evaluate condition (1). First, we define µ′ = (µ′i ) ∈ Fm and
ν ′ = νj′ ∈ Fn from A and B, respectively, as follows:

a′ij := trunc (µi aij ) ,

µ′i := 25−⌊log2 (max1≦h≦k |aih |)⌋ ,

b′ij := trunc (νj bij ) .

νj′ := 25−⌊log2 (max1≦h≦k |bhj |)⌋ .

Here, µ and ν are chosen such that A′ and B ′ satisfy

Here, ⌊log2 (·)⌋ can be computed directly from the exponent

the following condition:

field of a floating-point number. Then, we convert A and
B into INT8 matrices Ā = (āij ) ∈ Zm×k and B̄ = b̄ij ∈

A and B into integer matrices A

k X

a′ih b′hj ≤

h=1

k X

′

= (a′ij ) ∈ Zm×k

|a′ih | b′hj <

h=1

P 2

∀i, j.

(1)

Step 3. Compute C ′′ ≡ A′ B ′ mod P using the CRT as

Zk×n using µ′ and ν ′ , respectively, as follows: āij := ⌈µ′i |aij |⌉ ≤ 27 − 1,   b̄ij := νj′ |bij | ≤ 27 − 1.

follows: C ′ :=

N X

Using C̄ := ĀB̄ = (c̄ij ) ∈ Zm×n , which can be computed Pi qi · rmod (A′ , pi ) · rmod (B ′ , pi ) ,

i=1

′

(2) C ′′ := rmod (C ′ , P ) .

Prepared using sagej.cls

via an INT8 matrix engine, µ and ν are defined as follows:

(3)

µi := µ′−1 · 2⌊Paccu −0.5·log2 (max1≦h≦n c̄ih )⌋ , i ′ Paccu −0.5·log2

νj := νj′−1 · 2⌊

(max1≦h≦m c̄hj )⌋ .

(4) (5)

4

Journal Title XX(X)

′ Here, Paccu ∈ F32 is a precomputed single-precision

constant defined as ′ Paccu := single▽ (log2 (P − 1)/2 − 0.5) .

The log2 (·) in (4) and (5) is evaluated using the singleprecision fast base-2 logarithm function log2f provided by the CUDA Math API on NVIDIA GPUs (NVIDIA 2026) or the HIP Math API on AMD GPUs (AMD 2026). The log2f

Algorithm 1 Computation of µi in fast mode (row i of A)  P k 2 Require: anorm := h=1 aih computed with roundingup Require: amax := max1≤h≤k |aih | ′ Require: Pfast := single▽ (log2 (P − 1)/2 − 1.5) ∈ F32 Require: δ := single△ (0.5/(1 − 4u32 )) ∈ F32 Ensure: log2 (µi ) ∈ Z 1: s := ⌊log2 (anorm )⌋ 2: t := single△ (δ · ( log2f (anorm /2s ) + s)) ′ 3: return ⌊Pfast − max (1.0, t)⌋ − ⌊log2 (amax )⌋

in the CUDA Math API has an absolute error of at most 4u32 for inputs x ∈ [0.5, 2] and an absolute error of at most 2 ULP (unit in the last place) for other inputs (NVIDIA 2026). For the HIP Math API, the absolute error is at most 1 ULP for x ∈ [10−6 , 106 ] (AMD 2026). According to Uchino et al. (2026a), normalizing x > 0 to x′ := x/2⌊log2 (x)⌋ ∈ [1, 2)

g (·) µ and ν are defined by the following formulas, where log 2 and δ are as given in (7) and (6), respectively: ′ g µi := 2⌊Pfast −max(1,δ·log2 (

Pk

2 h=1 aih

))⌋−⌊log2 (max1≤h≤k |aih |)⌋ ,

yields log2f(x′ ) ≤ log2 (x′ )/(1 − 4u32 ). Thus, defining  δ := single△

0.5 1 − 4u32

 ,

(6)

and evaluating log2 (x) as g (x) := single log 2 △



log2f

(11) νj := 2⌊

′ Pfast −max

2 g (Pk (1,δ·log 2 h=1 bhj ))⌋−⌊log2 (max1≤h≤k |bhj |)⌋

(12) ′ Here, Pfast ∈ F32 is a precomputed single-precision constant

x





2⌊log2 (x)⌋

 + ⌊log2 (x)⌋ ,

defined as ′ Pfast := single▽ (log2 (P − 1)/2 − 1.5) .

(7)

(4) and (5) can be computed as follows: Pk

2 h=1 aih

and

Pk

2 h=1 bhj ,

(max1≦h≦n c̄ih )⌋ ,

(8)

respectively, are evaluated in single or double precision with

′ g Paccu −δ·log 2

(max1≦h≦m c̄hj )⌋ .

(9)

rounding-up; these sums, together with the row- and column-

µi := µ′−1 · 2⌊ i νj := νj′−1 · 2⌊

In (11) and (12), the sums

′ g Paccu −δ·log 2

wise maximum values, dominate the cost of computing µ and ′

With µ and ν defined as above, A and B

′

satisfy

detailed computation of µi is summarized in Algorithm 1.

condition (1) as follows: 2

k X

|a′ih |

b′hj

h=1

ν, requiring O(mk) and O(kn) operations, respectively. The

≤ 2µi

k X

Accurate mode provides a tighter estimate, allowing a

! |aih ||bhj | νj

larger bit width for A′ and B ′ than that in fast mode.

h=1

However, it incurs an additional computational cost for the

′−1 ≤ 2µi µ′−1 < P. i c̄ij νj νj

auxiliary INT8 GEMM. Fast mode avoids this overhead

The dominant cost in computing µ and ν is the auxiliary INT8 GEMM, which has a time complexity of O(mnk). A detailed error analysis of accurate mode is provided by Uchino et al. (2026c).

by using the Cauchy–Schwarz inequality, which is less expensive to evaluate, resulting in shorter computation time than that of accurate mode. However, since the estimate is looser, the bit width of A′ and B ′ may be smaller than that in accurate mode, potentially reducing accuracy.

In fast mode, the estimation is performed using the

3

Limitation of Fast Mode Scaling

Cauchy–Schwarz inequality. That is, µ and ν are chosen to

Here, we identify a potential limitation in the fast mode

satisfy the following:

scaling formulas (11). We analytically show that (11) lacks

2

k X

scale invariance and may cause CRT recovery failure. |a′ih |

b′hj

h=1

Prepared using sagej.cls

≤ 2µi ∥ai,: ∥2 ∥b:,j ∥2 νj < P.

(10)

We then present the results of numerical experiments to demonstrate that these issues occur in practice.

.

Kawakami and Takahashi

3.1

5

Theoretical Analysis

In the scaling step (Step 2) of Ozaki scheme II, the i-th row of the scaled integer matrix A′ has an effective precision of    log2 (µi ) + log2 max |aih | +1 1≤h≤k

bits. When fast mode is used, this precision is evaluated as follows:    log2 (µi ) + log2 max |aih | +1 1≤h≤k $ !!% k X ′ 2 g = P − max 1, δ · log a 2

fast

Figure 1. (A′ B ′ )ij in fast mode versus scalar α = 2s , where A = α and B = αB̂ with âij = b̂ij = 1 (all-ones matrices), and N = 20 moduli. The horizontal black line marks the CRT recovery threshold P/2. The dark-shaded region indicates where the CRT uniqueness condition is violated ((A′ B ′ )ij > P/2), and the light-shaded region indicates where it is satisfied. Thin and thick lines correspond to k = 1024 and k = 16384, respectively, where k is the inner dimension of the matrix product.

ih

h=1

      − log2 max |aih | + log2 max |aih | +1 1≤h≤k 1≤h≤k $ !!% k X ′ 2 g = P − max 1, δ · log a + 1. 2

fast

ih

h=1

Let A = α for some scalar α = 2s (s ∈ Z) and matrix  ∈ Fm×k . Ideally, the scaled integer matrix A′ =

the effective precision does not depend on α and is fixed

trunc(diag(µ)A) should be invariant with respect to α; that

′ − 1⌋ + 1 bits; however, A′ B ′ may exceed P/2, at ⌊Pfast

is, its entries and bit width should be independent of α.

causing the CRT recovery to fail. For example, suppose

This is because, to maximize accuracy, µ should be chosen

âij = 1 for all i, j and α ≤ k −1/2 . Then:

to maximize the bit width of the scaled integer matrix A

′

subject to the CRT condition (1). More precisely, for α0 =

δ·

2s0 and α1 = 2s1 (s0 , s1 ∈ Z), let A0 = α0 Â and A1 =

g log 2

k X

! â2ih

! + 2 log2 (α)

h=1

α1 Â, and let µ0 and µ1 be the scaling factors determined

= δ · (log2 (k) + 2 log2 (α))

by (11) for A0 and A1 , respectively. Then, ideally A′0 =

≤ δ · (log2 (k) − log2 (k)) = 0 < 1,

trunc(diag(µ0 )A0 ) and A′1 = trunc(diag(µ1 )A1 ) should be equal. In fast mode, however, this scale invariance does

so the scaled integer matrix A′ = trunc(diag(µ)A) does not

not hold.

′ depend on α, and its bit is fixed at ⌊Pfast − 1⌋ + 1  width 

For A = αÂ, the effective precision of the i-th row of A′ then becomes: $ ′ Pfast − max

bits. Similarly, let B̂ = b̂ij ∈ Fk×n with b̂ij = 1 and B =

αB̂; then, for α ≤ k −1/2 , B ′ = trunc(B diag(ν)) is also 1, δ ·

g log 2

k X

! â2ih

!!% + 2 log2 (α)

α-independent with the same fixed bit width. In this case, since max1≤h≤k |aih | = max1≤h≤k |bhj | = α for all i, j,

h=1

the scaling factors µ and ν obtained from A and B are given

+ 1.

(13)

When

by ′ ′ µi = νj = 2⌊Pfast −1⌋−log2 (α) = 2⌊Pfast −1⌋ · α−1 .

1<δ·

g log 2

k X

! â2ih

! + 2 log2 (α) ,

h=1

this expression depends on α, and the effective precision of A′ decreases as α increases; hence, the fast mode scaling lacks scale invariance.

As all entries of A and B are α, all entries of A′ = trunc(diag(µ)A) and B ′ = trunc(B diag(ν)) ′ P −1 become 2⌊ fast ⌋ . Consequently, each entry of A′ B ′ is ′

′ k · 22⌊Pfast −1⌋ . Since P and Pfast do not depend on k, for

sufficiently large k, the entries of A′ B ′ exceed P/2 and CRT recovery fails.

On the other hand, when

Here, for k ∈ {1024, 16384}, we set âij = 1 and b̂ij = 1≥δ·

g log 2

k X h=1

Prepared using sagej.cls

! â2ih

! + 2 log2 (α) ,

1, and evaluate (A′ B ′ )ij for input matrices A = α and B = αB̂ with α = 2s (s = −10, −9, . . . , 10); the results

6

Journal Title XX(X)

Figure 2. Maximum relative error of DGEMM with respect to double-double precision versus scalar α = 2s for cuBLAS, OS II-fast, and OS II-accu, where A and B are the input matrices given by A = α and B = αB̂ , with  and B̂ random matrices generated via (15). Each column corresponds to a different value of ϕ controlling the spread of element magnitudes. Matrix dimensions are m = n = k ∈ {1024, 16384} and N = 20 moduli. Thin and thick lines correspond to m = n = k = 1024 and m = n = k = 16384, respectively.

are shown in Fig. 1. We use N = 20 moduli as defined in

• OS II-fast: DGEMM computed by Ozaki scheme II with fast mode implemented in GEMMul8.

GEMMul8:

• OS II-accu: DGEMM computed by Ozaki scheme II

{pi }N i=1 = {256, 255, 253, 251, 247,

with accurate mode implemented in GEMMul8.

241, 239, 233, 229, 227, 223, 217, 211, 199, 197,

The number of moduli used in OS II-fast and OS II-accu

193, 191, 181, 179, 173}.

(14)

is set to N = 20, with the moduli {pi } as defined in (14). Note that for the scaling procedure in GEMMul8 v3.0.4 (the

′

′

As shown in Fig. 1, (A B )ij varies with α and exceeds the −2

CRT recovery threshold P/2 in the region α ≤ 2

.

latest version at the time of writing), fast mode differs from Algorithm 1. However, throughout all evaluations in this paper, to discuss the results based on Uchino et al. (2025),

3.2

we reverted the part to the procedure in Algorithm 1, which

Numerical Demonstration

corresponds to the implementation in GEMMul8 v2.0.19.

We verify the above behavior using numerical examples. All

We confirmed that the fast mode in GEMMul8 v3.0.4

experiments were conducted on an NVIDIA GH200 GPU

exhibits worse accuracy than that of the fast mode based

with CUDA Toolkit 13.2.

on Algorithm 1 and that it shows accuracy characteristics

We generate double-precision random matrices  =  

(âij ) ∈ Fm×k and B̂ = b̂ij ∈ Fk×n according to 64 64

similar to those of OS II-fast, as described below. The modified source code used in this evaluation is publicly available (Kawakami 2026).

âij , b̂ij = (rand −0.5) · exp (ϕ · randn) ,

(15)

The results are shown in Fig. 2. Both cuBLAS and OS IIaccu maintained constant accuracy across all values of α for

where rand denotes a uniform random variable on [0, 1),

each ϕ value. In contrast, for OS II-fast, accuracy decreased

randn denotes a standard normal random variable, and ϕ is a parameter that controls the spread of the entry magnitudes.

as α increased. This is consistent with the analysis above; in the regime where

We evaluate ϕ ∈ {0.5, 1, 2, 3, 4}. Fixing the matrices  and B̂, we generate the inputs A = α and B = αB̂ for α =

1<δ·

2s , s ∈ {−10, −9, . . . , 10}. The matrix dimensions are set

g log 2

k X

! â2ih

! + 2 log2 (α)

h=1

to m = n = k ∈ {1024, 16384}. We evaluate the maximum relative error with respect to the result computed in doubledouble precision. We compare the following methods:

or 1<δ·

g log 2

k X

! b̂2hj

! + 2 log2 (α) ,

h=1

the effective bit width of A′ or B ′ depends on α and • cuBLAS: DGEMM computed by NVIDIA cuBLAS. Prepared using sagej.cls

decreases as α increases. Furthermore, for ϕ = 0.5, the

Kawakami and Takahashi

7

Figure 3. Maximum relative error of DGEMM with respect to double-double precision versus scalar α = 2s for cuBLAS, OS II-fast, and OS II-accu, where A and B are the input matrices given by A = α and B = αB̂ , with âij = b̂ij = 1 (all-ones matrices), m = n = k ∈ {1024, 16384}, and N = 20 moduli. cuBLAS achieves a maximum relative error of exactly 0 in all cases.

maximum relative error reaches approximately 1 in the

Taking the base-2 logarithm of the bound on µi in (16) gives

region α ≤ 2−3 for k = 1024 and α ≤ 2−5 for k = 16384, indicating a complete breakdown of the computation. This is

s log2 (µi ) ≤ log2

likely because the entries of the scaled integer matrix product A′ B ′ exceed P/2, causing CRT recovery to fail.

P −1 1 · 2 ∥ai,: ∥2

!

= log2 (P − 1)/2 − 0.5 − 0.5 · log2

In addition, we evaluated the all-ones case (âij = b̂ij = 1). The results are shown in Fig. 3. For cuBLAS, the maximum relative error was exactly 0 in all cases. For both k = 1024 and k = 16384, OS II-accu maintained constant error, whereas OS II-fast reached a maximum relative error of approximately 1 in the region α ≤ 2−2 , indicating a complete breakdown of the computation. This is consistent −2

with the results shown in Fig. 1; in the region α ≤ 2

,

the entries of the scaled integer matrix product A′ B ′ exceed

k X

! a2ih

.

h=1

We therefore define µi as the largest power of two satisfying this bound: ′ k 2 g µi := 2⌊Pprop −0.5·(log2 ( h=1 aih )+4u32 )⌋ s P −1 1 ≤ · , 2 ∥ai,: ∥2

P

(17)

and analogously for νj :

P/2, causing CRT recovery to fail. k 2 ′ g νj := 2⌊Pprop −0.5·(log2 ( h=1 bhj )+4u32 )⌋ s 1 P −1 ≤ · . 2 ∥b:,j ∥2

P

These observations indicate that the fast mode scaling lacks scale invariance. In this paper, we propose revised definitions of µ and ν that still rely on the Cauchy–Schwarz

(18)

inequality and introduce no additional computational cost, ′ Here, Pprop is defined as

while ensuring scale invariance.

′ Pprop := single▽ (log2 (P − 1)/2 − 0.5) .

4

Pk Pk The sums h=1 a2ih and h=1 b2hj are evaluated in single or double precision with rounding-up. The existing fast and

Proposed Method

We define µ and ν directly from the Cauchy–Schwarz bound

accurate modes of GEMMul8 bound the error of

(10). With the requirement of

for x ∈ [1, 2) as

log2f(x) ≤ log2 (x)/(1 − 4u32 ), treating

the error as multiplicative. However, since 2µi ∥ai,: ∥2 ∥b:,j ∥2 νj ≤ P − 1 < P

log2f has an

absolute error of at most 4u32 for x ∈ [1, 2), the correct bound is | log2f(x) − log2 (x)| ≤ 4u32 , which gives

to hold for all i, j, a sufficient condition is s µi ≤

P −1 1 · , 2 ∥ai,: ∥2

Prepared using sagej.cls

s νj ≤

log2 (x) ≤ log2f(x) + 4u32 . P −1 1 · . (16) 2 ∥b:,j ∥2

log2f

8

Journal Title XX(X)

Algorithm 2 Computation of µi in proposed method (row i of A) P  k 2 Require: anorm := h=1 aih computed with roundingup ′ := single▽ (log2 (P − 1)/2 − 0.5) ∈ F32 Require: Pprop Ensure: log2 (µi ) ∈ Z 1: s := ⌊log2 (anorm )⌋ 2: t := single△ (0.5 · log2f (anorm /2s ) + 2u32 )
′ 3: w := single▽ Pprop − 0.5s 4: return ⌊w − t⌋

of the scaled integer matrix product A′ B ′ are always bounded by P/2, ensuring that CRT recovery never fails. The time complexity of computing µ and ν is O(mk) and O(kn), respectively, which is the same order as that of fast mode, introducing no additional computational cost and remaining faster than accurate mode.

5

Evaluation

We compare the following implementations for both DGEMM and SGEMM in terms of numerical accuracy and computational performance.

Accordingly, (17) and (18) use

• cuBLAS: computed by NVIDIA cuBLAS.

  log2 (x) ≤ log2f x/2⌊log2 (x)⌋ + ⌊log2 (x)⌋ + 4u32

• OS II-fast: computed by Ozaki scheme II with fast

g (x) + 4u32 . = log 2

mode implemented in GEMMul8.

The detailed computation of µi is summarized in Algo-

• OS II-accu: computed by Ozaki scheme II with

rithm 2.

accurate mode implemented in GEMMul8.

Using this µi , the effective bit width of the i-th row of the

• OS II-prop: computed by Ozaki scheme II with fast

scaled integer matrix A′ = trunc(diag(µ)A) is:

mode implemented in GEMMul8, with the scaling formula replaced by the proposed one.

   log2 (µi ) + log2 max |aih | +1 1≤h≤k $ ! !% k X ′ 2 g = P − 0.5 · log a + 4u32 2

prop

The source code of OS II-prop is available in (Kawakami 2026). As noted in Section 3.2, the scaling procedure of

ih

OS II-fast was reverted to the v2.0.19 implementation (the

h=1





+ log2

 max |aih |

procedure in Algorithm 1) throughout all evaluations in this

+ 1.

1≤h≤k

section. All experiments were conducted on an NVIDIA

As in the analysis in Section 3.1, by setting A = α (α = 2s , s ∈ Z), the effective bit width becomes: $ ′ Pprop − 0.5 ·

g log 2

k X

5.1

! â2ih

!% + 2 log2 (α) + 4u32

h=1





  + log2 max |âih | + log2 (α) + 1 1≤h≤k $ ! !% k X ′ 2 g = P − 0.5 · log â + 4u32 − log (α) 2

prop

ih

2

h=1





 + log2 max |âih | + log2 (α) + 1 1≤h≤k $ ! !% k X ′ 2 g = P − 0.5 · log â + 4u32 2

prop

GH200 GPU with CUDA Toolkit 13.2, as in Section 3.2.

ih

h=1

   + log2 max |âih | + 1. 1≤h≤k

Accuracy

Input matrices A ∈ Fm×k and B ∈ Fk×n were generated using (15) and the matrix product was computed by each method. The maximum relative error with respect to a double-double precision reference was measured. The matrix dimensions were set to m = n = k ∈ {1024, 16384}. The parameter ϕ in (15) was set to ϕ ∈ {0.5, 1, 2, 3, 4} for DGEMM and ϕ ∈ {0, 0.5, 1, 1.5} for SGEMM. For OS IIfast, OS II-accu, and OS II-prop, the number of moduli N was varied over 9 ≤ N ≤ 20 for DGEMM and 3 ≤ N ≤ 12 for SGEMM. The results are shown in Fig. 4. For all values of ϕ, while the maximum relative errors of OS II-accu, OS IIfast, and OS II-prop exceed that of cuBLAS for small N ,

Here, the second equality uses ⌊x ± s⌋ = ⌊x⌋ ± s, which s

they tend to decrease as N increases. A comparison of

holds since log2 (α) = log2 (2 ) = s ∈ Z. As log2 (α)

the maximum relative errors of OS II-fast and OS II-accu

cancels out in this way, the effective bit width of the

indicates that the gap is small for small values of ϕ but

scaled integer matrix is constant and independent of α.

grows larger as ϕ increases, for both DGEMM and SGEMM.

Furthermore, since µi is designed to satisfy (16), the entries

In particular, for DGEMM with the ill-conditioned cases

Prepared using sagej.cls

Kawakami and Takahashi

9

Figure 4. Maximum relative error with respect to double-double precision of DGEMM (top) and SGEMM (bottom) versus number of moduli N for cuBLAS, OS II-fast, OS II-accu, and OS II-prop, using random matrices with m = n = k and varying ϕ. Thin and thick lines correspond to m = n = k = 1024 and m = n = k = 16384, respectively.

ϕ ∈ {3, 4}, the maximum relative error of OS II-fast exceeds

were reached; timed iterations then continued until at least

that of cuBLAS. In contrast, OS II-prop maintains accuracy

5 iterations had run and either 12 seconds elapsed or 100

comparable to that of OS II-accu regardless of the value of

iterations were reached. The median execution time Tmed

ϕ. Furthermore, OS II-prop achieves higher accuracy than

was used to compute throughput in FLOPS as follows:

that of cuBLAS with a sufficient number of moduli, both for FLOPS :=

DGEMM with ϕ ∈ {3, 4}. Although OS II-prop determines µ and ν using the Cauchy–Schwarz inequality, it maintains accuracy comparable to that of OS II-accu. Although OS II-

2mnk . Tmed

The results are shown in Figs. 5 and 6.

(19) OS II-prop

fast also uses the Cauchy–Schwarz inequality, it exhibits

achieves higher throughput than that of OS II-accu for

accuracy degradation for large ϕ, suggesting that the root

both DGEMM and SGEMM. Since the accuracy evaluation

cause of the degradation is not the use of the Cauchy–

showed that OS II-prop maintains accuracy comparable to

Schwarz inequality itself.

that of OS II-accu, OS II-prop delivers higher throughput at equivalent accuracy.

5.2

Both OS II-prop and OS II-fast achieve similar throughput

Throughput

Input matrices A ∈ F

for ϕ = 0.5 across DGEMM and SGEMM. For ϕ = 4 in m×k

and B ∈ F

k×n

were gener-

DGEMM and ϕ = 1.5 in SGEMM, the two methods are

ated using (15). Matrix multiplication was performed

also comparable for small matrix sizes and for large matrix

with each method and the execution time was mea-

sizes with a large number of moduli. However, for large

sured. The matrix sizes were m = n = k, where k ∈

matrix sizes with a small number of moduli in these cases,

{1024, 2048, 4096, 8192, 16384}. The parameter ϕ in (15)

OS II-fast achieves higher throughput than that of OS II-

was set to ϕ ∈ {0.5, 4} for DGEMM and ϕ ∈ {0.5, 1.5} for

prop. We attribute this difference to the fraction of zero

SGEMM. For OS II-fast, OS II-accu, and OS II-prop, the

elements in the scaled integer matrices. For ϕ = 0.5, the

number of moduli was varied over 9 ≤ N ≤ 20 for DGEMM

spread of element exponents is small. The zero fraction

and 3 ≤ N ≤ 12 for SGEMM. For each method, warm-up

after scaling is therefore similar between OS II-prop and

iterations were performed until at least 3 iterations had run

OS II-fast. With a large number of moduli, the CRT bound

and either 3 seconds of total time elapsed or 100 iterations

P/2 grows with N , allowing a larger range of values to be

Prepared using sagej.cls

10

Journal Title XX(X)

Figure 5. Throughput (TFLOPS) of DGEMM for cuBLAS, OS II-fast, OS II-accu, and OS II-prop versus number of moduli N for random matrices with ϕ = 0.5 (top) and ϕ = 4 (bottom), for varying m = n = k .

Figure 6. Throughput (TFLOPS) of SGEMM for cuBLAS, OS II-fast, OS II-accu, and OS II-prop versus number of moduli N for random matrices with ϕ = 0.5 (top) and ϕ = 1.5 (bottom), for varying m = n = k .

represented after scaling; consequently, the zero fractions of

where the zero fraction of OS II-fast’s scaled integer matrices

the two methods become comparable. With a small number

is substantially higher than that of OS II-prop.

of moduli, however, the representable range is narrower for OS II-fast due to its smaller scaling factor, resulting in a higher fraction of zero elements. This effect is especially pronounced for ϕ = 4 in DGEMM and ϕ = 1.5 in SGEMM,

Prepared using sagej.cls

As noted in the accuracy evaluation, both OS II-prop and OS II-fast exhibit considerably higher error than that of cuBLAS when the number of moduli is small. Therefore, the throughput disadvantage of OS II-prop at small N is not a

Kawakami and Takahashi

11

Figure 7. Accuracy–throughput trade-off for DGEMM (top, ϕ ∈ {0.5, 1, 2, 3, 4}) and SGEMM (bottom, ϕ ∈ {0, 0.5, 1, 1.5}) at m = n = k = 16384. The horizontal axis shows throughput in TFLOPS and the vertical axis shows maximum relative error with respect to double-double precision. Each curve is traced from right to left as the number of moduli N increases. The black dot marks the throughput and accuracy of cuBLAS; the shaded region indicates where a method simultaneously achieves higher accuracy and higher throughput than those of cuBLAS.

practical concern when cuBLAS-level accuracy is required.

The curve of OS II-prop lies consistently to the right

The accuracy–throughput trade-off is analyzed in detail in

of that of OS II-accu for all values of ϕ. This indicates

the next subsection.

that OS II-prop achieves higher throughput than that of OS II-accu at any given accuracy level. OS II-prop therefore dominates OS II-accu in the accuracy–throughput trade-off.

5.3

Trade-off between Accuracy and Performance

To analyze the accuracy–throughput trade-off, we measured additional throughput data for m = n = k = 16384, specifically for DGEMM with ϕ ∈ {1, 2, 3} and SGEMM with ϕ ∈ {0, 1}, supplementing the throughput results for DGEMM with ϕ ∈ {0.5, 4} and SGEMM with ϕ ∈ {0.5, 1.5} shown in Figs. 5 and 6; we combined these results with the accuracy data from Fig. 4.

The curve of OS II-fast shifts upward (toward higher error) as ϕ increases. In particular, for DGEMM with ϕ ∈ {3, 4}, the curve of OS II-fast does not reach the shaded region, meaning that OS II-fast cannot simultaneously outperform cuBLAS in both accuracy and throughput under these conditions. In contrast, OS II-accu and OS II-prop reach the shaded region for all values of ϕ, meaning that they can simultaneously outperform cuBLAS in both accuracy and throughput.

The combined results are shown in Fig. 7, where the

With the exception of SGEMM with ϕ = 0, the curve

horizontal axis represents throughput and the vertical axis

of OS II-prop lies furthest to the right among all methods,

represents the maximum relative error. The black dot

confirming that OS II-prop achieves the highest throughput

indicates the throughput and accuracy of cuBLAS, and the

for any target accuracy level. In the SGEMM case with

shaded region marks the area where a method achieves

ϕ = 0, the spread of element exponents is small and OS II-

both higher accuracy and higher throughput than those of

fast also achieves comparable accuracy; consequently, the

cuBLAS. For each of OS II-fast, OS II-accu, and OS II-prop,

curves of OS II-fast and OS II-prop nearly coincide.

the data points trace a curve from the upper right to the lower left as the number of moduli increases, reflecting the trade-

throughput than that of OS II-accu at equivalent accuracy,

off between accuracy and throughput.

while resolving the accuracy degradation of OS II-fast for

Prepared using sagej.cls

In summary, OS II-prop achieves consistently higher

12

Journal Title XX(X)

large ϕ without sacrificing throughput. The proposed method thus overcomes both the accuracy limitation of OS II-fast and the throughput constraint of OS II-accu, offering a superior balance of accuracy and performance.

Funding This research was supported by JST SPRING (grant number JPMJSP2124) and used the computational resources of Miyabi provided by the Multidisciplinary Cooperative Research Program in the Center for Computational Sciences, University of Tsukuba.

6

Conclusion

In this paper, we proposed an improved fast mode scaling

Supplemental material Not applicable.

for DGEMM and SGEMM via Ozaki scheme II with INT8 matrix engines, achieving better accuracy without sacrificing throughput. We identified a potential limitation in the scaling formula of the existing fast mode. We showed that the formula lacks scale invariance: when the input matrix is scaled by a factor α, the effective bit width of the scaled integer matrix varies with α, leading to accuracy degradation for large α and CRT recovery failure for small α.

References AMD (2026) AMD ROCm documentation: HIP Math API. URL https://rocm.docs.amd.com/projects/HIP/e n/latest/reference/math_api.html. Brower C, Rodriguez Bernabeu S, Hammond J, Gunnels J, Xantheas SS, Ganahl M, Menczer A and Legeza O (2026) Mixed-precision ab initio tensor network state methods adapted

To address this limitation, we proposed a revised scaling formula derived directly from the CRT uniqueness condition via the Cauchy–Schwarz inequality. The proposed formula is scale-invariant by construction and guarantees that the entries of the scaled integer matrix product A′ B ′ never exceed the CRT recovery threshold P/2. Its computational complexity is the same as that of the original fast mode, introducing no additional overhead. Experiments on the NVIDIA GH200 GPU confirmed that the proposed method achieves accuracy comparable to that of accurate mode for all values of ϕ. Its throughput also approaches that of the original fast mode as the number of moduli increases. In the accuracy–throughput trade-off, the proposed method thus overcomes both the accuracy limitation of OS II-fast and the throughput constraint of OS II-accu, offering a superior balance of accuracy and performance. As future work, it would be interesting to extend the proposed scaling formula beyond DGEMM and SGEMM with INT8 matrix engines to complex matrix multiplication (ZGEMM and CGEMM) (Uchino et al. 2026a), as well as to methods that exploit FP8 matrix engines (Uchino 2026). A rigorous error analysis of the proposed method, analogous to that performed for accurate mode in (Uchino et al. 2026c), is also left as future work.

for nvidia blackwell technology via emulated fp64 arithmetic. J. Chem. Theory Comput. DOI:10.1021/acs.jctc.6c00203. Dawson W, Ozaki K, Domke J and Nakajima T (2024) Reducing Numerical Precision Requirements in Quantum Chemistry Calculations. URL https://arxiv.org/abs/2407 .13299. Dongarra J, Gunnels J, Bayraktar H, Haidar A and Ernst D (2024) Hardware Trends Impacting Floating-Point Computations In Scientific Applications. URL https://arxiv.org/ab s/2411.12090. Kawakami S (2026) GEMMul8. URL https://github.c om/kotatsumuri/GEMMul8. Fork of GEMMul8 with improved fast mode scaling. Mukunoki D (2026) DGEMM using FP64 Arithmetic Emulation and FP8 Tensor Cores with Ozaki Scheme.

In: Proc.

Supercomputing Asia and International Conference on High Performance Computing in Asia Pacific Region Workshops. pp. 303–311. DOI:10.1145/3784828.3785017. Mukunoki D, Ozaki K, Ogita T and Imamura T (2020) DGEMM Using Tensor Cores, and Its Accurate and Reproducible Versions. In: Proc. ISC High Performance 2020. pp. 230–248. DOI:10.1007/978-3-030-50743-5 12. NVIDIA (2024) NVIDIA Blackwell Architecture Technical Overview. URL https://resources.nvidia.com /en-us-blackwell-architecture/blackwell-a rchitecture-technical-brief.

Declaration of conflicting interests

NVIDIA (2025) NVIDIA GH200 Grace Hopper Superchip Datasheet. URL https://nvdam.widen.net/s/rr

The authors declared no potential conflicts of interest with respect

gqqnpbz8/grace-datasheet-gh200-grace-hop

to the research, authorship, and publication of this article.

per-superchip-3773000.

Prepared using sagej.cls

Kawakami and Takahashi

NVIDIA (2026) CUDA Programming Guide: 5.5. Floating-Point Computation. URL https://docs.nvidia.com/cuda /cuda-programming-guide/05-appendices/ma thematical-functions.html. Ootomo H, Ozaki K and Yokota R (2024) DGEMM on integer matrix multiplication unit. Int. J. High Perform. Comput. Appl. 38(4): 297–313. DOI:10.1177/10943420241239588. Ozaki K, Ogita T, Oishi S and Rump SM (2012) Errorfree transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numer. Algorithms 59(1): 95–118. DOI:10.1007/s11075-011-9478-1. Ozaki K, Uchino Y and Imamura T (2025) Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique. URL ht tps://arxiv.org/abs/2504.08009. Uchino Y (2026) GEMMul8. URL https://github.com/R IKEN-RCCS/GEMMul8. Uchino Y, Ma Q, Imamura T, Ozaki K and Gutsche PL (2026a) Emulation of Complex Matrix Multiplication based on the Chinese Remainder Theorem. In: Proc. ISC High Performance 2026. pp. 1–12. DOI:10.23919/ISC.2026.11520500. Uchino Y, Ozaki K and Imamura T (2025) High-Performance and Power-Efficient Emulation of Matrix Multiplication using INT8 Matrix Engines. In: Proc. SC ’25 Workshops of the International Conference for High Performance Computing, Networking, Storage and Analysis. p. 1824–1831. DOI:10.1 145/3731599.3767539. Uchino Y, Ozaki K and Imamura T (2026b) Double-Precision Matrix Multiplication Emulation via Ozaki-II Scheme with FP8 Quantization. URL https://arxiv.org/abs/26 03.10634. Uchino Y, Ozaki K and Imamura T (2026c) Error Analysis of Matrix Multiplication Emulation Using Ozaki-II Scheme. URL https://arxiv.org/abs/2602.02549.

Prepared using sagej.cls

13