You Only Convert Twice: Ozaki Scheme II for Chained Tensor Mode Products

arXiv:2610.05059v1 [cs.DC] 4 Oct 2026

YOU O NLY C ONVERT T WICE : O ZAKI S CHEME II FOR C HAINED T ENSOR M ODE P RODUCTS

Xuanzhengbo Ren Information Technology Center The University of Tokyo Kashiwanoha, Chiba, 2770882 Japan xuanzhengbo.ren@cc.u-tokyo.ac.jp

October 6, 2026

A BSTRACT Ozaki scheme II enables high-precision matrix multiplication using low-precision Tensor Cores through an integer residue representation. Applying it to successive tensor mode products offers a way to exploit the low-precision performance of recent GPUs. This work proposes OzII-RescaleBE, a method that performs rescaling and base extension directly on residue representations. It enables persistent execution of chained tensor mode products, keeping intermediate tensors in residue form. Conversion to residues is required only at the beginning of the chain, followed by reconstruction at the end. An implementation using CuTe DSL and INT8 Tensor Cores is provided for verification and evaluation. The method is compared with per-mode composition of Ozaki scheme II (GEMMul8Composed) and a dense Kronecker-product formulation (GEMMul8-KRON), both implemented using the GEMMul8 library. Across three datasets with different input distributions, OzII-RescaleBE achieves normwise relative errors ranging from 2.5 × 10−16 to 2.6 × 10−15 for chain depths d = 3 through 7. The corresponding errors range from 10−16 to 10−15 for GEMMul8-KRON and from 1 × 10−16 to 3 × 10−16 for GEMMul8-Composed and the FP64 chain. At d = 8, the errors of OzII-RescaleBE increase to between 7 × 10−14 and 6 × 10−13 . For mode sizes n = 8 to 256, OzII-RescaleBE achieves throughputs ranging from 1.47 to 4.42 GDoF/s. Its throughput is within 12 % of GEMMul8-KRON at n = 8 and 12 and exceeds it at larger tested sizes. Compared with GEMMul8-Composed, OzII-RescaleBE achieves 2.1–60× the throughput for n ≤ 128, but has 6–12 % lower throughput at n = 256. It also uses less device memory than the GEMMul8 baselines in the measured memory comparisons. These results demonstrate that rescaling and base extension enable accurate residue-domain tensor chains while avoiding repeated intermediate conversions to floating point. Keywords tensor mode product · Ozaki scheme II · residue number systems · rescaling and base extension

1

Introduction

The growing demand for artificial intelligence training and inference is shaping GPU architectures around highthroughput, low-precision matrix arithmetic. However, improvements in these capabilities do not necessarily extend to native double-precision (FP64) computation. For example, NVIDIA lists 34 TFLOP/s of conventional FP64 throughput for the H100 SXM, whereas an individual Blackwell Ultra GPU in the HGX B300 is specified at 1.2 TFLOP/s of FP64 throughput and 2.25 PFLOP/s of dense FP16/BF16 Tensor Core throughput [1, 2]. For scientific applications that depend on FP64 accuracy, this disparity motivates algorithms that exploit low-precision matrix units while retaining higher-precision numerical behavior. Mixed-precision methods exploit faster lower-precision arithmetic by combining computations at different precisions, reserving higher precision for operations that are more sensitive to numerical error [3]. Precision emulation provides a complementary approach by constructing higher-precision arithmetic from multiple low-precision operations.

A PREPRINT - O CTOBER 6, 2026

The Ozaki schemes offer a promising foundation for such emulation. Originally developed using error-free transformations of matrix products, they have since been adapted to floating-point and integer Tensor Cores [4, 5, 6]. Among these approaches, Ozaki scheme II uses an integer residue representation to perform matrix multiplication on low-precision integer hardware [7]. It first scales and truncates the floating-point operands to integers, then converts them to residues with respect to several pairwise coprime moduli. This conversion process is referred to as encoding in this paper. Matrix products are then computed independently for each modulus, and the results are reconstructed using the Chinese remainder theorem (CRT) and scaled back to floating point. GEMMul8 provides a public implementation of this approach [8]. We consider its application to chained tensor mode products, a central building block of sum factorization in high-order finite element methods. For tensor-product elements, sum factorization applies one-dimensional operators successively along the tensor modes, avoiding the formation and application of a full multidimensional matrix [9]. A representative three-dimensional operation is Y = X ×1 A ×2 B ×3 C,

X ∈ Rn×n×n ,

A, B, C ∈ Rn×n .

(1)

Here, ×j denotes multiplication along mode j. Each mode product can be expressed as a matrix multiplication through an appropriate unfolding. For the square case in (1), this reduces the arithmetic cost from O(n6 ) for a dense operator application to O(n4 ). The workload therefore appears well suited to GEMM-based emulation, but its contraction dimension is the small element-local size n. Our principal target is a batch of independent small tensors, such as those arising in an element loop, rather than a single large dense matrix product. Applying an emulated GEMM to each mode product requires encoding the floating-point input through scaling and conversion to integer residues, then reconstruct the result. Even when the operator encodings are cached, these conversions must be repeated for the changing tensor. The chain in (1) therefore requires three tensor encodings and reconstructions, or d of each for a chain of depth d. Large matrix multiplications can amortize these conversions over substantial matrix arithmetic, whereas small mode products provide much less arithmetic per converted entry. For a fixed residue count, their product work scales as O(n4 ), while tensor conversion traverses O(n3 ) entries. Repeated encoding and reconstruction can consequently consume a substantial fraction of execution time. Fusion within an emulated GEMM reduces intermediate traffic and launch overhead [10], but composing such GEMMs still requires a floating-point intermediate between mode products. This motivates keeping the intermediate tensor in residue form, with encoding only at the beginning of the chain and reconstruction only at the end. To avoid repeated conversion, we can keep the intermediate tensor in residue form and reindex its residue planes for successive mode products without reconstructing an FP64 tensor between them. However, the represented integer can grow as the chain progresses, even though each stored residue remains small. Ozaki scheme II already discusses products of more than two matrices by allocating Qs a range budget for the complete product [7]. For a fixed residue basis containing s moduli, with product P = ℓ=1 pℓ , a sufficient condition for unambiguous signed reconstruction of an integer v is 2|v| < P . The available INT8 modulus pool contains 49 pairwise coprime moduli whose product is approximately 2342 . Even if all these moduli were selected, the centered reconstruction range would accommodate integer magnitudes only up to approximately 2341 . A chain of tensor mode products can require a much larger range, because each mode product multiplies the current integer entries by another encoded operator and sums the products along one mode. For an retained significand bits b, also termed mantissa bits in Ozaki schemes related work, suppose that the encoded input and all effective integer operators have entries bounded in magnitude by 2b . Each mode product can then increase the largest magnitude by a factor of at most n2b . Applying this estimate over d mode products gives ∥X̄d ∥max ≤ nd 2(d+1)b ,

(2)

where ∥X̄d ∥max is the largest absolute entry after the final product. This conservative estimate follows from the triangle inequality. Keeping the bound below P/2 is sufficient to guarantee unambiguous reconstruction. For example, let b = 53 keeping full FP64 significand, at n = 8, d = 6 the bound reaches 2389 , which exceeds the centered range of the entire modulus pool. The actual integer magnitudes can be smaller because this bound uses only entrywise magnitude limits and does not account for cancellation. Nevertheless, the bound can no longer certify reconstruction within the available range. Enlarging the basis also increases both residue storage and the number of residue GEMMs required for each mode product. These observations lead to a central question: How can a chain of tensor mode products reuse an integer residue representation while controlling the magnitude of the represented integers and retaining the tensor-product structure? To address this question, we draw on established techniques in residue-number arithmetic and approximate homomorphic encryption [11, 12, 13] and propose a persistent execution method that combines rescaling and base extension between 2

A PREPRINT - O CTOBER 6, 2026

mode products. Rescaling divides the represented integer by selected moduli, with rounding computed directly from its residues. Base extension then computes residues of the rescaled integer for the removed moduli, restoring the working basis for the next product. Together, these operations allow the intermediate tensor to remain in residue form throughout the chain, with encoding only at entry and reconstruction only at exit. 1. A rescaling and base-extension algorithm for persistent Ozaki scheme II execution of chained tensor mode products. The algorithm combines operator-dependent range planning, separable scale handling, and mixedradix base extension to reuse a fixed working residue basis across the chain. 2. A GPU implementation using CuTe DSL and INT8 tensor cores for verification and evaluation. 3. An experimental evaluation of numerical accuracy, execution time, and resident device memory, including comparisons with GEMMul8-based implementations that apply successive mode products or a cached dense operator [8]. The remainder of the paper is organized as follows. Section 2 introduces tensor mode products, sum factorization, Ozaki scheme II, and the residue arithmetic used in this work. Section 3 discusses related work. Section 4 explains how rescaling and base extension control integer magnitude and manage the available representation range. Section 5 applies these operations to persistent Ozaki scheme II execution of tensor mode products and describes the GPU implementation. Section 6 evaluates numerical accuracy, execution time, and memory usage, including comparisons with GEMMul8. Section 7 concludes the paper.

2

Background

2.1

Tensor Mode Products and Sum-factorization

A tensor mode product applies a matrix along one index of a tensor. For a third-order tensor X ∈ Rn1 ×n2 ×n3 and a matrix A ∈ Rm1 ×n1 , the mode-1 product Z = X ×1 A has entries Zajk =

n1 X

Z ∈ Rm1 ×n2 ×n3 .

Aai Xijk ,

(3)

i=1

Products along modes 2 and 3 are defined in the same way. A mode-r fiber is the vector obtained by fixing every tensor index except the index of mode r. The same operator matrix acts on every fiber of that mode [14]. Figure 1 illustrates three successive mode products for square operators A, B, C ∈ Rn×n and an input X ∈ Rn×n×n . We denote the intermediate tensors by T (1) = X ×1 A,

T (2) = T (1) ×2 B,

Y = T (2) ×3 C.

(4)

The colored fibers in panel (a) indicate the direction of each product. Panel (b) gives its matrix form. The mode-r unfolding X(r) arranges the mode-r fibers as matrix columns. With consistent column ordering, the relation between a mode product and its unfolding is Z = X ×r A ⇐⇒ Z(r) = AX(r) . (5) In the square case, X(r) has dimensions n × n2 . Each mode product is therefore a GEMM with contraction dimension n and n2 fiber columns. The next mode product unfolds the resulting intermediate tensor along a different mode. This is a mathematical reindexing and need not require a separate data copy in an implementation [14]. These mode products form the contraction kernels of sum factorization in high-order finite element methods. On a tensor-product hexahedral element, a basis function is a product of three one-dimensional basis functions. Evaluating a field at tensor-product quadrature points can therefore be performed by three successive one-dimensional interpolation operations. For example, if Vai = ℓi (ξa ) is the one-dimensional interpolation matrix, evaluation of a coefficient tensor U takes the form Uq = U ×1 V ×2 V ×3 V. (6) Here ℓi is a basis function and ξa is a quadrature point. The matrix V may be rectangular when the number of quadrature points differs from the number of basis functions. Differentiation uses the corresponding derivative matrix along the required direction. Mass and stiffness operator applications combine these contractions with quadrature weights, geometric factors, and transpose operations [9]. The same three-direction operation can also be expressed using the Kronecker product of the one-dimensional operator matrices, with their order determined by the tensor vectorization convention [14]. This assembles the directional 3

A PREPRINT - O CTOBER 6, 2026

(a) Tensor modes Mode 2

Mode 1

B

A

(b) Successive GEMMs

A ×

X(1)

B ×

(1)

T (2)

C ×

(2)

=

T (1) (1)

reindex

=

T (2) (2)

Mode 3 reindex

C X

n×n n×n×n

T (3)

=

n × n²

Y(3) n × n²

Contraction length n

Figure 1: Three mode products and their GEMM representations. Panel (a) highlights one fiber in each mode of an n × n × n tensor. Operators A, B, and C act on all fibers of modes 1, 2, and 3, respectively. Panel (b) shows the successive unfolded products using the intermediate tensors. Each GEMM applies an n × n matrix to n2 fibers of length n. operators into a single matrix that acts on the vectorized tensor, applying all three directions in one matrix operation. It is the assembled alternative to sum factorization, which applies the factors successively without forming the multidimensional matrix. For square factors, the assembled matrix has dimensions n3 × n3 and requires O(n6 ) storage and arithmetic for a dense application. Sum factorization instead costs O(n4 ) arithmetic operations and stores only O(n2 ) operator entries, together with O(n3 ) tensor intermediates. 2.2

The Ozaki Scheme II

The Ozaki scheme constructs higher-precision matrix products by splitting the operands into lower-precision components and combining their products [4]. This approach has been implemented using both floating-point and integer Tensor Cores [5, 6]. Ozaki scheme II follows a different approach based on integer residues [7]. After scaling and truncating the operands to integers, it represents them by residues with respect to several pairwise coprime moduli. The resulting residue matrices can be multiplied independently using INT8 Tensor Cores, and their products are combined through the Chinese remainder theorem (CRT) to reconstruct the integer result. Consider a matrix multiplication H = F G, where F ∈ Rm×k and G ∈ Rk×t . Ozaki scheme II computes an approximation to H in four steps. First, the operands are scaled to retain their significant information in an integer representation. Positive diagonal matrices DL and DR , with power-of-two entries, scale the rows of F and the columns of G, respectively. The scaled entries are then truncated toward zero to obtain F̄ = trunc(DL F ),

Ḡ = trunc(GDR ).

(7)

The choice of scaling determines both the magnitudes of these integers and the precision retained from the original operands. Second, these integers are converted into residue matrices that fit the low-precision multiplication units. A basis P = {p1 , . . . , ps } of pairwise coprime moduli is selected, and each operand is reduced modulo every modulus. For INT8 computation, the moduli are at most 256, with a typical basis beginning with 256, 255, 253, . . .. The residues are stored as signed representatives centered around zero, so they lie within the INT8 range [−128, 127]. This produces s residue matrices for each operand, even when the encoded integers themselves are too large to fit in INT8. Third, each residue matrix of F̄ is multiplied by the residue matrix of Ḡ associated with the same modulus. These products use INT8 Tensor Cores with INT32 accumulation. Reducing the accumulated entries modulo that modulus gives the residues of the integer product. The conversion and multiplication steps can therefore be written as   Rℓ = (F̄ mod pℓ )(Ḡ mod pℓ ) mod pℓ , ℓ = 1, . . . , s. (8) Here, each modular reduction is applied entrywise using centered representatives. The products for different moduli are independent. 4

A PREPRINT - O CTOBER 6, 2026

Qs Finally, the output residue matrices are combined using CRT. Let P = ℓ=1 pℓ and Pℓ = P/pℓ . Pairwise coprimality ensures that Pℓ has a multiplicative inverse modulo pℓ , allowing reconstruction in the form ! s X
H̄ = Rℓ Pℓ Pℓ−1 mod pℓ mod P. (9) ℓ=1

The factor Pℓ (Pℓ−1 mod pℓ ) multiplying Rℓ is congruent to one modulo pℓ and zero modulo every other modulus in the

basis. The sum therefore reproduces each output residue matrix when reduced modulo its corresponding modulus. The reconstructed matrix therefore matches all the output residues. The final reduction selects the centered representative modulo P , which recovers H̄ = F̄ Ḡ when the intended integer entries lie within that reconstruction interval. Removing the original scaling then gives b = D−1 H̄D−1 ≈ F G. H (10) L R 2.3

Residue Arithmetic Used in This Work

The residue operations used in this work build on residue number system (RNS) arithmetic. An RNS represents an integer by its remainders with respect to a basis P of pairwise coprime moduli [15]. This representation allows addition and multiplication to be carried out independently on smaller residues, without carry propagation between components. Ozaki scheme II uses this property to compute matrix products Q through low-precision integer operations. The residue tuple determines an integer modulo the basis product P = p∈P p. Recovering the intended integer therefore also requires a known reconstruction interval. To store these residues as signed low-precision integers, we use centered representatives. For an integer v and modulus p ≥ 2, first compute the nonnegative remainder r = v − p⌊v/p⌋. This remainder lies in {0, . . . , p − 1}. Subtracting p from remainders in the upper part of this range gives the centered residue  r, r < ⌈p/2⌉, [v]p = (11) r − p, r ≥ ⌈p/2⌉. The resulting representative satisfies −⌊p/2⌋ ≤ [v]p ≤ ⌈p/2⌉ − 1.

(12)

For example, the nonnegative remainder 5 modulo 7 becomes 5 − 7 = −2 in the centered interval [−3, 3]. For moduli at most 256, this interval lies within the INT8 range [−128, 127], allowing the residues to serve as operands in the integer matrix products of Ozaki scheme II [7]. Changing to centered representatives preserves the rules of modular addition and multiplication. In particular, reducing the operands before an operation gives the same residue as reducing its result, [uv]p = [[u]p [v]p ]p ,

[u + v]p = [[u]p + [v]p ]p .

(13)

These identities also hold for sums of products, which explains why tensor contractions can be computed independently for each modulus. Since reduction acts on individual entries, reindexing a tensor also preserves its residue representation. The same residue planes can therefore be used for successive tensor modes. Although modular reduction keeps each stored residue small, it does not reduce the magnitude of the integer represented by the complete tuple. Managing that magnitude requires a different operation. We use rescaling to replace the represented integer by a rounded quotient, followed by base extension to represent that quotient on a larger basis. 2.3.1

Rescaling

Rescaling reduces the magnitude of the represented integer through rounded division. Related operations are used in approximate homomorphic encryption and its RNS implementations [11, 12]. Suppose that v is represented on P and that a modulus q ∈ P is selected for removal. Its centered residue rq = [v]q provides the remainder needed to form the integer quotient v − rq v′ = . (14) q Because |rq | ≤ q/2, this quotient satisfies v 1 v′ − ≤ . q 2 Thus, subtracting the centered remainder before division produces a nearest-integer approximation to v/q. 5

A PREPRINT - O CTOBER 6, 2026

Moreover, the residues of this quotient can be obtained directly from the residues of v. For every remaining modulus p ̸= q, pairwise coprimality ensures that the modular inverse q −1 mod p exists. Reducing the quotient relation modulo p therefore gives   [v ′ ]p = ([v]p − rq )(q −1 mod p) p . (15) This update applies only to the remaining moduli because q has no multiplicative inverse modulo itself. The quotient is consequently represented on the reduced basis P ′ = P \ {q}, whose product is P ′ = P/q. Rescaling therefore reduces both the integer magnitude and the available representation range. To restore the original basis, we must compute the missing residue [v ′ ]q , which is the role of base extension. 2.3.2

Base Extension

Base extension computes residues of the same represented integer on additional coprime moduli [13]. After the rescaling step above, it restores P by computing [v ′ ]q from the residues on P ′ . The integer remains v ′ , while the additional modulus enlarges the range available for subsequent operations. For this extension to represent the intended quotient, the remaining residues must identify v ′ within a known reconstruction interval. For centered recovery, 2|v ′ | < P ′ is a sufficient condition. Integers differing by a multiple of P ′ have identical residues on P ′ but can have different residues modulo q. Without a known reconstruction interval, the remaining residues alone cannot distinguish those possibilities. Mixed-radix conversion provides one way to compute the additional residue. Order the remaining moduli as P ′ = {e p1 , . . . , pes−1 }. The nonnegative CRT representative u ∈ [0, P ′ ) of the quotient residues has the mixed-radix expansion u=

s−1 X

ai Wi ,

W1 = 1,

Wi =

i=1

Y

pej ,

0 ≤ ai < pei .

(16)

j<i

The digits ai can be obtained by successive modular elimination using Garner’s method [15]. Once these digits are known, the expansion can also be evaluated in nested form using Horner’s method [16]. Starting from the last digit, each step multiplies the current value by the preceding modulus and adds the corresponding digit. Performing this evaluation modulo q gives u mod q. The nonnegative representative u may differ from the intended signed quotient. To account for this, write v ′ = u − δP ′ , where δ = 0 for u < P ′ /2 and δ = 1 otherwise. Reducing this expression modulo q gives # "s−1 X [v ′ ]q = ai (Wi mod q) − δ(P ′ mod q) . (17) i=1

q

Appending this residue restores the basis P. The extended tuple still represents the rescaled integer v ′ , so base extension restores representation capacity but does not recover precision lost during rescaling.

3

Related Work

3.1

Tensor-product Kernels in Finite Element Methods

GPU implementations of tensor-product finite element kernels have progressed from optimizing data movement on conventional GPU cores to exploiting Tensor Cores. Świrydowicz et al. [9] study mass and stiffness operator applications and examine how shared-memory use, intermediate storage, and data movement affect performance. Their work provides a foundation for organizing successive element-local contractions within GPU kernels. More recent work maps these contractions directly to Tensor Cores. Cui [17] investigates tensor-product finite element operations on NVIDIA A100 GPUs, emphasizing direct instruction-level programming and shared-memory layouts that avoid bank conflicts. The study considers both double-precision kernels and a mixed-precision multigrid-preconditioned Poisson solver. Tu et al. [18] extend FP64 Tensor Core acceleration to high-order finite element operators in MFEM. They combine tensor-core programming with kernel fusion and data-layout optimizations, and evaluate the resulting kernels within a production wave-propagation application at large GPU scale. These studies demonstrate the importance of matching small tensor contractions to the hardware and reusing intermediate data. Our implementation also targets these tensor contractions, while addressing the representation changes introduced by Ozaki-II emulation. We retain intermediate tensors as integer residues across mode products and use rescaling and base extension to manage their representation range. The focus is therefore on reducing repeated conversion while preserving the factored contraction structure. 6

A PREPRINT - O CTOBER 6, 2026

3.2

Precision Emulation on Matrix Engines

The Ozaki scheme builds on error-free transformations of matrix multiplication developed by Ozaki et al. [4]. Subsequent work adapted this approach to the low-precision matrix units available on GPUs. Mukunoki et al. [5] investigated DGEMM using floating-point Tensor Cores, including accurate and reproducible variants, while Ootomo et al. [6] developed DGEMM emulation on integer matrix multiplication units. Together, these studies established how exact low-precision component products can support higher-precision matrix arithmetic. Ozaki scheme II introduced integer modular arithmetic into this GEMM-oriented emulation framework [7]. GEMMul8 implements this approach through a floating-point matrix interface and allows preprocessing to be reused for unchanged operands [8]. Such caching is useful for constant finite element operators. However, composing these GEMM calls still requires the changing tensor to pass through the floating-point interface between mode products, motivating our investigation of persistent residue representations. Recent work has improved both the efficiency and the capabilities of Ozaki-based emulation. To reduce intermediate data movement, EmuGEMM fuses operations within an emulated GEMM [10]. Its Scheme-I implementation combines component products and reconstruction in a persistent kernel. For Scheme II, it places modular reduction in the GEMM epilogue, avoiding a separate global-memory write of the INT32 accumulators before producing output residues. It also provides a fused complex kernel. These optimizations can support efficient residue products within our method, while our rescaling and base extension manage the representation passed between successive mode products. Another direction expands the precision formats available for emulation. Mukunoki [19] investigated FP8 Tensor Cores for DGEMM, and Uchino et al. [20] developed an FP8 quantization approach for Ozaki-II. Hayashi et al. [21] extended both Ozaki schemes to FP4 Tensor Cores by representing integers with base-13 FP4 limbs and keeping the relevant partial products exact under FP32 accumulation. Alongside these hardware-oriented developments, Schwarz et al. [22] addressed accuracy selection through Automatic Dynamic Precision. Their method selects the decomposition size using an exponent-span estimate and incorporates exception handling and fallback to native FP64, improving the reliability of emulated GEMM across different inputs. Ozaki-II has also been extended beyond real matrix multiplication. Uchino et al. [23] developed a Karatsuba-based formulation for complex matrix multiplication that combines real partial products in the residue domain before CRT reconstruction. This provides an existing example of combining several matrix products before reconstructing one logical output. Our setting extends the lifetime of the residue representation across successive tensor mode products, where each intermediate becomes an operand of the next contraction. A related application is the three-dimensional FFT studied by Matsuoka [24]. The proposed design uses Bailey decompositions and combines emulated DFT GEMMs, fused twiddle factors, and CRT reconstruction using Tensor Cores and fixed-point accumulation. Its analysis identifies conversion and integer-epilogue costs as important performance constraints. In the described algorithm, the intermediate between the two Bailey GEMM stages is reconstructed and then converted to residues again. This shares our motivation to reduce the cost of precision emulation in structured multidimensional computations. We instead investigate carrying the residue representation across tensor mode products through rescaling and base extension. The FFT study presents projected performance bounds, whereas our evaluation reports hardware measurements. Our method builds on the scaling rules, residue products, and CRT reconstruction of Ozaki scheme II. The original paper already considers products of more than two matrices and derives a sufficient CRT range condition for the combined integer result [7]. Our contribution concerns maintaining a reusable residue representation throughout a chain of tensor mode products. We combine established RNS rescaling and base-extension operations with operator-dependent range planning and separable scale handling to control intermediate integer growth while preserving the factored tensor computation. 3.3

Residue Arithmetic in Exact Computation and Homomorphic Encryption

Residue representations and mixed-radix conversion have a long history. Garner investigated residue number systems as a means of accelerating integer arithmetic through independent addition and multiplication on smaller residues, without carry propagation between components [15]. He also described a method for converting a residue tuple into a mixed-radix representation. This conversion, commonly known as Garner’s method, is identified in Ozaki scheme II as an alternative to direct CRT reconstruction [7]. Modular representations have also been used in other computational settings. FFLAS-FFPACK implements dense linear algebra over word-size prime fields using optimized numerical kernels [25]. Residue-based accelerators such as Mirage use these representations to address precision limitations in photonic DNN computation [26]. These studies illustrate how modular arithmetic can adapt numerical computations to the capabilities of the underlying hardware. 7

A PREPRINT - O CTOBER 6, 2026

Direct Direct multiplication by 2

v = 12

2v = 24

Ozaki scheme II Encode Basis

Multiply by 2

CRT

and reduce

reconstruction

Basis

(7, 5, 3) Centered residues

(−2, 2, 0)

v = 12

CRT weights eᵢ (P = 105)

(7, 5, 3) ×2 gives (−4, 4, 0) Centered residues

(15, 21, 70) all three moduli kept: |24| < P/2 = 52.5

Σ rᵢ eᵢ (mod P = 105)

45 − 21 + 0 = 24

(3, −1, 0) v = 12

α=1

2v = 24

α=1

2v = 24

2v = 24

α=1

Rescale + base extension (ours) Encode Basis

and reduce Basis

(7, 5, 3) Centered residues

(−2, 2, 0)

v = 12

Multiply by 2

(7, 5, 3) ×2 gives (−4, 4, 0) Centered residues

Drop modulus 3

Base extension

Basis

Basis

(7, 5) Centered residues

(1, −2)

Reconstruct, apply scale Σ rᵢ eᵢ (mod 105)

(7, 5, 3)

15 − 42 − 70 ≡ 8

Centered residues

Rescaled integer × scale

(1, −2, −1)

3 · 8 = 24

v′ = 8

v′ = 8

(3, −1, 0) v = 12

v′ is the rescaled integer,

α=1

2v = 24

α=1

v′ = 8

α′ = 3

(24 − [24]₃) / 3 = 8

α′ = 3

[1 + 7 · 1]₃ = −1

2v = 24

α′ = 3

same CRT weights (15, 21, 70)

α′ = 3 its reconstructing scale.

Figure 2: Three execution paths for multiplying v = 12 by two. The upper path computes 2v = 24 directly. The middle path illustrates Ozaki scheme II, which encodes the input, multiplies its residues, and reconstructs 24 using CRT. The lower path additionally applies rescaling and base extension before final reconstruction. Each residue tuple follows the displayed basis order. Dropping modulus 3 changes the represented integer from 24 to 8 and the reconstruction scale from α = 1 to α′ = 3. Base extension preserves this integer and scale. Final reconstruction and scaling recover 24 exactly in this example.

Managing the scale and basis of residue representations is particularly important in homomorphic encryption. The CKKS scheme uses rescaling to control scale growth during approximate arithmetic [11], and its full RNS variant develops modulus-switching procedures that operate on residue components [12]. The RNS variants developed by Bajard et al. [27] and Halevi et al. [13] for the Fan–Vercauteren (FV) scheme, also known as Brakerski–Fan–Vercauteren (BFV), provide further techniques for scaling and basis conversion. These works provide the arithmetic foundations for our range management. We apply rounded modulus division and mixed-radix base extension to the intermediate integers of a tensor mode-product chain. Our contribution lies in organizing these operations around the chain’s magnitude bounds and scaling requirements, allowing the residue representation to persist between contractions.

4

Representation Range Management through Rescaling and Base Extension

4.1

A Scalar Example

To illustrate how rescaling and base extension work together, Figure 2 compares three paths: direct multiplication, Ozaki scheme II, and residue-domain execution with rescaling and base extension. All three paths start from v = 12 and produce 2v = 24. The residue paths also track a reconstruction scale α. Multiplying the integer recovered from the residues by this scale gives the represented numerical value. The initial scale normally accounts for the scaling applied during input encoding. 8

A PREPRINT - O CTOBER 6, 2026

Here, the input is already an integer, so we begin with α = 1. When rescaling divides the represented integer by a modulus q, the reconstruction scale is multiplied by q to compensate for that division, apart from any rounding error. Both residue paths first encode v = 12 on the ordered basis (7, 5, 3), whose product is 105. Its centered residues are (−2, 2, 0). Multiplying each residue by two gives (−4, 4, 0), which becomes (3, −1, 0) after centered reduction. This tuple represents the doubled integer 24, with the reconstruction scale still equal to one. In the Ozaki scheme II path, all three moduli are retained, and CRT reconstruction recovers the integer directly. For the ordered basis (7, 5, 3), the CRT weights are (15, 21, 70). Each weight is congruent to one modulo its associated modulus and zero modulo the other two. The reconstructed result is therefore [3 · 15 + (−1) · 21 + 0 · 70]105 = [45 − 21]105 = 24. Since 2|24| < 105, centered reconstruction recovers the intended integer. The reconstruction scale remains α = 1, so the numerical result is also 24. The rescaling and base-extension path instead reduces the represented integer before final reconstruction. It first drops the modulus q = 3. Since the corresponding residue is [24]3 = 0, the rescaled integer is 2v − [2v]3 24 − 0 = = 8. (18) 3 3 This expression describes the change in the represented integer. In residue arithmetic, the remaining residues are computed directly using (15), without reconstructing 24. The result is (1, −2) on the reduced basis P ′ = {7, 5}, whose product is P ′ = 35. At the same time, the reconstruction scale changes from α = 1 to α′ = αq = 3. The smaller integer and the larger scale therefore still represent the numerical value 3 · 8 = 24. v′ =

The reduced basis is sufficient to identify v ′ = 8 because 2|v ′ | = 16 < P ′ . Base extension can therefore recover its missing residue modulo 3. Using the centered mixed-radix expansion 8 = 1 + 7 · 1 gives [1 + 7 · 1]3 = −1. Appending this residue produces (1, −2, −1) on the original basis (7, 5, 3). The represented integer remains 8, and the reconstruction scale remains 3. As shown in the figure, the same CRT weights now recover [1 · 15 + (−2) · 21 + (−1) · 70]105 = [−97]105 = 8. Applying the reconstruction scale then gives α′ v ′ = 3 · 8 = 24 = 2v. The recovery is exact in this example because 24 is divisible by the dropped modulus 3, so rescaling introduces no rounding error. This is not true for every input. Starting instead from v = 13 gives a doubled integer of 26, whose centered residue modulo 3 is −1. Rescaling therefore produces (26 − (−1))/3 = 9. With the updated scale of 3, the decoded result is 27 rather than 26. Base extension preserves the rescaled integer and cannot recover the precision lost in this rounding step. 4.2

Integer Growth and Reconstruction Range

In the scalar example, the represented integers remain within the centered reconstruction range at every step. For a chain of mode products, this condition must be maintained as the integer magnitudes change. The representation range is an essential constraint in Ozaki scheme II. For a residue basis with product P , CRT determines an integer only modulo P . Centered reconstruction selects a representative from an interval around zero, with a magnitude limit of approximately P/2. If the intended integer falls outside this interval, reconstruction returns a different integer with the same residues. The stored residues can therefore remain valid even when the intended value can no longer be recovered through centered reconstruction. This constraint must hold throughout a persistent mode-product chain. The input is encoded once, and each subsequent product updates the represented integers without reconstructing an intermediate floating-point tensor. Multiplication by an encoded operator and summation along a tensor mode can increase their magnitudes, so a basis that accommodates the first product may be insufficient for later products. Reducing each stored residue modulo its own modulus does not prevent this growth. Figure 3 illustrates how rescaling and base extension manage the relationship between integer magnitude and representation range. Let B be an upper bound on the absolute value of every represented integer entry. Let P denote the Q working residue basis and L ⊆ P the currently available basis, whose product is PL = p∈L p. A sufficient condition for recovering the intended integers is 2B < PL . (19) 9

A PREPRINT - O CTOBER 6, 2026

Magnitude bound log₂ B

Centered range log₂(PL/ 2) bits

0

20

40

60

80

100

120 116.8

Encode 53.0

15 moduli

116.8

Mode product 108.6

15 moduli 62.6

Rescale 54.4

8 moduli

116.8

Base extension 54.4

15 moduli

116.8

Next mode product 109.8

15 moduli Rescale reduces the integer. Extension restores the range.

Figure 3: Changes in the integer magnitude bound and the available centered reconstruction range. Filled bars show log2 B, and capped lines show log2 (PL /2) for the currently available basis. Rescaling removes seven moduli and reduces both quantities. Base extension restores the original range without changing the rescaled integer. The displayed exponents are rounded illustrative values, not measurements.

The figure shows both quantities on a logarithmic scale. The filled bars represent the magnitude bound log2 B, while the capped lines indicate the centered range limit log2 (PL /2). The working basis contains 15 moduli from the INT8 Ozaki-II table, which begins with 256, 255, 253, . . .. Their product satisfies log2 P ≈ 117.8, giving an initial centered range limit of approximately 116.8 on this scale. The initial magnitude bound is determined separately by the significand bits b, which bounds the encoded entries by 2b . With b = 53, the figure therefore starts from log2 B = 53. A mode product can increase the magnitude bound while leaving the representation range unchanged. Rescaling reduces this magnitude by dividing the represented integer by a selected modulus q, with rounding. It also removes q from the available basis, reducing the basis product from PL to PL /q. The figure illustrates seven successive modulus drops, which reduce both the magnitude bound and the available range. Base extension then computes residues of the rescaled integer for the removed moduli, restoring the working basis P. The integer remains unchanged, but the restored basis provides additional range for the next product. Together, rescaling and base extension reduce the integer magnitude while restoring representation capacity between successive mode products. Base extension does not reverse the rounding introduced by rescaling. 4.3

Rescaling and Drop planning

The procedure RnsDrop in Algorithm 1 applies the rescaling update from Section 2.3.1 to every tensor entry. For a selected modulus q, each remaining residue plane is updated using the corresponding entry from the q plane and a precomputed modular inverse. If an entry represents an integer v before the drop, the update represents v ′ = (v − rq )/q on L \ {q}, where rq = [v]q . The reconstruction scale changes from α to α′ = αq, giving the local error α′ v ′ − αv = −αrq ,

|α′ v ′ − αv| ≤ |α|q/2.

(20)

To determine the modulus-dropping schedule for the rescaling phase before execution, we first estimate how much each mode product can increase the integer magnitude. Let Āi denote the effective integer operator for the i-th mode product, including the scaling incorporated into that operator. Each output entry of the mode product is an inner product between one row of Āi and a tensor fiber. Let j index the output entries along that fiber and k index the entries being contracted. The maximum absolute row sum of the operator is ρi = max

n X

1≤j≤n

k=1

10

|(Āi )jk |.

(21)

A PREPRINT - O CTOBER 6, 2026

Algorithm 1 RnsDrop for rescaling by one modulus Require: Centered residue planes R representing integer entries v on L, drop modulus q ∈ L, scale α Require: The planned quotient satisfies 2B ′ < PL /q Ensure: Residues of v ′ = (v − [v]q )/q on L \ {q}, with decoding scale α′ = αq 1: procedure R NS D ROP(R, L, q, α) 2: S ← Rq 3: for all p ∈ L \ {q} do 4: Rp ← [(Rp − S)(q −1 mod p)]p 5: end for 6: Remove plane Rq and modulus q from L 7: Record q in the scale metadata for α′ = αq 8: return (R, L, α′ ) 9: end procedure

Algorithm 2 DropPlan from operator row-sum bounds, independent of the input tensor Require: Encoding parameter b, row-sum bounds (ρi )di=1 , basis P, restoration flag τ ∈ {0, 1} Ensure: Drop counts (gi )di=1 or failure ⊥ 1: procedure D ROP P LAN(b, (ρi ), P, τ ) 2: B ← 2b 3: L ← PQ 4: PL ← p∈L p 5: for i = 1, . . . , d do 6: B ← Bρi 7: if 2B ≥ PL then 8: return ⊥ 9: end if 10: gi ← 0 11: if i < d then 12: while |L| > 1 and ⌊B/ min L⌋ ≥ 2b do 13: q ← min L 14: B ′ ← ⌊B/q⌋ + 1 15: if 2B ′ ≥ PL /q then 16: return ⊥ 17: end if 18: (L, PL , B) ← (L \ {q}, PL /q, B ′ ) 19: gi ← gi + 1 20: end while 21: if τ = 1 then 22: L ← PQ 23: PL ← p∈P p 24: end if 25: end if 26: end for 27: return (g1 , . . . , gd ) 28: end procedure

For an integer fiber x whose entries satisfy |xk | ≤ B, the triangle inequality gives |(Āi x)j | ≤

n X

|(Āi )jk | |xk | ≤ B

k=1

n X

|(Āi )jk | ≤ Bρi .

(22)

k=1

Since the same operator acts on every fiber along the contracted mode, Bρi bounds every entry of the output tensor. The planner can therefore update the magnitude bound by B ← Bρi at each mode product. This is the standard induced infinity-norm bound, with ρi = ∥Āi ∥∞ . Using these operator bounds, DropPlan in Algorithm 2 determines how many moduli to drop after each mode product. It starts from B = 2b , where b is the significand bits, so that every encoded input integer has magnitude at most 2b , and updates B ← Bρi for each mode product. The updated bound must satisfy 2B < PL on the currently available basis. 11

A PREPRINT - O CTOBER 6, 2026

Algorithm 3 Planner for precision selection on a fixed working basis Require: Operator row-sum bound functions (ρi ), fixed working basis P with s = |P| Ensure: Selected encoding parameter b and drop counts g, or failure ⊥ 1: procedure P LANNER((ρi ), P) 2: b ← 53 3: g ← D ROP P LAN(b, (ρi (b)), P, 1) 4: while g = ⊥ do 5: b←b−1 6: if b < 1 then 7: return ⊥ 8: end if 9: g ← D ROP P LAN(b, (ρi (b)), P, 1) 10: end while 11: return (b, g) 12: end procedure

After each nonfinal product, the planner considers dropping the smallest available modulus q. It schedules a drop when ⌊B/q⌋ ≥ 2b and uses the conservative bound B ′ = ⌊B/q⌋ + 1

(23)

for the rescaled integer. The additional unit accounts conservatively for rounding during division. Before accepting the drop, the planner checks 2B ′ < PL /q so that the quotient remains recoverable on the reduced basis. It then updates the bound and the available basis and considers the next smallest modulus. This process continues while the threshold is satisfied and more than one modulus remains. A failed range check causes the candidate schedule to be rejected. This selection is greedy. At each step, the planner chooses the smallest available modulus rather than comparing alternative drop sequences. For our proposed method, the restoration flag is set to τ = 1. The planner therefore assumes that base extension restores the working basis after the scheduled drops and uses that basis when checking the next product. The flag allows the same planner to support the rescale-only variant in Appendix A.3, which uses τ = 0. When applying only the proposed method, τ remains fixed at one. It returns a drop count gi for each mode product. During execution, the same sequence is reproduced by calling RnsDrop gi times, each time removing the smallest available modulus, before applying base extension. The final product is decoded directly, so gd = 0. The remaining decision is the significand bits b. The procedure Planner in Algorithm 3 selects it on the supplied working basis. Because changing b changes the encoded integer operators, we write ρi (b) for the corresponding maximum absolute row sum for the i-th mode product. Starting from b = 53, which is the FP64 significand, the procedure prepares the effective integer operators, computes their bounds, and calls DropPlan with restoration enabled. If the schedule fails a range check, it decreases b by one and repeats these steps. It returns the first successful significand bits and its drop counts. If no positive b passes the checks, it reports failure. The working basis remains fixed throughout this search. Planning depends on the operators, their contraction order, and the working basis. It does not depend on the individual input tensor because the encoded input entries share the bound 2b . The plan and prepared operator residues can therefore be reused across tensors with the same operator chain. The selected b specifies the encoding precision accepted by this planner, rather than guaranteeing b accurate bits in the final result. In particular, rescaling can lose accuracy when individual entries are much smaller than the propagated bound, including when cancellation occurs. 4.4

Base Extension Using Shared Mixed-Radix Digits

After the scheduled drops, base extension restores the removed moduli while preserving the rescaled Qm integer. Let P ′ = {q1 , . . . , qm } be the reduced basis, processed in the order q1 , . . . , qm , with product P ′ = i=1 qi . For each rescaled integer entry v ′ , we require 2|v ′ | < P ′ so that its remaining residues uniquely determine its value. The currently available basis is L = P ′ , and the moduli to be restored form the target set D = P \ P ′ . To compute the missing residues, we first convert the remaining residues into mixed-radix digits. Algorithm 4 provides the procedures MrcDigits and BaseExtend for these two operations. For each tensor entry, BaseExtend calls MrcDigits once and reuses the resulting digits for every modulus in D. The conversion follows the successive modular elimination underlying Garner’s method [15]. Unlike the nonnegative digits introduced in Section 2.3, we use centered digits. If the basis contains an even modulus, we place it first, followed by the odd moduli. The descending basis used here satisfies this ordering when it contains modulus 256. 12

A PREPRINT - O CTOBER 6, 2026

Algorithm 4 MrcDigits and BaseExtend for centered residue conversion Require: Centered residues R of v ′Qon reduced basis P ′ = {q1 , . . . , qm }, processed in that order Require: 2|v ′ | < P ′ , where P ′ = i qi , with the even modulus first if present Require: Removed moduli D = P \ P ′ Q Require: Precomputed inverses and weights ωi,p = [ j<i qj ]p Ensure: Base extension represents the same integer v ′ on P ′ ∪ D = P 1: procedure M RC D IGITS(R, P ′ ) 2: U ←R ▷ working copy of the input residues 3: for i = 1, . . . , m do 4: di ← U i 5: for j = i + 1, . . . , m do 6: Uj ← [(Uj − di )(qi−1 mod qj )]qj 7: end for 8: end for 9: return (d1 , . . . , dm ) 10: end procedure 11: procedure BASE E XTEND(R, P ′ , D) 12: (d1 , . . . , dm ) ← M RC D IGITS(R, P ′ ) 13: for all p ∈ D do 14: t←0 15: for i = 1, . . . , m do 16: t ← [t + di ωi,p ]p 17: end for 18: Rp ← t 19: end for 20: return R 21: end procedure

In this representation, each digit di has a positional weight Wi , analogous to the powers of ten in decimal notation. For mixed-radix digits, the weight is the product of all preceding moduli, Y W1 = 1, Wi = qj . (24) j<i

The rescaled integer is then expressed as v′ =

m X

di W i ,

−⌊qi /2⌋ ≤ di ≤ ⌈qi /2⌉ − 1.

(25)

i=1

Successive modular elimination produces digits whose weighted sum is congruent to v ′ modulo every modulus in P ′ . With the stated ordering, the possible weighted sums cover exactly the centered interval [−⌊P ′ /2⌋, ⌈P ′ /2⌉ − 1]. Since the intended integer also lies in this interval, uniqueness guarantees that the expansion equals v ′ . The modulus ordering is therefore required for this centered representation. Once the digits are available, each missing residue can be computed by reducing the same expansion modulo its target modulus. For every removed modulus p ∈ D, we precompute the reduced positional weights ωi,p = [Wi ]p and evaluate "m # X Rp = di ωi,p . (26) i=1

p

Thus, the digit conversion is performed only once per tensor entry, rather than repeated for each restored modulus. It requires m(m − 1)/2 modular elimination updates, followed by m weighted terms for each target modulus. Only the reduced weights are needed during extension, so the computation does not materialize the full integer v ′ or reconstruct an intermediate FP64 tensor. Appending the new residues restores L = P and its product P . The existing residues, the reconstruction scale α′ , and the rescaled magnitude bound B ′ remain unchanged. The restored residue planes are then available for the next mode product, subject to the range condition checked by the planner. 13

A PREPRINT - O CTOBER 6, 2026

Z = T ×ₘᵢ Aᵢ Fiber scales

T FP64

Encode

Residue GEMM

Reconstruct

Fiber scaling

(n × n)(n × n²)

MRC + Horner

To residues

Then reduce mod p

Apply scales

Rₚ

Z FP64

Product residues

Cached operator Residues Āᵢ,ₚ and row scales Row scales

Next mode starts from the FP64 result

Figure 4: One tensor mode product in OzII-Composed. The input tensor is encoded, multiplied independently on each residue plane, and decoded to FP64. Operator residues can be cached when the operator and encoding settings remain unchanged. Each modulus requires an (n × n)(n × n2 ) GEMM. The next mode product repeats the conversion for the FP64 intermediate.

5

Applying Ozaki Scheme II to Tensor Mode Products

This section combines the residue operations with Ozaki scheme II to execute chained tensor mode products. It first describes the composition of separate Ozaki scheme II mode products, then explains how rescaling and base extension allow intermediate tensors to remain in residue form across the chain. Finally, it describes the GPU implementation of these algorithms. 5.1

Composed Ozaki-II Execution

Figure 4 illustrates one mode product in OzII-Composed. For the i-th mode product, unfolding along mode mi arranges the n2 fibers as columns of an n × n2 matrix. Each fiber receives a power-of-two encoding scale, and each operator row receives its own scale. After scaling and rounding to integers, both integer operands are converted to centered residue planes. For each modulus, the operator residue matrix multiplies the corresponding tensor unfolding, and the accumulated result is reduced to centered residues. Reconstruction and inverse scaling produce the FP64 tensor used by the next mode product. This follows the GEMM-based organization of Ozaki scheme II [7]. Algorithm 5 defines the conversion procedures. For a group of entries G that share an encoding scale, the exponent eG is chosen so that the largest scaled magnitude is less than 2b before rounding. An all-zero group is encoded as zero with exponent zero. Unlike the truncation used in the original Ozaki scheme II formulation, our implementation rounds scaled operands to the nearest integer, with ties to even. The resulting integer magnitudes are bounded by 2b . For the tensor input, an encoding group is one fiber in OzII-Composed and the complete tensor in the persistent method. The procedure ToRns then converts an integer tensor T̄ into one centered residue plane Rp = [T̄ ]p for each modulus p in the supplied basis L, with reduction applied entrywise. For the i-th mode product, let ei,c be the encoding exponent of fiber column c and fi,a the encoding exponent of operator row a. If T denotes the input unfolding and k indexes its contracted dimension, the integer operands are (T̄ )kc = round(2ei,c Tkc ),

(Āi )ak = round(2fi,a (Ai )ak ).

(27)

Every term contributing to output entry (a, c) therefore has the same combined scale. After reconstruction, that entry is multiplied by σi,ac = 2−fi,a −ei,c to undo the encoding scales. Together, these factors form an entrywise reconstruction scale field determined by the operator rows and tensor fibers. Algorithm 6 repeats these steps over the chain. The supplied Q basis P contains the first s entries in GEMMul8’s table of pairwise coprime moduli, starting with 256. With P = p∈P p, the sufficient condition n22b < P/2 ensures that each length-n integer dot product remains within the centered reconstruction range. The algorithm selects the largest positive integer b ≤ 53 satisfying this condition and uses it to encode both operands. The basis remains unchanged throughout execution. 14

A PREPRINT - O CTOBER 6, 2026

Algorithm 5 Encoding and reconstruction primitives 1: procedure E NCODE(T, b, G) 2: for all groups G ∈ G do 3: if max |G| = 0 then 4: (Ḡ, eG ) ← (0, 0) 5: else 6: Choose EG with 2EG −1 ≤ max |G| < 2EG 7: eG ← b − EG 8: Ḡ ← round(2eG G) 9: end if 10: end for 11: return (T̄ , e) 12: end procedure 13: procedure T O R NS(T̄ , L) 14: return {Rp = [T̄ ]p }p∈L 15: end procedure 16: procedure R ECONSTRUCT FP64(R, L, σ) 17: Order L = {q1 , . . . , qm } with the even modulus first 18: (d1 , . . . , dm ) ← M RC D IGITS(R, L) 19: Vb ← fl64 (dm ) 20: for j = m − 1, . . . , 1 do 21: Vb ← fl64 (qj Vb + dj ) 22: end for 23: return fl64 (σ ⊙ Vb ) 24: end procedure

For fixed operators and encoding settings, the operator residues are prepared once and reused. The intermediate tensor, however, is encoded and decoded for every mode product. A three-direction pass therefore requires three tensor encodings and three reconstructions. To perform the reconstruction, ReconstructFP64 uses centered mixed-radix conversion followed by FP64 Horner evaluation and scaling. This reconstructs the integer represented by the residues without evaluating the direct CRT formula used in Ozaki scheme II. Mixed-radix conversion is also used for base extension in our persistent method, so this choice allows both operations to share the same conversion procedure. The composed variant uses the same reconstruction procedure for consistency. The digit conversion is exact, while FP64 Horner evaluation and subsequent scaling can introduce rounding errors. The conversion in Algorithm 4 operates on the ordered basis supplied to the procedure, which is the full working basis for OzII-Composed. In the pseudocode, fl64 denotes floating-point evaluation in FP64, and ⊙ denotes entrywise multiplication by the supplied reconstruction scale field σ. The range condition guarantees a unique integer for reconstruction, although its FP64 evaluation and subsequent scaling may introduce rounding. The algorithms assume finite operands and representable scales and intermediate results during floating-point evaluation. 5.2

Persistent Ozaki-II Execution with Rescaling and Base Extension

Algorithm 7 combines tensor mode products with the PlanRescaleBE, RnsDrop, and BaseExtend procedures. Given the operators, their contraction order, and the working basis P, the planner selects the significand bits b and the number of drops gi after each mode product. The effective integer operators Āi are then prepared using the same scaling and rounding rules used to compute the bounds ρi (b). Operators already prepared for the accepted b can be reused. The input tensor is encoded once using a single exponent e. Its initial reconstruction scale is therefore the scalar α = 2−e . The operator row scales are tracked separately and combined into a final reconstruction scale field, denoted by Γ. Unlike α, this field can vary between tensor entries. Their treatment is described below. For the i-th mode product, the encoded operator multiplies each tensor residue plane along the contracted mode. The products use centered INT8 operands with INT32 accumulation, followed by centered modular reduction. The accumulations must remain within the INT32 range, in addition to satisfying the separate CRT range conditions. When gi > 0, the algorithm starts from the working basis L = P and calls RnsDrop gi times. Each call removes the smallest available modulus, following the planner’s order, and updates the scalar reconstruction scale to account for the division. The returned scale α′ becomes α for the next operation. After the scheduled drops, the remaining basis is P ′ = L. A 15

A PREPRINT - O CTOBER 6, 2026

Algorithm 6 OzII-Composed for successive tensor mode products Require: Tensor X, operators (Ai )di=1 , modes (mi )di=1 Require: Fixed residue basis P with s pairwise coprime moduli 1: procedureQO Z IIC OMPOSED(X, P, (Ai ), (mi )) 2: P ← p∈P p 3: Choose the largest b ∈ {1, . . . , 53} with n22b < P/2 4: if no such b exists then 5: return ⊥ 6: end if 7: for i = 1, . . . , d do 8: (Āi , fi ) ← E NCODE(Ai , b, Rows(Ai )) 9: Cache Āi,p = [Āi ]p for p ∈ P 10: end for 11: T ←X 12: for i = 1, . . . , d do 13: (T̄ , ei ) ← E NCODE(T, b, Fibersmi (T )) 14: R ← T O R NS(T̄ , P) 15: for all p ∈ P do 16: Rp ← [Rp ×mi Āi,p ]p 17: end for 18: σi ← 2−ei 2−fi with output-index broadcasting 19: T ← R ECONSTRUCT FP64(R, P, σi ) 20: end for 21: return T 22: end procedure

single call to BaseExtend restores the removed moduli in D = P \ P ′ , reusing the mixed-radix digits across all targets. The next mode product therefore uses the full working basis P again. Since gd = 0, the final mode product proceeds directly to reconstruction. Figure 5 compares this execution with OzII-Composed. Both paths apply the same sequence of tensor mode products, but handle the intermediate tensor differently. OzII-Composed reconstructs an FP64 tensor and encodes it again before the next product. OzII-RescaleBE instead reduces the integer magnitude through rescaling and restores the representation range through base extension. The intermediate tensor remains in residue form, accompanied by its reconstruction scales, until the final reconstruction. For the final reconstruction to recover the appropriate numerical scale, the operator row scales must also be tracked throughout the chain. In a three-direction pass, each operator acts along a different mode. Its row scales depend only on the corresponding output index, so they can be applied after all three products. Using the row encoding exponents fi,a introduced for OzII-Composed, the final operator reconstruction scale field is Γa1 a2 a3 = 2−f1,a1 2−f2,a2 2−f3,a3 .

(28)

The complete reconstruction scale is αΓ. Here, α accounts for the initial tensor encoding and all modulus divisions, while Γ accounts for the operator row scales. A longer chain requires an additional step when a mode is contracted again. The pending row scales from the previous product on that mode now belong to the entries being summed. They must therefore be incorporated into the next operator before its integer encoding. This is done by multiplying each column of that operator by the pending reconstruction scale for the corresponding contracted entry. Encoding this effective operator produces new row scales, which replace the pending scales on that mode. The pending scales on the other modes remain unchanged. At the end of the chain, Γ is formed from the remaining scale on each mode. An earlier scale that has already been incorporated into a later operator is therefore not applied again during final reconstruction. The planner computes ρi (b) from these effective integer operators so that its bounds describe the products actually executed. For each candidate b, operator preparation starts with all pending scales equal to one and processes the complete operator sequence again. Once b is selected, the resulting operator residues and final scale field can be cached and reused for tensors with the same operator chain. 16

A PREPRINT - O CTOBER 6, 2026

Algorithm 7 OzII-RescaleBE for persistent tensor mode products Require: Tensor X, operators (Ai )di=1 , modes (mi )di=1 Require: Fixed residue basis P with s pairwise coprime moduli 1: procedure O Z IIR ESCALE BE(X, P, (Ai ), (mi )) 2: Define (ρi (b)) from effective operator encoding at precision b 3: π ← P LAN R ESCALE BE((ρi ), P) 4: if π = ⊥ then 5: return ⊥ 6: end if 7: (b, (gi )) ← π 8: Prepare effective operators (Āi ) and scale field Γ at b 9: Cache Āi,p = [Āi ]p for p ∈ P 10: (X̄, e) ← E NCODE(X, b, {X}) 11: α ← 2−e 12: R ← T O R NS(X̄, P) 13: for i = 1, . . . , d do 14: for all p ∈ P do 15: Rp ← [Rp ×mi Āi,p ]p 16: end for 17: if gi > 0 then 18: L←P 19: for j = 1, . . . , gi do 20: q ← min L 21: (R, L, α′ ) ← R NS D ROP(R, L, q, α) 22: α ← α′ 23: end for 24: P′ ← L 25: D ← P \ P′ 26: R ← BASE E XTEND(R, P ′ , D) 27: L←P 28: end if 29: end for 30: return R ECONSTRUCT FP64(R, P, αΓ) 31: end procedure

5.3

Implementation

The implementations serve as a proof of concept and use Python, selected PyTorch utilities, and CuTe DSL version 4.5.2 for ease of development. Python handles operator preparation, range-bound calculations, drop planning, and the construction of modular inverse and base-extension tables. PyTorch provides timing and memory-measurement utilities. Its einsum operation is also used to implement the FP64 tensor mode-product baseline, using cuBLAS internally for the tested contractions. The GPU kernels are written in CuTe DSL and use the selected working basis, number of significand bits, and drop schedule. Plans and prepared operator data are reused across tensors with the same operators and configuration. The target GPU is the NVIDIA RTX PRO 6000 Blackwell, whose Tensor Cores support INT8 matrix multiplyaccumulate (MMA) instructions. However, CuTe DSL 4.5.2 does not expose a warp-level INT8 MMA operation for the sm_120 architecture [28]. A user-defined operation using the DSL’s @dsl_user_op decorator therefore issues the required MMA instruction through LLVM inline assembly [29]. The implementations use custom residue GEMM kernels to handle the different unfoldings required by the three tensor mode products. In particular, the middle-mode unfolding is not directly available as a single contiguous matrix in the tensor’s storage layout. The custom kernels implement the indexing for each contracted mode and combine INT8 matrix multiplication with centered modular reduction, reducing intermediate data traffic. Alongside OzII-Composed and OzII-RescaleBE, the implementation includes OzII-NoRescale and OzII-Rescale as comparison variants. Their algorithms and the rescale-only planner are provided in Appendix A. 17

A PREPRINT - O CTOBER 6, 2026

OzII-Composed

OzII-RescaleBE

3 encodings · 3 reconstructions

1 encoding · 1 reconstruction

X (FP64)

X (FP64)

Encode

Encode

Mode 1

Mode 1

15 residue GEMMs

15 residue GEMMs

Reconstruct

Rescale 15 → 8 moduli

FP64 intermediate

Base extension

Encode

8 → 15 moduli

Mode 2

Mode 2

15 residue GEMMs

15 residue GEMMs

Reconstruct

Rescale 15 → 8 moduli

FP64 intermediate

Base extension

Encode

8 → 15 moduli

Mode 3

Mode 3

15 residue GEMMs

15 residue GEMMs

Reconstruct and apply scales

Reconstruct and apply scales

Y (FP64)

Y (FP64)

FP64 between products

Residues between products Illustrative bases and schedule

Figure 5: Three successive tensor mode products. OzII-Composed reconstructs an FP64 intermediate and encodes it again between products. OzII-RescaleBE retains the intermediate residues and applies rescaling and base extension between products. Table 1: Hardware and software environment. Hardware GPU CUDA target Streaming multiprocessors Device memory Peak memory bandwidth Peak FP64 throughput Peak INT8 Tensor Core throughput

NVIDIA RTX PRO 6000 Blackwell Max-Q Workstation Edition sm_120 188 96 GB GDDR7 1792 GB/s ≈ 1.71 TFLOP/s 877.9 TOP/s

Software CuTe DSL PyTorch CUDA cuBLAS

6

Evaluation

6.1

Experimental Environment

4.5.2 2.11.0 13.3 13.6

Table 1 summarizes the hardware and software environment. The theoretical peak performance values are calculated from NVIDIA’s architecture specifications using the specified 2.28 GHz boost clock and counting each multiply–add as two operations [30]. Tensor Core rates refer to dense operations without structured sparsity. The FP64 rate is calculated assuming two FP64 units per SM. The experiments are characterized by the tensor mode size n, chain depth d, modulus count s, and number of retained significand bits b. Each input tensor has dimensions n × n × n, and each mode product applies an n × n operator along 18

A PREPRINT - O CTOBER 6, 2026

one direction. The chain depth d counts the successive mode products, so a pass over all three directions has d = 3. Longer chains revisit these directions while the tensor remains three-dimensional. Two baselines implement the tensor contraction using GEMMul8 [8]. GEMMul8-Composed applies the d mode products through separate GEMMul8 calls, following the same composition as OzII-Composed. GEMMul8-KRON uses the Kronecker-product representation introduced in Section 2.1. For a single tensor, it multiplies a dense n3 × n3 operator by the vectorized input of length n3 . For a batch sharing the same operator, the vectorized inputs are collected as columns and processed in one matrix multiplication. This baseline produces a larger dense GEMM even at small tensor mode sizes, providing a workload better suited to GEMMul8’s optimized matrix multiplication. The dense operator is prepared once and reused across evaluations. All residue methods start from the same basis containing s = 15 moduli, with one residue plane per modulus. In OzII-RescaleBE, rescaling temporarily reduces the available basis to P ′ , and base extension restores all s moduli before the next mode product. OzII-Rescale does not restore removed moduli. Its final modulus counts are therefore reported alongside the initial count. For the OzII implementations, b is selected according to each method’s range conditions on the supplied basis. It sets the bound 2b on the encoded input integer magnitudes and, where rescaling is used, the nominal threshold in the drop planner. For GEMMul8, the reported b is the nominal bit count determined by the range rule, while the actual operand scales are selected internally. The number of significand bits b is distinct from the modulus count s and does not by itself determine the final numerical accuracy. 6.2 6.2.1

Accuracy Evaluation Datasets

Three datasets with different input distributions are used for the accuracy evaluation. Uniform. Each tensor entry is sampled independently and uniformly between −1 and 1. This dataset provides an unstructured baseline without an imposed wide exponent distribution. Smooth. Each tensor samples the field X(x, y, z) = a sin(3x + ϕ1 ) cos(5y + ϕ2 ) exp(cz)

(29)

on an n × n × n unit grid. For each tensor, the amplitude is sampled uniformly between 0.5 and 2, the phases uniformly between 0 and 2π, and the growth rate uniformly between 0.5 and 1.5. This construction provides a simple model of smooth nodal fields in high-order finite element calculations. The analytic field is separable and has tensor rank one. In exact arithmetic, mode products preserve this structure, and each output entry factors into three one-dimensional contributions. Cancellation in any contribution can produce a small output entry and amplify its relative error. Rounding the sampled field to FP64 can perturb exact separability. Wide160. For each entry, we sample a mantissa uniformly between 1 and 2 and an integer exponent uniformly from −80 to 80, inclusive. The mantissa is multiplied by two raised to the sampled exponent. The result is then assigned a positive or negative sign at random. The name refers to the 160-bit span between the smallest and largest exponents. This dataset tests the effect of rounding during encoding when entries sharing a scale have widely different magnitudes. With a fixed b-bit encoding, sufficiently small entries relative to the group’s largest entry can be rounded to zero. This applies to the per-fiber scaling in the composed methods and the single scale per tensor in the persistent methods. A small normwise error can therefore coexist with larger relative errors in small output entries. For each dataset and chain depth d ∈ {3, 4, 5, 6, 7, 8}, 32 tensors with n = 8 are evaluated. Each of the d operator matrices has dimensions n × n, with entries sampled independently and uniformly between −1 and 1. The contracted modes cycle through 1, 2, 3, 1, 2, 3, . . .. Thus, a depth-eight chain applies three mode products along each of modes 1 and 2 and two along mode 3. All three datasets use the same operator distribution and mode sequence, differing only in how the input tensor X is constructed. 6.2.2

Error Metrics

For each dataset and depth, the index i below ranges over all 32n3 = 16384 output entries in the test batch. The reference Y ref is computed by evaluating the chain using exact rational arithmetic on the FP64 inputs and operators with Python’s Fraction class. It serves as the ground truth for the evaluation. Therefore, the FP64 results computed using PyTorch’s torch.einsum are also evaluated against this reference. 19

A PREPRINT - O CTOBER 6, 2026

Table 2: Selected encoding parameter b for the accuracy tests with n = 8 and initial modulus count s = 15 for every residue method. Method

d=3

d=4

d=5

d=6

d=7

d=8

— 53 27 27 53 53 53

— 53 21 21 53 53 53

— 53 17 17 53 53 53

— 53 14 14 53 53 53

— 53 12 12 53 53 53

— 53 10 10 53 53 53

FP64 einsum OzII-Composed OzII-NoRescale OzII-Rescale OzII-RescaleBE GEMMul8-Composed GEMMul8-KRON

The normwise relative error is εnorm =

maxi |Yi − Yiref | . maxi |Yiref |

We also compute an error for each output entry,  |Yi − Yiref |/|Yiref |, Yiref ̸= 0, εi = |Yi − Yiref |, Yiref = 0,

(30)

(31)

and report its maximum and median, εmax = max εi , i

6.2.3

εmed = mediani εi .

(32)

Modulus and Encoding Configurations

All residue methods start from the first s = 15 moduli in GEMMul8’s ordered modulus table. Table 2 reports the significand bits selected for this fixed basis at n = 8. The selected counts are the same for Uniform, Smooth, and Wide160. Ozaki scheme II discusses successive matrix multiplications and derives a range condition for reconstructing their combined integer result [7]. As more encoded operands are multiplied, the integer result can require a larger representation range. On a fixed residue basis, satisfying this condition can therefore require fewer encoding bits per operand. The same constraint applies to a chain of tensor mode products whose intermediate integers remain in residue form. OzII-Composed and GEMMul8-Composed reconstruct the result after every mode product, so their basis only needs to accommodate one product at a time. OzII-Composed uses b = 53 throughout these tests. The nominal significand bits is also 53 for GEMMul8-Composed and GEMMul8-KRON. For KRON, the selected basis accommodates a single dense product with inner dimension n3 = 512. GEMMul8 selects the actual operand scales internally. GEMMul8-Composed uses its standard configuration in this accuracy comparison. OzII-NoRescale instead retains the integer representation across the complete chain without rescaling. To accommodate its growing magnitude within the same 15-modulus basis, the selected significand bits decreases from 27 at d = 3 to 10 at d = 8. OzII-Rescale reduces the integer magnitude between products, but each drop also reduces the available representation range because the removed moduli are not restored. Its planner selects the same significand bits as OzII-NoRescale at the tested depths. Its final modulus counts are listed below the table. OzII-RescaleBE reduces the integer magnitude and then restores the working basis between products. Its planner accepts b = 53 at every tested depth. For d = 3 through 7, seven moduli are dropped and restored after each nonfinal mode product. At d = 8, seven moduli are dropped after each of the first six products, and eight after the seventh. Base extension restores all 15 moduli before the next product. These range checks permit the selected significand bits, but do not eliminate the rounding introduced by rescaling. 6.2.4

Accuracy Results

Figure 6 compares the normwise errors across the three datasets. The upper row shows OzII-NoRescale and OzIIRescale, whose significand bits decrease with depth. The lower row shows the other residue methods together with the FP64 einsum. Different logarithmic vertical ranges are used for the two rows to show their respective error levels. OzII-NoRescale and OzII-Rescale already have normwise errors of approximately 10−8 at d = 3, where their selected significand bit is only b = 27. As depth increases, the available range must accommodate a longer chain, and the 20

A PREPRINT - O CTOBER 6, 2026

FP64 einsum OzII-Composed

OzII-NoRescale OzII-Rescale

(a) Uniform

OzII-RescaleBE GEMMul8-Composed

(b) Smooth

GEMMul8-KRON

(c) Wide160

100 10−2

Normwise relative error

10−4 10−6 10−8

bits b = 27, 21, 17, 14, 12, 10

(d)

bits b = 27, 21, 17, 14, 12, 10

(e)

bits b = 27, 21, 17, 14, 12, 10

(f)

10−12

10−13

10−14

10−15

10−16 3

4

5

6

Chain depth d

7

8

3

4

5

6

Chain depth d

7

8

3

4

5

6

7

8

Chain depth d

Figure 6: Normwise relative error against the reference at n = 8. Columns correspond to Uniform, Smooth, and Wide160. Panels (a)–(c) show OzII-NoRescale and OzII-Rescale, with significand bits for d = 3, . . . , 8. Panels (d)–(f) show the FP64 einsum and the emulation-based methods.

selected count falls to b = 10 at d = 8. The errors of OzII-NoRescale consequently increase to approximately 10−3 on all three datasets. OzII-Rescale uses the same significand bits but also introduces rounding between products. Its depth-eight errors reach 0.423, 1.32, and 1.00 on Uniform, Smooth, and Wide160, respectively. Although dropping moduli reduces the integer magnitude, it also reduces the representation range. These results show that rescaling without base extension does not maintain high accuracy for the tested chains on the fixed 15-modulus basis. OzII-Composed maintains normwise errors below 8 × 10−16 across all tested depths and datasets. OzII-RescaleBE also achieves small errors through d = 7, while retaining the intermediate tensor in residue form. Its errors range from 2.51 × 10−16 to 5.66 × 10−16 on Uniform. They are larger on Smooth and Wide160, reaching 1.69 × 10−15 and 2.60 × 10−15 , respectively. Unlike OzII-Composed, OzII-RescaleBE rounds the represented integers during intermediate rescaling. Entries that are small relative to the rescaling divisor are more sensitive to this rounding. Cancellation in Smooth and the wide input magnitude range in Wide160 can make such effects more pronounced. The observed differences are consistent with this sensitivity, although the tests do not separately measure the contributions of encoding and rescaling errors. At d = 8, the normwise errors of OzII-RescaleBE rise to 6.74 × 10−14 , 3.37 × 10−13 , and 6.00 × 10−13 on Uniform, Smooth, and Wide160. This increase coincides with the additional modulus drop after the seventh mode product, which introduces another rounded division. Base extension restores the representation range but does not reverse this rounding. Retaining b = 53 at encoding therefore does not guarantee the same accuracy at every chain depth. GEMMul8-Composed remains below 8 × 10−16 in normwise error across all tests. GEMMul8-KRON generally has larger errors on Uniform and Smooth, although its errors remain below 2 × 10−15 across all three datasets. Both configurations have a nominal significand bit of b = 53 at n = 8. A possible explanation is the scaling used when the dense operator is cached through skip_scalA. Reusing the operator fixes its encoding scales, and the cached execution path may select more conservative scales for the tensor operand, increasing its encoding error. Table 3 reports maximum and median elementwise errors at depths 3 and 8. OzII-RescaleBE maintains substantially smaller median errors than OzII-NoRescale and OzII-Rescale on all three datasets. At d = 3, its median errors range 21

A PREPRINT - O CTOBER 6, 2026

Table 3: Maximum and median elementwise errors at depths 3 and 8, using 32 tensors per dataset with n = 8. All residue methods start with s = 15. Method

d

Uniform Max εmax

Smooth

Median εmed

Max εmax

Wide160

Median εmed

Max εmax

Median εmed

−12

FP64 einsum 3 8.26 × 10 OzII-Composed 3 3.95 × 10−12 OzII-NoRescale 3 1.55 × 10−5 OzII-Rescale 3 4.58 × 10−5 OzII-RescaleBE 3 9.49 × 10−13 GEMMul8-Composed 3 3.09 × 10−12 GEMMul8-KRON 3 2.37 × 10−12

−16

1.93 × 10 2.53 × 10−16 7.82 × 10−9 8.03 × 10−9 1.87 × 10−16 1.51 × 10−16 7.28 × 10−16

−11

9.35 × 10 8.47 × 10−10 1.44 × 10−1 1.72 × 10−1 4.52 × 10−9 1.24 × 10−10 7.54 × 10−9

−16

4.33 × 10 6.53 × 10−16 4.00 × 10−8 4.26 × 10−8 1.47 × 10−15 2.09 × 10−16 2.67 × 10−15

−12

1.09 × 10 9.18 × 10−13 1.44 × 10−3 1.99 × 10−3 3.90 × 10−11 9.17 × 10−13 4.55 × 10−12

1.61 × 10−16 3.05 × 10−16 1.72 × 10−8 2.54 × 10−8 1.54 × 10−15 1.64 × 10−16 4.49 × 10−16

FP64 einsum 8 2.56 × 10−12 OzII-Composed 8 1.29 × 10−11 OzII-NoRescale 8 8.27 OzII-Rescale 8 4.06 × 103 OzII-RescaleBE 8 6.79 × 10−10 GEMMul8-Composed 8 1.13 × 10−12 GEMMul8-KRON 8 5.27 × 10−12

3.82 × 10−16 5.39 × 10−16 1.48 × 10−3 1.00 1.56 × 10−13 2.67 × 10−16 8.04 × 10−16

4.67 × 10−12 5.28 × 10−12 7.10 × 101 4.54 × 101 5.35 × 10−8 2.45 × 10−12 2.39 × 10−11

8.56 × 10−16 1.24 × 10−15 6.18 × 10−3 1.00 2.00 × 10−12 5.58 × 10−16 1.81 × 10−15

5.79 × 10−11 2.84 × 10−11 3.53 × 101 1.00 9.96 × 10−8 1.13 × 10−11 1.86 × 10−11

3.78 × 10−16 6.21 × 10−16 2.27 × 10−3 1.00 2.27 × 10−12 3.15 × 10−16 5.23 × 10−16

from 1.87 × 10−16 to 1.54 × 10−15 , compared with approximately 10−9 to 10−8 for the two other persistent variants. For Uniform, its median error is comparable to those of FP64 einsum and the composed implementations. The larger medians on Smooth and Wide160 indicate greater sensitivity to the input distribution. At d = 8, the median errors of OzII-RescaleBE increase to 1.56 × 10−13 , 2.00 × 10−12 , and 2.27 × 10−12 on Uniform, Smooth, and Wide160, respectively. Although larger than those of FP64 einsum and the composed implementations, these errors remain many orders of magnitude below those of OzII-NoRescale and OzII-Rescale. The former has median errors between 1.48 × 10−3 and 6.18 × 10−3 , while the latter has a median of 1.00 on every dataset. Thus, base extension substantially improves the accuracy of persistent execution on the fixed residue basis, although rescaling still limits accuracy as the chain becomes deeper. Maximum elementwise errors are substantially larger than the medians for every method. For example, on Smooth at d = 3, OzII-RescaleBE has a maximum error of 4.52 × 10−9 and a median of 1.47 × 10−15 . FP64 einsum also shows a large separation, with corresponding errors of 9.35 × 10−11 and 4.33 × 10−16 . These results are consistent with the sensitivity of elementwise relative error to small reference values. Consequently, high normwise accuracy can coexist with much larger relative errors in individual output entries. 6.3

Performance Evaluation

Execution times are measured using CUDA events over multiple consecutive evaluations of the complete chain to reduce the effect of timing noise. The total elapsed time is divided by the number of repetitions to obtain the mean execution time. The OzII implementations and FP64 einsum use CUDA events through PyTorch, while GEMMul8 provides an internal CUDA-event timer for its library operations. Device memory usage is measured through PyTorch’s interface to cudaMemGetInfo. 6.3.1

Throughput

The throughput tests keep the total number of degrees of freedom (DoFs) approximately constant while varying the tensor mode size. Each scalar tensor entry counts as one element-local DoF. The tests use Uniform inputs, chain depth d = 3, and n ∈ {8, 12, 16, 32, 64, 128, 256}. The batch size is ⌊224 /n3 ⌋, giving 32768, 9709, 4096, 512, 64, 8, and 1 tensors, respectively. Each evaluation of the complete chain therefore processes 224 DoFs, except at n = 12, where it processes 224 − 64. Throughput is calculated by dividing this DoF count by the elapsed time for the complete chain, including tensor encoding and reconstruction. Operator preparation and host–device transfers are excluded. Figure 7 shows the throughput of all compared methods. FP64 einsum serves as the native-precision baseline. Its throughput increases with tensor mode size up to n = 32 because the smaller mode products provide insufficient work to keep the GPU busy. At these sizes, execution is dominated by launch and latency costs and by the permutation copies introduced by einsum. As n increases, cuBLAS processes larger, more efficient GEMMs. From n = 64, FP64 execution 22

A PREPRINT - O CTOBER 6, 2026

FP64 einsum OzII-Composed

OzII-RescaleBE GEMMul8-Composed

GEMMul8-Composed (fast) GEMMul8-KRON

Throughput (GDoF/s)

5

2 1 0.5

0.2 0.1 8

12

16

32

64

128

256

Tensor mode size n

Figure 7: Throughput for three successive mode products on Uniform inputs with approximately 224 element-local DoFs per application.

approaches the GPU’s double-precision throughput limit, so throughput decreases approximately in proportion to the floating-point operation count. OzII-Composed has throughput similar to that of the FP64 baseline at n = 8, but falls below it for n = 12 through 64. Its three encodings and three reconstructions impose a substantial conversion cost per DoF, while the FP64 DGEMMs become more efficient as n grows. The two methods have similar throughput at n = 128. At n = 256, OzII-Composed achieves 1.37× the throughput of FP64 einsum. This advantage is consistent with the increasing importance of arithmetic throughput, because the residue GEMMs use INT8 Tensor Cores rather than native FP64 arithmetic. At n = 8, OzII-RescaleBE is 1.44× higher than the FP64 baseline, whose small mode products incur substantial launch and latency costs. It is also 1.53× higher than OzII-Composed because it performs one encoding and one reconstruction instead of three. From n = 12 to n = 32, its throughput remains between 3.5 and 4.1 GDoF/s, while the FP64 DGEMMs become more efficient and reach their highest throughput at n = 32. The FP64 baseline is therefore 1.2–1.5× higher over these sizes. From n = 64, OzII-RescaleBE is 1.11–1.55× higher than FP64 einsum as the latter becomes increasingly constrained by FP64 arithmetic throughput. Across all tested tensor mode sizes, OzII-RescaleBE is 1.13–1.53× higher than OzII-Composed, showing the benefit of removing intermediate encodings and reconstructions. GEMMul8-KRON outperforms the FP64 baseline at n = 8 and OzII-RescaleBE at n = 8 and 12. It applies the complete chain through a single emulated GEMM with inner dimension n3 , providing a longer contraction that is better suited to the library’s GEMM implementation. Its throughput then falls rapidly because the dense operator requires n3 multiply–adds per output DoF, compared with 3n for the three mode products. The operator’s O(n6 ) storage requirement also limits its use at larger mode sizes. At n = 32, an FP64 copy of the operator occupies approximately 8.6 GB. At n = 64, it would require approximately 550 GB, exceeding the GPU’s 96 GB memory. Therefore the measurements ends up at n = 32. The throughput of GEMMul8-Composed increases with mode size over the tested range. Larger n increases the operator dimensions and contracted length, making the mode products better suited to the library’s dense GEMM implementation. Nevertheless, its throughput remains substantially below that of the FP64 baseline and the OzII implementations through n = 128. One reason is that the tested GEMMul8 configuration pads internal buffer dimensions to multiples of 256. At n = 8, padding the operator row count and contracted dimension considerably increases storage and data movement. At n = 256, these dimensions require no padding, and GEMMul8-Composed achieves 1.07× the throughput of OzII-RescaleBE, 1.21× that of OzII-Composed, and 1.65× that of the FP64 baseline. GEMMul8-Composed (fast) follows the same trend and achieves 1.07–1.10× the throughput of the standard configuration at every tested size. OzII-RescaleBE outperforms both GEMMul8-Composed configurations through n = 128, achieving 55–60× their throughput at n = 8 and 2.1–2.3× at n = 128. Its throughput is 6–12 % lower at n = 256. It is also competitive with GEMMul8-KRON at small tensor mode sizes and supports sizes beyond n = 32, where the dense operator becomes 23

A PREPRINT - O CTOBER 6, 2026

Encoding

Residue products

Rescaling + base extension

Reconstruction

Other

7 250

Time per application (ms)

Time per application (ms)

6.08 6 5 3.91

4

3.32

3 2 1 0

OzII-Composed OzII-RescaleBE GEMMul8-KRON

233.62 210.48

200

150

100

50

0

GEMMul8Composed (fast)

GEMMul8Composed

Figure 8: Execution-time breakdown at n = 8, d = 3, and s = 15 for a batch of tensors containing a total of 224 element-local DoFs.

impractical. Compared with the FP64 baseline, OzII-RescaleBE achieves higher throughput at n = 8 and from n = 64 onward. 6.3.2

Execution-time Breakdown

To better understand the preceding throughput results, execution-time breakdowns are examined at n = 8 and 256. Figure 8 shows the breakdown for the five compared methods at n = 8. The measurements are grouped into five categories. Encoding includes scale discovery and conversion to residues. Residue products include matrix multiplication and modular reduction, which are fused in the OzII product kernels. For GEMMul8, these operations are measured separately but combined in the figures. Rescaling and base extension form one category. Reconstruction includes recovery from residues and inverse scaling to FP64. For GEMMul8, Other includes unfolding copies and the difference between the measured library-call duration and the sum of its internally measured phases. Cached operator preparation is excluded. For OzII-Composed, encoding and reconstruction dominate, accounting for 4.73 ms of the 6.08 ms total and leaving approximately 22% for the residue products. OzII-RescaleBE removes the intermediate encodings and reconstructions, reducing their combined cost from 4.73 to 1.38 ms. However, rescaling and base extension add 1.21 ms. Together with the residue products, these operations account for approximately 65% of its execution time. GEMMul8-KRON also encodes and reconstructs only once, with these two steps taking 1.85 ms in the dense formulation. Because it requires no rescaling or base extension, its total execution time is 0.59 ms shorter than that of OzII-RescaleBE. The two GEMMul8-Composed configurations are placed in a separate panel with its own scale because their execution times are much larger. Residue products dominate the standard configuration, accounting for approximately 156 ms of its 233.62 ms total. Encoding is the second-largest cost at approximately 60 ms, which the fast scaling path reduces to 43 ms. Figure 9 compares the four remaining methods at n = 256. Residue products dominate both OzII implementations, accounting for approximately 9.2 ms of the 13.2 ms total for OzII-Composed and 9.2 ms of the 11.6 ms total for OzII-RescaleBE. OzII-Composed still spends approximately 31% (4.05 ms) of its execution time on three encodings and three reconstructions. In contrast, OzII-RescaleBE spends only 21% on encoding, reconstruction, rescaling, and base extension combined. For the two GEMMul8-Composed configurations, the distribution is reversed. Encoding and reconstruction now dominate, accounting for 5.4–5.9 ms, or approximately 52–53% of the total execution time. Although encoding and reconstruction are substantially less costly in OzII-RescaleBE, the custom residue GEMM kernels in the current implementation become the bottleneck at large n. This suggests that, for large contracted dimensions, using cuBLAS INT8 GEMMs followed by a separate modular reduction pass, as GEMMul8 does, could improve performance. 24

A PREPRINT - O CTOBER 6, 2026

Encoding

Time per application (ms)

14

Residue products

Rescaling + base extension

Reconstruction

Other

13.19 11.62

12

10.26

11.11

10 8 6 4 2 0

OzII-Composed

OzII-RescaleBE

GEMMul8Composed (fast)

GEMMul8Composed

Device memory (GB)

Figure 9: Execution-time breakdown at n = 256, d = 3, and s = 15 for one tensor containing 224 element-local DoFs. FP64 einsum

GEMMul8-Composed

OzII-Composed

GEMMul8-Composed (fast)

OzII-RescaleBE

GEMMul8-KRON

4 3 2 1 0

n=8

n = 256

Figure 10: Device-memory footprint for tensor mode products on 224 element-local DoFs with d = 3, n ∈ {8, 256}, and s = 15.

6.3.3

Memory Usage

Figure 10 reports the device-memory footprint, measured as the increase in used device memory from before input allocation to after setup, one warm-up execution, and synchronization. It includes the input and output tensors (0.268 GB together), the method’s retained buffers, and cached allocations and workspaces, but does not capture transient peaks during execution. At n = 8, FP64 einsum occupies 0.539 GB, including approximately 0.27 GB for two FP64 intermediates retained by PyTorch’s caching allocator. OzII-Composed occupies 0.975 GB, while OzII-RescaleBE occupies 0.656 GB, a reduction of 32.7 %. GEMMul8-KRON occupies 0.965 GB, and the two GEMMul8-Composed configurations occupy 4.40 and 4.42 GB. The GEMMul8-Composed footprint at n = 8 results from the library’s padding and workspace allocation. Each mode product is a GEMM with dimensions M = K = n and a column count equal to the batch size times n2 . At n = 8, this gives M = K = 8 and 2,097,152 columns. GEMMul8 rejects GEMMs with 218 or more columns, so the implementation splits the unfolding into 16 chunks of 131,072 columns, corresponding to 2,048 tensors each, and calls 25

A PREPRINT - O CTOBER 6, 2026

the library once per chunk. The library pads M and K from 8 to 256. For a fixed chunk, a buffer that scales linearly with one of these padded dimensions is therefore 32 times larger than its unpadded counterpart. One workspace per mode product is retained to preserve the cached operator scaling. The library reports 1.18 GB of workspace per mode product at n = 8, compared with 0.62 GB at n = 256. At n = 256, the footprint of FP64 einsum is unchanged. OzII-Composed occupies 0.797 GB, while OzII-RescaleBE occupies 0.656 GB, a reduction of 17.7 %. Both GEMMul8-Composed configurations occupy 2.743 GB, which is 4.18× the footprint of OzII-RescaleBE. The footprints of OzII-Composed and GEMMul8-Composed decrease with increasing mode size for different reasons. OzII-Composed scales each contracted fiber separately, so its scale arrays and scale-discovery temporaries grow in proportion to the number of fibers. At n = 8, there are approximately two million fibers across 32,768 tensors, each containing 64 fibers, requiring about 17 MB per FP64 scale array. At n = 256, there are only 65,536 fibers, requiring approximately 0.5 MB per array. OzII-RescaleBE uses one scale per tensor, and its measured footprint remains unchanged. For GEMMul8-Composed, the footprint decreases because the GEMM dimensions at n = 256 are already multiples of 256 and require no padding.

7

Conclusion

This work presented OzII-RescaleBE, a persistent execution method for applying Ozaki scheme II to chained tensor mode products. The method combines rescaling and base extension to control integer magnitude and restore representation capacity between products. It preserves the tensor mode-product structure while encoding the input once, retaining intermediate tensors in residue form, and reconstructing only the final result. Successive mode products can therefore be evaluated without converting intermediate tensors back to FP64. An implementation using CuTe DSL and INT8 Tensor Cores was provided for verification and performance evaluation. Across three datasets with different input distributions, OzII-RescaleBE achieved normwise relative errors ranging from 2.5 × 10−16 to 2.6 × 10−15 for chain depths d = 3 through 7. The corresponding errors ranged from 10−16 to 10−15 for GEMMul8-KRON and from 1 × 10−16 to 3 × 10−16 for GEMMul8-Composed and the FP64 chain. At d = 8, the errors of OzII-RescaleBE increased to between 7 × 10−14 and 6 × 10−13 . OzII-RescaleBE achieved throughputs ranging from 1.47 to 4.42 GDoF/s across mode sizes n = 8 to 256. Its throughput was within 12 % of GEMMul8-KRON at n = 8 and 12 and higher at the larger sizes where both methods were evaluated. OzII-RescaleBE supported every tested size, whereas GEMMul8-KRON was limited to n ≤ 32 by its O(n6 ) operator storage requirement. OzII-RescaleBE also achieved 2.1–60× the throughput of the GEMMul8Composed configurations through n = 128. The largest gains reflect the library’s padding overhead for small contracted dimensions. At n = 256, its throughput was 6–12 % lower because the custom residue GEMM kernels were slower than the library’s cuBLAS-based implementation. Compared with the FP64 baseline, OzII-RescaleBE achieved higher throughput at n = 8 and from n = 64 onward. The execution-time breakdown explains these differences. At n = 8, rescaling and base extension introduced additional costs that outweighed the conversion savings relative to GEMMul8-KRON. At n = 256, residue products accounted for approximately 79 % of the OzII-RescaleBE execution time. Its device-memory footprint at n = 8 was approximately 32 % smaller than that of GEMMul8-KRON, while GEMMul8-Composed required approximately 6.7× as much memory. At n = 256, GEMMul8-Composed required approximately 4.2× as much memory, with its retained workspaces for individual mode products contributing to the difference even without padding. Beyond the performance limitations of the current implementation, the method has two main limitations. First, rescaling introduces rounding errors, as reflected in the increased errors at depth eight. A larger modulus basis may permit less aggressive rescaling, but also increases storage and residue-product work. Second, the current formulation uses a single encoding scale for each input tensor. When tensor entries span a wide range of magnitudes, small entries can lose significant bits or round to zero during the initial encoding. Unlike per-fiber encoding in OzII-Composed, this shared scale cannot adapt to the magnitude distribution of individual fibers. Any resulting information loss occurs before rescaling and cannot be recovered by base extension. Future work includes improving the residue-product kernels and exploring shared-memory execution that retains intermediate residues on chip for small tensors. A native CUDA implementation will support integration into DG simulations, followed by evaluation of representative operators over a broader range of tensor sizes, chain depths, and modulus counts. Alternatives to greedy modulus dropping will also be investigated to account for later mode products and balance rounding error, representation capacity, and execution cost. 26

A PREPRINT - O CTOBER 6, 2026

Acknowledgements This work is partially supported by MEXT as "Project for Establishment of a Center for Advanced HPC-AI Development Support". Generative AI tools were used under the authors’ direction. Claude Code (Fable 5.1) assisted with the implementation, the evaluation scripts and measurements, the figures, and the drafting of parts of the text. GPT-6 Astra generated initial manuscript drafts and assisted with revision. The authors directed the investigation, verified the AI-generated material against exact references and independent measurements, and take responsibility for the scientific content.

27

A PREPRINT - O CTOBER 6, 2026

References [1] NVIDIA. NVIDIA H100 GPU: Product specifications. Manufacturer product page. H100 SXM specifications; accessed 2026-09-11. [2] NVIDIA. NVIDIA Blackwell Ultra: Datasheet. Manufacturer datasheet. Technical specifications, p. 5; accessed 2026-09-11. [3] Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. Acta Numerica, 31:347–414, 2022. [4] Katsuhisa Ozaki, Takeshi Ogita, Shin’ichi Oishi, and Siegfried M. Rump. Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numerical Algorithms, 59(1):95– 118, 2012. [5] Daichi Mukunoki, Katsuhisa Ozaki, Takeshi Ogita, and Toshiyuki Imamura. DGEMM using tensor cores, and its accurate and reproducible versions. In High Performance Computing, pages 230–248. Springer, 2020. [6] Hiroyuki Ootomo, Katsuhisa Ozaki, and Rio Yokota. DGEMM on integer matrix multiplication unit. The International Journal of High Performance Computing Applications, 2024. [7] Katsuhisa Ozaki, Yuki Uchino, and Toshiyuki Imamura. Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique. The International Journal of High Performance Computing Applications, 2026. [8] RIKEN R-CCS. GEMMul8: GEMM emulation using low-precision matrix engines based on Ozaki Scheme II. Software repository. Accessed 2026-09-11. [9] Kasia Świrydowicz, Noel Chalmers, Ali Karakus, and Tim Warburton. Acceleration of tensor-product operations for high-order finite element methods. The International Journal of High Performance Computing Applications, 33(4):735–757, 2019. [10] Denghui Lu, Alexander Maeder, Mathieu Luisier, and Alexandros Nikolaos Ziogas. EmuGEMM: Fused tensor core kernels for precision emulation in matrix multiplication, 2026. [11] Jung Hee Cheon, Andrey Kim, Miran Kim, and Yongsoo Song. Homomorphic encryption for arithmetic of approximate numbers. In Advances in Cryptology – ASIACRYPT 2017, volume 10624 of Lecture Notes in Computer Science. Springer, 2017. [12] Jung Hee Cheon, Kyoohyung Han, Andrey Kim, Miran Kim, and Yongsoo Song. A full RNS variant of approximate homomorphic encryption. Cryptology ePrint Archive, Paper 2018/931, 2018. The ePrint record identifies SAC 2018 publication. [13] Shai Halevi, Yuriy Polyakov, and Victor Shoup. An improved RNS variant of the BFV homomorphic encryption scheme. Cryptology ePrint Archive, Paper 2018/117, 2018. Revised 2019-03-06; published at CT-RSA 2019, DOI 10.1007/978-3-030-12612-4_5. [14] Tamara G. Kolda and Brett W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):455–500, 2009. [15] H. L. Garner. The residue number system. IRE Transactions on Electronic Computers, EC-8(2):140–147, June 1959. [16] W. G. Horner. A new method of solving numerical equations of all orders, by continuous approximation. Philosophical Transactions of the Royal Society of London, 109:308–335, 1819. [17] Cu Cui. Acceleration of tensor-product operations with tensor cores. ACM Transactions on Parallel Computing, 11(4), 2024. [18] Jiqun Tu, Ian Karlin, John Camier, Veselin Dobrev, Tzanio Kolev, Stefan Henneking, and Omar Ghattas. Accelerating high-order finite element simulations at extreme scale with FP64 tensor cores, 2026. [19] Daichi Mukunoki. DGEMM without FP64 arithmetic – using FP64 emulation and FP8 tensor cores with Ozaki scheme, 2025. [20] Yuki Uchino, Katsuhisa Ozaki, and Toshiyuki Imamura. Double-precision matrix multiplication emulation via Ozaki-II scheme with FP8 quantization, 2026. [21] Shun-ichiro Hayashi, Daichi Mukunoki, Tetsuya Hoshino, and Takahiro Katagiri. DGEMM with Ozaki Scheme I/II on FP4 tensor cores: A base-13 E2M1 limb representation, 2026. 28

A PREPRINT - O CTOBER 6, 2026

[22] Angelika Schwarz, Anton Anders, Cole Brower, Harun Bayraktar, John Gunnels, Kate Clark, RuQing G. Xu, Samuel Rodriguez, Sebastien Cayrols, Paweł Tabaszewski, and Victor Podlozhnyuk. Guaranteed DGEMM accuracy while using reduced precision tensor cores through extensions of the Ozaki scheme, 2025. Also published at SCA/HPCAsia 2026, DOI 10.1145/3773656.3773670. [23] Yuki Uchino, Qianxiang Ma, Toshiyuki Imamura, Katsuhisa Ozaki, and Patrick Lars Gutsche. Emulation of complex matrix multiplication based on the Chinese Remainder Theorem, 2025. Inspected preprint version 3. The arXiv record also identifies ISC High Performance 2026 proceedings, DOI 10.23919/isc.2026.11520500. [24] Satoshi Matsuoka. FP8 is all you need (part 2): Full-FP64 3-d FFT on FP8-generation tensor cores: The integer-epilogue wall and the minimal hardware that would remove it, 2026. Version 3, 8 September 2026. [25] Jean-Guillaume Dumas, Pascal Giorgi, and Clément Pernet. Dense linear algebra over word-size prime fields: The FFLAS and FFPACK packages. ACM Transactions on Mathematical Software, 35(3):1–42, 2008. [26] Cansu Demirkiran, Guowei Yang, Darius Bunandar, and Ajay Joshi. Mirage: An RNS-based photonic accelerator for DNN training. In 2024 ACM/IEEE 51st Annual International Symposium on Computer Architecture (ISCA), 2024. [27] Jean-Claude Bajard, Julien Eynard, Anwar Hasan, and Vincent Zucca. A full RNS variant of FV like somewhat homomorphic encryption schemes. Cryptology ePrint Archive, Paper 2016/510, 2016. Extended version of the SAC 2016 paper. [28] NVIDIA Corporation. CuTe DSL 4.5.2 warp submodule, 2026. Accessed September 14, 2026. [29] NVIDIA Corporation. PTX ISA, warp-level matrix multiply-accumulate instructions. Accessed September 16, 2026. [30] NVIDIA. NVIDIA RTX PRO Blackwell GPU Architecture. Architecture white paper, version 1.1. FP64 throughput ratio on p. 8 and Max-Q specifications in Table 4 on pp. 46–47. Accessed 2026-09-15.

29

A PREPRINT - O CTOBER 6, 2026

A

Additional comparison methods

This appendix describes the two additional persistent Ozaki-II variants and the FP64 einsum used in the evaluation. OzII-NoRescale performs the chain without intermediate rescaling. OzII-Rescale drops moduli through rescaling but does not restore them. Both encode the input tensor once and decode only the final result. They follow the effective-operator preparation and reconstruction conventions of the main persistent method. A.1

FP64 einsum baseline

The FP64 einsum applies one torch.einsum contraction per mode product to CUDA tensors stored in torch.float64. Each call applies the corresponding n × n operator along the specified mode, preserving the batch dimension and the other two tensor indices. Its FP64 output becomes the input to the next mode product, following the same operator order and mode sequence as the residue implementations. The baseline does not form a dense Kronecker operator or perform residue encoding. It provides a finite-precision execution baseline, distinct from the exact rational-chain reference whose final result is converted to FP64 for error measurement. A.2

OzII-NoRescale

For a candidate significand bits b, let Āi (b) denote the effective integer operator for the i-th mode product. Its maximum absolute row sum is n X ρi (b) = max |(Āi (b))jk |. (33) 1≤j≤n

k=1

Each output entry is a weighted sum along one contracted fiber, so multiplying the current magnitude bound by ρi (b) gives a bound for the output entries. The effective operators and their bounds are recomputed for each candidate b. The input uses one encoding scale per tensor, giving the initial reconstruction scale α = 2−e . The field Γ contains the remaining operator reconstruction scales. When a mode is revisited, its pending scales are folded into the next effective operator before encoding, as in the main methodology. Thus, Γ contains only the scales still to be applied at final reconstruction. Q OzII-NoRescale retains the initial basis P throughout the chain, with P = p∈P p. Without rescaling, the integer magnitude after the j-th mode product is bounded by Bj (b) = 2

b

j Y

ρi (b),

B0 (b) = 2b .

(34)

i=1

The planner selects the largest b ∈ {1, . . . , 53} satisfying 2Bj (b) < P for every j = 0, . . . , d. When all ρi (b) ≥ 1, Qd the bounds are nondecreasing, and this condition reduces to the whole-chain test 2b i=1 ρi (b) < P/2. The range requirement can reduce the selected encoding bit count as depth increases, even though the number of residue planes remains constant. The method retains the supplied s moduli and does not require the entire available modulus pool. A.3

OzII-Rescale and its planner

OzII-Rescale starts from the same basis but removes moduli after nonfinal mode products. Algorithm 9 calls DropPlan with restoration disabled. The flag τ = 0 corresponds to the S HRINK mode in the implementation. In contrast, τ = 1, used by PlanRescaleBE, assumes that base extension restores the working basis before the next mode product. The planner starts at b = 53 and decreases b until the range checks pass for every product and planned drop. It returns the selected significand bits and the drop counts (gi ), where gi specifies the number of drops after the i-th mode product and gd = 0. If no positive b passes the checks, it reports failure. No base extension is performed, so removed moduli remain unavailable for subsequent products. Algorithm 10 maintains the available basis L, initialized to P. At each planned drop, RnsDrop divides the represented integer by the smallest available modulus q with rounding, removes that modulus, and updates the reconstruction scale to α′ = αq. The returned scale becomes α for the next operation. At the end of the chain, α therefore contains the initial tensor reconstruction scale multiplied by every dropped modulus. The final reconstruction uses the remaining basis and the scale αΓ. Rescaling reduces both the integer magnitude and the available representation range. Without base extension, that lost representation capacity is not restored.

30

A PREPRINT - O CTOBER 6, 2026

Algorithm 8 OzII-NoRescale with a fixed basis throughout the chain Require: Tensor X, operators (Ai )di=1 , modes (mi )di=1 Require: Initial basis P with s = |P| 1: procedureQO Z IIN O R ESCALE(X, P, (Ai ), (mi )) 2: P ← p∈P p 3: Define (ρi (b)) from effective operator encoding 4: Choose the largest b ∈ {1, . . . , 53} with 2Bj (b) < P for all j = 0, . . . , d 5: if no such b exists then 6: return ⊥ 7: end if 8: Prepare effective operators (Āi ) and scale field Γ at b 9: Cache Āi,p = [Āi ]p for p ∈ P 10: (X̄, e) ← E NCODE(X, b, {X}) 11: α ← 2−e 12: R ← T O R NS(X̄, P) 13: for i = 1, . . . , d do 14: for all p ∈ P do 15: Rp ← [Rp ×mi Āi,p ]p 16: end for 17: end for 18: return R ECONSTRUCT FP64(R, P, αΓ) 19: end procedure

Algorithm 9 PlanRescale for precision selection without base extension Require: Operator row-sum bound functions (ρi ), initial basis P Ensure: Encoding parameter b and drop counts g, or failure ⊥ 1: procedure P LAN R ESCALE((ρi ), P) 2: for b = 53, 52, . . . , 1 do 3: g ← D ROP P LAN(b, (ρi (b)), P, 0) 4: if g ̸= ⊥ then 5: return (b, g) 6: end if 7: end for 8: return ⊥ 9: end procedure

31

A PREPRINT - O CTOBER 6, 2026

Algorithm 10 OzII-Rescale with a shrinking residue basis Require: Tensor X, operators (Ai )di=1 , modes (mi )di=1 Require: Initial basis P with s = |P| 1: procedure O Z IIR ESCALE(X, P, (Ai ), (mi )) 2: Define (ρi (b)) from effective operator encoding 3: π ← P LAN R ESCALE((ρi ), P) 4: if π = ⊥ then 5: return ⊥ 6: end if 7: (b, (gi )) ← π 8: Prepare effective operators (Āi ) and scale field Γ at b 9: Cache Āi,p = [Āi ]p for p ∈ P 10: (X̄, e) ← E NCODE(X, b, {X}) 11: α ← 2−e 12: L←P 13: R ← T O R NS(X̄, L) 14: for i = 1, . . . , d do 15: for all p ∈ L do 16: Rp ← [Rp ×mi Āi,p ]p 17: end for 18: for j = 1, . . . , gi do 19: q ← min L 20: (R, L, α′ ) ← R NS D ROP(R, L, q, α) 21: α ← α′ 22: end for 23: end for 24: return R ECONSTRUCT FP64(R, L, αΓ) 25: end procedure

32

Related documents

• You Only Convert Twice: Ozaki Scheme II for Chained Tensor Mode Products#1201072
• DGEMM with Ozaki Scheme I/II on FP4 Tensor Cores: A Base-13 E2M1 Limb Representation#440342
• Improved Scaling for Fast Mode of Ozaki Scheme II#321809
• MoRo: From One-Sample Ring-LWE Rounding Key Exchange to Module-LWE IND-CCA KEM.#317971
• A simpler kernel for stratified Mukai flops#1296985