1 Introduction↩︎

Multi-precision integer arithmetic is a basic component of computer algebra, cryptography, exact scientific computation, and symbolic-numeric software. Its performance matters not only for isolated large computations, but also for applications that require many independent integer operations at the same precision. This makes GPUs attractive: they offer high arithmetic throughput and massive parallelism, provided that the computation can be organized to keep data movement and inter-thread communication under control.

Existing GPU support for multi-precision integers is strongest at relatively small precisions, where one arithmetic instance can be mapped to a small cooperative group of threads. The nvidia cgbn library (Cooperative Groups Big Numbers) is an important example of this approach [1]. Such libraries provide very high performance in their intended range, but they do not cover all practically interesting sizes. In particular, there is a middle range of integer precisions—large enough that warp-level methods become strained, but still small enough that a complete arithmetic instance can fit within the fast memories of a single GPU thread block. This paper is concerned with this midsize regime.

In previous work we studied GPU implementations of multi-precision addition and multiplication in this regime [2]. That work showed that the classical algorithms, when scheduled carefully, can be made effective on GPUs by assigning one multi-precision operation to a cuda block, keeping operands and intermediate values in registers or shared memory, and minimizing global-memory traffic. It also showed that such algorithms can be expressed in a high-level data-parallel language such as Futhark, although some low-level transformations needed for peak performance remain beyond the current compiler.

Division is a more demanding operation. It is not a simple local operation on digits, and the usual high-performance approach is to reduce division to multiplication by first computing an approximation to the reciprocal of the divisor. In conventional Newton iteration this typically requires working in a domain where such reciprocals exist, which can introduce multiple precision floating-point approximations and interactions between different arithmetic domains. For exact integer arithmetic this is undesirable: the implementation must preserve exactness while still exposing enough parallelism for the GPU.

The algorithm of Watt [3] addresses this problem by replacing the reciprocal with a whole shifted inverse. Instead of computing \(1/v\), it computes an integer approximation to \(B^h/v\), where \(B\) is the digit base and \(h\) is an appropriate precision. Newton iteration can then be formulated using integer multiplication and shift operations. This is attractive for GPUs because the main operations are data-parallel: multi-precision multiplication, addition, subtraction, comparison, and shifts. Moreover, the algorithm is parameterized by the multiplication method, so it can in principle benefit from either classical multiplication or faster multiplication algorithms.

This paper investigates how this shifted-inverse division algorithm can be implemented efficiently on GPUs for midsize integers. Our present implementation is written directly in cuda rather than generated from Futhark. This is intentional: the goal is to identify the low-level scheduling, storage, and specialization issues that a high-level compiler would ultimately have to handle. In particular, division stresses the compiler more than addition or multiplication because the Newton refinement uses operations whose effective precisions change during the iteration. The main contributions of the paper are as follows:

• We present a cuda implementation of multi-precision integer division based on the whole-shifted-inverse algorithm of Watt [3], targeting integer precisions from \(2^{15}\) to \(2^{18}\) bits.

• We describe implementation refinements needed to make the algorithm robust in an unsigned-integer setting, including explicit sign handling in the close-product computation and quotient correction when the computed shifted inverse may overestimate by one.

• We give a cost model for the implementation in terms of the number of full multi-precision multiplications required. For the classical multiplication used here, the model predicts that the full division operation should require at least five and at most seven full multiplication costs.

• We show how the supporting operations—shifts, comparisons, subtractions, close products, and variable-size multiplications—are mapped to cuda blocks using registers and shared memory.

• We evaluate the implementation on an nvidia A100 GPU and compare it with cgbn’s corresponding operation. The results show that, at the largest tested precisions, the measured division time is close to the five-multiplication lower bound predicted by the cost model, while also covering precision ranges not supported by cgbn division in our experiments.

The rest of the paper is organized as follows. Section 2 reviews the shifted-inverse division algorithm, including the refinements needed for the implementation and the resulting multiplication-based cost model. Section 3 describes the cuda implementation, focusing on the use of registers, shared memory, scans, shifts, subtraction, and variable-size multiplication. Section 4 presents the experimental evaluation and compares the implementation with cgbn where the corresponding operations are supported.

Section 5 discusses related work on GPU multi-precision arithmetic and high-level data-parallel programming. Section 6 summarizes the results and outlines future directions, including clipped products and compiler support for this class of exact arithmetic kernels.

2 Division Algorithm↩︎

2.1 Intuition↩︎

Parallel implementations for division of multi-precision integers typically rely on Newton’s method to compute the reciprocal of the divisor, by applying it to function \(f(x)=1/x-v=0\), where \(v\) is the divisor and \(x=\frac{1}{v}\) its reciprocal: \[x_{(i+1)}=x_{(i)} - \frac{f(x_{(i)})}{f'(x_{(i)})}=x_{(i)}+x_{(i)} \cdot \left(1-v\cdot x_{(i)}\right) \label{eq:newton}\tag{1}\] However, Newton’s method in its general form requires working in a related domain where the reciprocal exists. This can lead to a complex library structure in which the arithmetic domains are interdependent. As well, internal floating point representation can lead to potential loss of precision and overhead when converting between domains. The algorithm proposed by Watt [3] addresses these shortcomings and allows computation to be carried out in the integral domain, essentially by applying the Newton method to \(f(x)=B^h/x-v=0\), where \(B\) is the base of multi-precision integer \(v\) and \(h\) its precision, i.e., \(B^h \geq v\). This computes the “whole shifted inverse” \(x = \lfloor B^h/v \rfloor\) with the Newton iteration: \[w_{(i+1)} = w_{(i)}+\left\lfloor w_{(i)} \cdot (B^h-v \cdot w_{(i)}) \cdot B^{-h}) \right\rfloor, \quad w_{(i)} \in \mathbb{Z} \label{eq:tmp}\tag{2}\] 2 provides the intuition that the resulting Newton iteration can be written in terms of multiplications and efficient shift operations: \[w_{(i+1)} = w_{(i)}+\mathop{\mathrm{shift}}_{-h}(\mathop{\mathrm{shift}}_h(w_{(i)})-v \cdot w_{(i)}^2), \quad w_{(i)} \in \mathbb{Z} \label{eq:newton95shift}\tag{3}\] i.e., (1) shift operations are used to scale the input to a degree that captures all the necessary information from the fractional counterpart, (2) then computations are performed at this higher scale, and (3) a corresponding inverse shift is applied at the end to restore the original magnitude, avoiding potential loss of information. The shift and shift inverse operations are formally defined as:

When \(n \geq 0\), \(\mathop{\mathrm{shift}}_{n,B}(u)\) corresponds to integer multiplication, i.e., \(u \cdot B^n\). When \(n < 0\) it is instead a specialized quotient operation, with \(u\) as dividend and \(B^n\) as divisor. Using a multi-precision integer representation, this is equivalent to an arithmetic shift, e.g. \(\mathop{\mathrm{shift}}_1([1,2,3])=[0,1,2]\) and \(\mathop{\mathrm{shift}}_{-1}([1,2,3])=[2,3,0]\).

In contrast, the whole shifted inverse \(\mathop{\mathrm{shinv}}_n(v)\) corresponds to a specialized reciprocal, i.e., a reciprocal that has been shifted into our domain, e.g., \(\mathop{\mathrm{shinv}}_3(8)=\mathop{\mathrm{shift}}_3(0.125)=125\).

Theorem 1 shows how these operations are used to derive the quotient: applying a reverse shift on the result of multiplying \(u\) with the whole shifted inverse of \(v\) produces a result that is at most one unit away from the correct quotient.

2.2 Original Algorithm and New Revisions↩︎

The pseudocode for computing the whole-shifted inverse is recounted in Algorithm 1, which uses multi-precision addition, subtraction and shift operations, and is defined in terms of a generic multi-precision multiplication method, denoted by mult. As well, it uses the PowDiff function, defined in Algorithm 2.

The algorithm’s correctness, fast convergence and work asymptotics are proven in [3]. In short, the algorithm has asymptotic work equal to one multi-precision multiplication; this holds even when Strassen’s \(O(n\cdot \text{log}~n)\) algorithm [4] for multiplication is used. The algorithm consists of three main stages: special case handling, initial approximation and iterative refinement, which we summarize below:

Special Case Handling (lines [g1-start]-[g1-end]) ensures that the easy cases—corresponding to \(B<v\leq B^h/2\)—are handled, and that the base is sufficiently large (\(B>=16\)), such that the prerequisites for a good initial-value choice are met.

Initial Approximation (lines [g2-start]-[g2-end]). The original algorithm uses a three-digit approximation of the original \(v\) value, namely \(V = \sum_{i=0}^{\ell} v_{k - \ell + i} \cdot B^i\), where \(\ell = min(k,2)\). The shifted inverse is approximated to \(w = \left(B^{2\ell} - V\right) \mathbin{\mathrm{quo}}V + 1\), which is more convenient and faster to compute than the equivalent \(B^{2l} \mathbin{\mathrm{quo}}V\), albeit both take constant time. If sufficiently many digits are considered correct, i.e., \(h - k \leq \ell\), then the original algorithm shifts \(w\) to the appropriate magnitude and performs an early return, i.e., \(\mathrm{\small shift}_{h - k - \ell}(w)\), otherwise \(w\) is refined.

We perform the revision shown in lines [g2-start]-[g2-end] of Algorithm 1 because there are rare corner cases in which the early return refers to an overestimated (unsafe) value (see  [5] for an example). Instead, we always pass the initial approximation through Refine, which guarantees correctness. This also allows a less-precise initial approximation of \(v\), namely \(V = v_{k-1} + v_k\cdot B\) that uses only two (instead of three) digits, i.e., \(\ell=2\), which is proven in [5] to preserve fast convergence and is convenient since it promotes machine arithmetic.

Iterative Refinement (functions Refine, Step, PowDiff). The initial approximation is refined at least once with the fastest-convergence routine—named Refine3 in [3] and Refine here—that employs both shorter-iterates and divisor-prefixes techniques to achieve optimal work. Refine iteratively calls the Step function, which performs a single Newton iteration. Step invokes PowDiff (shown in Algorithm 2), which computes \(B^h - v \cdot w\) using the close-product strategy for improved efficiency. Since PowDiff can return negative integers and since our implementation assumes unsigned integers, we refine the original implementation to explicitly keep track of the integers’ sign in PowDiff and Step.

With divisor prefixes [3], overestimation can also arise during refinement: when Refine uses only a prefix of \(v\), rare edge cases allow low digits to affect higher ones. It is proven [5] that Algorithm 1 can then overestimate the shifted inverse by at most one, so \(\texttt{shift}_{-h}(u \cdot \texttt{shinv}_{h}(v))\) can be one below or one above \(u \mathbin{\mathrm{quo}}v\). We therefore revise Theorem 1 and amend Algorithm 3:

2.3 Algorithm Cost in Number of Full Multiplications↩︎

This section assumes multi-precision integers consisting of \(M\) digits in base \(B\) and approximates the cost of the division algorithm in terms of the minimal and maximal number of full multiplications that are performed with the classical/quadratic algorithm. We consider that a full multiplication is performed whenever the result requires to compute more than \(\frac{M}{2}\) digits. In summary, the computation of the shifted inverse requires at least two and at most four full multiplications. Once the shifted inverse is known, the straightforward computation of the quotient and remainder shown in Algorithm 3 requires:

• a full multiplication \(v\cdot q\) at line [line:rem-init], just before the quotient adjustment, and

• a multiplication \(u\cdot \mathop{\mathrm{shinv}}(v,h,B)\) at line [line:quot-init] in the computation that approximates the quotient. This multiplication has to be computed in double precision \(2\cdot M\) because the result is shifted back by \(h\) digits. Assuming classical/quadratic multiplication, its cost is thus equal to two full multiplications.

It follows that the presented division algorithm requires at least five and at most seven full multiplications. The remainder of this section justifies the lower and upper bound of the cost of the whole shift inverse algorithm. Essentially, the loop inside the Refine function of Algorithm 1 exhibits at least one and at most two full multiplications inside PowDiff (called from Step and shown in Algorithm 2) and similarly for the computation of Step excluding PowDiff.

A full multiplications inside PowDiff requires that \(k > \frac{h}{2}\), where \(k\) is the precision of (the original) \(v\). The precision of the \(v\) parameter of PowDiff is \(\text{prec}_B(v) = \text{min}(2\cdot \ell_i, k)\), which is also a good approximation of \(L\). It follows that a full multiplication is performed whenever \(\lfloor \frac{h}{2} \rfloor < \text{prec}_{B}(v) \approx 2\cdot \ell_i\). However, we also know that the loop in Refine executes as long as \(\ell_i < h - k\). Since \(k > \frac{h}{2}\) it follows that \(h-k \leq \lceil \frac{h}{2} \rceil\) and the loop terminates whenever \(\ell_i\) reaches \(\lceil \frac{h}{2} \rceil\).

The condition for performing a full multiplication was \(\lfloor \frac{h}{2} \rfloor < 2\cdot \ell_i\), which is equivalent to \(\lfloor \frac{h}{2} \rfloor - 1 < 2\cdot \ell_i - 1 \approx \ell_{i+1}\), since in most relevant cases the update formula for \(\ell\) is \(\ell_{i+1} = 2\cdot \ell_i - 1\). It follows that it is possible, albeit unlikely, to be in the case \(\lfloor \frac{h}{2} \rfloor - 1 < \ell_{i+1} < \lceil \frac{h}{2} \rceil\) that requires the loop in Refine to execute another iteration performing a full multiplication.

Overall, the focus of an efficient implementation of this algorithm is to achieve a runtime close to that of five full multiplications, which critically requires that the computation outside said multiplication does not introduce bottlenecks.

3 GPU Considerations↩︎

3.1 High-Level Rationale of the Implementation↩︎

We report a cuda implementation of the division algorithm of 2 for multi-precision integers whose computation fits inside one cuda block. Following prior work on addition and multiplication [2], operands are copied once from global memory to registers, results are copied back, and the remaining computation uses fast memory. Intermediate arrays stay in registers when possible and are materialized transiently in shared memory only when communication or performance requires it, for example for overlapping access in classical multiplication or as a staging buffer for coalesced global-memory transfers. In addition we apply classical techniques such as efficient sequentialization of excess parallelism to minimize inter-thread communication and thus maximize throughput.

More detailed many-core models can account explicitly for memory transactions, synchronization, occupancy, and parallelism overheads [6], [7]. Our cost model is deliberately coarser: it uses the measured cost of our multi-precision multiplication kernel as the architecture-aware unit of cost, since that kernel uses the same representation, memory hierarchy, and block-level execution strategy as the division kernel.

This section uses the following notation:

uint:

the word size representing a digit of the multi-precision integer; we use word sizes of \(16\), \(32\) or \(64\) bits, since these are hardware supported.

M:

The total number of digits in the big integer. For example, an integer with \(2^{17}\) bits could be represented using \(M = 2048\) and a 64-bit word size.

Sequentialization factor (\(Q\)):

The amount of sequential work each thread performs. For simplicity, we assume that \(Q\) evenly divides \(M\).

Memory Limitations. Currently, our cuda implementation supports integer division on integers as large as \(2^{18}\) bits. These sizes are limited by the maximum amount of shared memory available per cuda block. For example, the cuda implementation of classical multiplication [2] requires manifesting both input arrays in shared memory, which sums up to \(64\)KB. Since the current practical maximal amount of shared memory per-cuda block of our nvidia A100 GPU is about \(100\)KB, this does not permit a multi-precision size of \(2^{19}\) bits.

Another limiting factor is the amount of register memory: currently cuda supports a maximum of 255 registers per thread or \(64\)K registers per cuda block, whichever is lowest. When the register demand exceeds these bounds, the nvcc compiler resorts to register spilling [8], which allocates the excess registers in a higher level cache that is, however, significantly less efficient to access.

Our implementation uses a maximal thread-sequentialization factor \(Q=4\), since the implementation of multiplication [2] is optimized for this value—i.e., each thread computes four elements of the multiplication result. Using \(Q=4\) for the biggest multi-precision size (\(2^{18}\) bits) results in spilling \(40\) registers (\(160\) bytes) to slower storage, each of them being accessed just under three times (\(392\) bytes of spilled storage are accessed). In principle, suitably increasing \(Q\) eliminates the spilling of registers—because it decreases the number of threads in a cuda block and allows each thread to use more registers—but this did not improve the overall performance of the division implementation.

3.2 Overview of the cuda Implementation↩︎

4 gives the gist of the cuda implementation by presenting several parallel components of the division algorithm. [lst:cpy] demonstrates the manner in which the input multi-precision integers are efficiently loaded from global to register memory using shared memory as a staging buffer: The first step uses consecutive threads to copy consecutive memory locations to shared memory (shmem), thus achieving coalesced reads from global memory. Each thread performs \(Q\) such copies in the loop at lines [line:cpy-loop1-beg]-[line:cpy-loop1-end], where the cuda block size is dimensioned such that M = Q * blockDim.x. Since shared-memory accesses are not affected by un-coalesced accesses, the loop at lines [line:cpy-loop2-beg]-[line:cpy-loop2-end] copies \(Q\) consecutive elements from shared memory to each of the thread’s private (register) memory.

[lst:shift] presents the operation that shifts a multi-precision integer \(U\) by \(n\): the result of the shift is manifested in shared memory by the first loop, and once all threads have terminated work (i.e., reached the __syncthreads() barrier), the second loop loads \(Q\) consecutive elements of the result to the register (private) memory of each thread. Please note that the shared-memory buffer is only transiently utilized, and the same buffer is reused for following operations.

Several operations of the whole-shifted inverse algorithm use arguments that are powers of the corresponding base \(B\), e.g., subtracting or comparing with \(B^{bpow}\). Such specialized cases accept a more efficient implementation than when the arguments are arbitrary integers. [lst:specialized_sub] demonstrates the case of subtraction \(U - B^{bpow}\): The first loop finds the lowest index \(n\) of a non-zero digit of \(U\) that is greater than or equal to \(bpow\) across the \(Q\) elements of each thread, and the atomicMin operation uses hardware-supported atomics to efficiently extend this computation of \(n\) across all threads. The second loop finalizes the implementation by subtracting one from each digit whose index is between \(bpow\) and \(n\). This specialization significantly reduces the inter-thread communication required by the general implementation of subtraction, presented in the next section 3.3.

[lst:less-than] implements the general case of the less-than operator \(U < V\) on arbitrary multi-precision integers \(U\) and \(V\). Essentially, the loop computes the less than and equal to relations among the \(Q\) digits of each thread. This is encoded in an int rather than a tuple of booleans: the first (least-significant) bit being set denotes less than (U[i] < V[i]), and the second bit being set denotes equality (U[i]==V[i]). The computation is extended across digits by means of the LTop operator (line [line:lt-apply-seq]) whose semantics on tuple-of-boolean arguments is:

def LTop (lt1, eq1) (lt2, eq2) = (lt2 || (eq2 && lt1), eq1 && eq2)

where lt1/2 denote whether the first/second digits are in a less-than relation and eq1/2 similarly treats equality. The operator implements comparison across two digits: The less-than relation holds either when the second digits are in a less-than relation (lt2) or when the second digits are equal and the first digits conform with the less-than relation. Equality holds if both digits are equal.

Finally, since LTop is associative, the computation is extended across all digits by performing a reduce with the LTop operator across the partial result of all threads. This is implemented as a prefix-sum (scanBlk at line [line:lt-scan]) followed by selecting the last element of the shared-memory result at line [line:lt-reduce]. Next section presents the cuda implementation of prefix sum and the manner in which it is applied to implement subtraction.

3.3 Subtraction of Multi-Precision Integers↩︎

We start by recounting the semantics of the map and exclusive scan (prefix sum) parallel skeletons. Scans are a classical primitive for parallel prefix computations [9]: map applies its function argument to each corresponding element of its input arrays and exclusive scan computes all prefixes under an associative operator \(\odot\) having neutral element \(e_\odot\):

\[\mathrm{map} \;f \;[a_0,\ldots,a_{m-1}] \;= \;[f~ a_0,\ldots, f ~a_{m-1}]\] \[\mathrm{map2} \;f \;[a_0,\ldots,a_{m-1}] \;[b_0,\ldots,b_{m-1}] \;= \;[f~ a_0~b_0,\ldots, f ~a_{m-1}~b_{m-1}]\] \[\mathrm{scan}^{\text{\textit{exc}}} \;\odot \;e_\odot \;[a_0,\ldots,a_{m-1}] \;\equiv \;[e_\odot, \;a_0, \;a_0 \odot a_1, \;\ldots, \;a_0 \odot \ldots \odot a_{m-2}]\]

Subtraction follows a similar procedure as the one used for addition in [2]. Assuming two multi-precision unsigned integers \(x > y\), each having \(m\) digits: \(x = x_0\cdot B^0 + \ldots + x_{m-1}\cdot B^{m-1}\) and \(y = y_0\cdot B^{0} + \ldots + y_{m-1}\cdot B^{m-1}\), their subtraction \(x - y\) is computed by means of a map-scan-map composition:

def $\ominus_{1}$ a b = (a - b > a, a - b == 0)
def $\odot$ (uf1, z1) (uf2, z2) = ( (uf1 && z2) || uf2, z1 && z2 ) 
def $\ominus_{3}$ a b (c, _) = a - b - c
def sub x y = map2 $\ominus_{1}$ x y |> scan$^{exc}$ $\odot$ (false,true) |> map3 $\ominus_{3}$ x y

where the |> operator pipes the result of the preceding computation as the last argument of the following computation. Essentially, \(\texttt{\bfseries{map2}}~ \ominus_{1} ~x ~y\) computes whether the (independent) per-digit subtraction results in underflow (a-b > a) or in 0. The result is passed to \(\texttt{\bfseries{scan}}^{exc} ~\odot~ (\texttt{\bfseries{false}},\texttt{\bfseries{true}})\), which computes and propagates the carry for each digit, and finally, the last map performs the per-digit subtraction to which it applies the carry. The mapped operators \(\ominus_{1}\) and \(\ominus_3\) differ from the ones used for addition, but the combine operator \(\odot\) is the same.

Listings [lst:scan] and [lst:sub] show the cuda implementation: [lst:scan] shows the implementation of the inclusive-scan operator at warp and block level, respectively, for a generic associative operator that accepts arguments of type . For efficiency, the warp level uses register shuffling to avoid accessing shared memory.

[lst:sub] shows the implementation for subtracting as - bs, where inputs as and bs and the result rs are maintained in register memory and as is assumed larger than bs. The combine operator \(\odot\), represented by class CarryOP, is optimized as in [2] to encode the two boolean values as the least-significant two bits of a \(32\)-bit integer. Each thread performs the scan sequentially across its \(Q\) elements (the loop at lines [line:sub-loop1-beg]-[line:sub-loop1-end]) and publishes the final carry in shared memory (line [line:sub-thd-res]). A block level scan (scanBlk) computes and propagates the carries (line [line:sub-scan]), and the subtraction result is updated accordingly at lines [line:sub-update-beg]-[line:sub-update-end]. For very small numbers of digits, hardware-inspired carry-lookahead or carry-select schemes may have lower constant overhead, but our target range is larger multi-precision integers fitting within a CUDA block; accordingly, the implementation uses sequential propagation within each thread and a block-level scan only across the per-thread summaries.

3.4 Implementation of Multi-Precision Multiplication↩︎

The shift-inverse Algorithm 1 uses operations of different (increasing) sizes inside the loop in the Refine function. The operations that have cheap (linear) cost—e.g., comparison, subtraction—are implemented for simplicity by padding to the multi-precision size of the input. However, the multiplications have quadratic cost and need to be adjusted to the actual size of the arguments in order to allow the whole computation to have the same asymptotic as one full multiplication.

Our implementation of size-variant multiplication builds on the one proposed in [2] that addresses the fixed multi-precision case. 6 visually depicts the scheduling proposed in [2] that minimizes inter-thread communication and ensures a load-balanced execution. Denoting by \(m\) the multi-precision size and by \(q\) the sequentialization factor and assuming \(m\) and \(q\) powers of two, the schedule assigns to each thread the computation of \(q\) digits from the first half of the result and another \(q\) symmetrical-opposite digits from the second half. It follows that each thread performs a fixed number \(\frac{Q}{2} \cdot m\) of base (digit) multiplications. 6 shows the cases \(q=2\) and \(q=4\) on the left- and right-hand sides, respectively.

Our approach dispatches at runtime variant-sized multiplications to a number of statically-specialized instantiations of the multi-precision size:

{#lst:dynamic_mult .c++ language="c++" style="cuda" numbers="none" label="lst:dynamic_mult"} if (m <= BLOCK) { smallMul<uint>(m, Ash, Bsh, Csh); } else if (m <= 2 * BLOCK) { effMul<uint, BLOCK, 2>(Ash, Bsh, Csh); } else if (m <= 4 * BLOCK) { effMul<uint, BLOCK, 4>(Ash, Bsh, Csh); } else if (m <= 8 * BLOCK) ...

The pseudocode above assumes that input Ash and Bsh and the result Csh are allocated in shared memory (where the result Csh overlaps with an input) and are relocated to registers inside the specializations. A small multiplication smallMul refers to the case when its precision is smaller than or equal to the cuda block size BLOCK. This is essentially computed by allocating each digit of the result to a different thread or possibly by a procedure that maximizes thread utilization by using atomic accumulations in shared memory. The other cases call effMul that uses the efficient scheduling depicted in 6: the second and third template arguments denote the cuda block size and the sequentialization factor \(q\); it follows that the result is computed in precision \(q\cdot \texttt{BLOCK} \geq m\). In practice we use \(q=4\) as the maximal sequentialization factor.

We recognize that a polished library would not use just one multiplication method throughout — that at different sizes one would transition between classical multiplication, Karatsuba or Toom-Cook multiplication and a number theoretic transform (NTT/FFT) method. For this work, however, we are content to study division in terms of whatever multiplication is used.

4 Empirical Evaluation↩︎

4.0.0.1 Hardware.

The evaluation was run under the Red Hat Enterprise Linux 8.10 operating system on an nvidia A100 GPU, which has 6912 cores, a peak global memory bandwidth of 1,555 GB/sec and FP32 peak performance of 19.5 TFLOP/s.

4.0.0.2 Methodology.

The benchmarking setup is available at:

https://github.com/aske0778/midint-arithmetic-division

and extends the framework presented in [2] with the implementation of division \(\frac{u}{v}\). Each problem instance is initialized with randomly generated integers. The precision of \(u\) is fixed to \(M-2\), such as to account for the two guard digits in the Refine function of Algorithm 1. The precision of the divisor \(v\) is randomly selected between \(2\) and \(\frac{M}{2}\). This configuration ensures that the refinement loop (in Refine) always performs the maximum number of iterations.

The evaluation instantiates the digit type (\(uint\)) to \(64\)-bit unsigned integer and the (maximal) sequentialization factor to \(Q=4\). The results are averaged across \(25\) runs to minimize statistical variance. We compare the performance of our implementation with that of the cgbn library [1] (namely cgbn_div_rem), which is also averaged over \(25\) runs. cgbn does not allow customization of the digits’ word size or of the sequentialization factor.

4.0.0.3 Datasets.

The datasets are chosen to vary the integer precision in number of bits (Num Bits) from \(2^{13}\) to \(2^{18}\) and the number of (parallel) division instances (Num Insts) from \(2^{19}\) down to \(2^{14}\), such that \(\texttt{Num Bits} \cdot \texttt{Num Insts} = 2^{32}\).

4.0.0.4 Validation.

The results of our cuda implementation match those of the GNU Multi-Precision Library (gmp) on all considered datasets. Beyond the performance datasets, correctness testing used \(1000\) randomized instances with random values and effective precisions, plus manual tests of known edge cases, including the shifted-inverse over-approximation that motivated Theorem 2. The tests also covered all multiplication sizes. All results agreed with gmp, providing strong evidence for correctness.

4.0.0.5 Performance Comparison with cgbn.

1 compares our multiplication and division with cgbn. For multiplication, our implementation is faster on the two highest-precision datasets by factors of \(27\times\) and \(3.7\times\), roughly breaks even on the next two datasets, and is slower on the last two by factors of \(1.5\times\) and \(2.1\times\). This is consistent with the different design points: cgbn uses a warp-level organization, whereas our implementation targets larger precisions with one arithmetic instance per cuda block. A precise attribution of the gap to register pressure, scheduling, memory transactions, or occupancy would require a separate counter-level study, which is outside the scope of this paper.

The fifth column shows the factor by which our implementation of division is slower than one full multiplication that is used in its implementation. The highest-precision dataset (\(2^{18}\)) exhibits close to optimal performance for division, being only \(5.2\times\) slower than a full multiplication—recall that Section 2.3 argued that the target division algorithm requires at least \(5\) full multiplications.

The next two datasets slightly enlarge the gap to \(5.4\times\) and \(5.7\times\) slower than one full multiplication, which is still close to optimal. The last three datasets exhibit significant loss of performance as high as \(2\times\) away from the optimal, i.e., from \(7.8\times\) to \(10\times\) slower than one full multiplication. The performance loss is due to:

• assigning a suboptimal number of threads per cuda block—the current implementation computes one division instance with a cuda block of threads, which, for example, has suboptimal size \(32\) for the last dataset.

• the fact that the weight of the five full multiplications in the total runtime decreases with the integer precision, i.e., the other operations become significant.

The sixth column compares the performance of our division with that of cgbn. Importantly, cgbn does not support division at precisions higher than \(2^{15}\) bits, although multiplication is supported suboptimally up to \(2^{18}\) bits. On the last three lower-precision datasets, cgbn division is faster by significant factors ranging from \(3.6\times\) to \(7.4\times\). In particular cgbn’s division is only around \(2.5\times\) slower than one cgbn full multiplication. This indicates that, for the lower precision range up to \(2^{15}\) bits, the cgbn implementation is better matched to the problem size than the block-level division method studied in this paper,¹ while the block-level method covers precisions higher than \(2^{15}\) and up to at least \(2^{18}\) bits. The last column shows the speedup of our division vs. GMP.

The present comparison is therefore intended to test the multiplication-based cost model and to identify the precision range where the shifted-inverse method is promising, rather than to provide a complete production-library benchmark. Further engineering work could combine a counter-level analysis with refinements such as clipped products [10], which could reduce the cost of the double-sized multiplication used in line [line:quot-init] of Algorithm 3.

5 Related Work↩︎

This paper reports a cuda implementation of the multi-precision division algorithm proposed in [3], which extends the Newton method to be carried out entirely in the integral domain by means of data-parallel operations such as shifting and multiplying multi-precision integers. This algorithm can in principle be further improved with the refinement of clipped products [10], which would enable cheaper computation of the middle digits of a multiplication result.

The implementation presented in this paper uses the cuda implementations for addition and multiplication reported in [2], which also inspires the parallelization strategy used for division. The work also reports matching implementations in the high-level data-parallel Futhark language. The Futhark compiler provides a set of useful optimizations, e.g., related to autotuning the degree of exploited parallelism [11], mapping computations to cuda block level with allocation of intermediate results in shared memory [12], reusing memory buffers [13], automatic differentiation [14], [15], static [16] or dynamic [17] verification of memory safety. These have enabled efficient acceleration of applications from various domains [18]–[20], which were exported as python libraries [21]. However, Futhark’s multiplication suffers significant loss of performance because Futhark lacks code transformations aimed at mapping parallel arrays in register rather than shared memory. We have implemented the target division algorithm in Futhark as well; its performance is poor because the mentioned shortcomings are exacerbated by the nature of the division algorithm that requires parallel operations whose sizes vary through a loop, e.g., shifts and multiplications.

The closest related work to the one in this paper is the “Cooperative Groups Big Numbers” (cgbn) library [1], authored by NvidiaLab, that offers high-performance implementation for multi-precision integers up to \(32\)K bits, including multiplication and division. The key technique enabling high performance in cgbn is to map an instance of integer computation on at most one warp of threads in order to leverage specialized nvidia hardware that allows values to be communicated directly between registers (of the warp), i.e., avoiding the overhead of accessing the shared memory, which has significantly higher latency. However, cgbn does not support division on precisions higher than \(32\)K bits. Other less related works targeting parallel multi-precision arithmetic include:

• Previous work has investigated the feasibility of retargeting the multi-precision algorithms used by gmp for GPU execution. The conclusion has been that the architectural differences do not allow easy and/or efficient porting, which motivated development of new, inherently parallel algorithms [22].

• Efficient parallel algorithms for division have been devised to cover specialized input [23]. For example, the case where the divisor is a digit has been efficiently solved with a modified version of Jebelean’s exact division algorithm, which has parallel complexity \(O(n/p + \log p)\), where \(n\) denotes the precision and \(p\) the number of processors.

• CAMPARY [24], [25] supports sequential CPU and GPU execution of multi-precision floating-point arithmetic up to a few hundred bits, e.g., addition, multiplication, division, square root. The key idea is to represent real numbers as unevaluated sums of multiple machine-precision floating point values.

• Older libraries aimed at supporting the functionality of gmp in cuda include GARPREC [26] and CUMP [27], [28]. The latter offers precision up to about \(62\) decimal digits and its objective was to outperform GARPREC. Similar works were ancestors of cgbn, with which we compare directly.

• Isupov [29] uses interval techniques to augment residue number arithmetic for operations that rely on magnitude for numbers with up to 4096 bits.

• A rich body of work [30]–[33] was aimed at accelerating multi-precision integer multiplication on GPUs using adaptations of the classical quadratic algorithm and of Strassen’s log-linear algorithm [4].

Libraries supporting multi-precision arithmetic expose highly-parameterized components (tower of generic types). For example, the division algorithm is parameterized over the underlying multiplication algorithm—e.g., classical, Karatsuba or Strassen—which in turn imposes different representations for the multi-precision integers. Interoperating such libraries with various DSLs and mainstream computer-algebra systems [34], [35], which support different mechanisms for parametrization, remains a direction of interest [36]–[38].

We note that algorithms for multi-precision integer arithmetic are closely related to algorithms for univariate polynomials with coefficients in \(\mathbb{Z} \mod m\). Of particular interest are those with \(m\) a word-sized prime, for which there are good GPU algorithms. An important complication with multi-precision integers is the need to handle carries and borrows.

6 Future Work and Conclusions↩︎

This paper has examined multi-precision integer division on GPUs in a midsize precision range: larger than the range where warp-level cooperative libraries are most effective, but still small enough that a complete division instance can be assigned to a single cuda block. In this regime, the whole-shifted-inverse algorithm of Watt [3] is a good structural match for GPU execution. Its Newton refinement is expressed using integer multiplication, shifting, comparison, subtraction, and related data-parallel operations, avoiding the need to move into a floating-point reciprocal domain.

The implementation maps one division instance to a cuda block, keeps active operands in registers where possible, uses shared memory for communication and staging, and dispatches variable-size multiplications to statically specialized implementations. These choices reduce global memory traffic and expose the compiler transformations needed for a high-level data-parallel language such as Futhark to generate similar code.

The empirical results support the cost model developed in Section 2.3. For the largest tested precision, division is only slightly more expensive than the predicted lower bound of five full multiplications. This indicates that, at these sizes, the overheads from shifts, comparisons, subtractions, and control structure have been kept small relative to the dominant multiplication costs. At smaller precisions these overheads become more visible, and cgbn remains faster where its division implementation is available. Thus the present implementation should be viewed not as a replacement for warp-level libraries in their strongest range, but as evidence that a block-level strategy can effectively cover larger midsize precisions.

Several directions remain open. The most immediate algorithmic improvement is to use clipped products [10]. In the current implementation, the quotient approximation requires a double-precision multiplication even though only a selected range of product digits is ultimately needed. A clipped-product implementation could reduce this cost and bring the observed runtime closer to the idealized multiplication count.

A second direction is better specialization across operand sizes. The current implementation uses a fixed block-level strategy and a maximal sequentialization factor chosen to work well for the largest tested inputs. This is not optimal for smaller inputs, where the number of active threads per division instance can be too low and the relative cost of non-multiplication operations increases. More aggressive autotuning of block size, sequentialization factor, and multiplication specialization would likely improve performance across the full tested range.

A third direction is integration with high-level GPU programming systems. The cuda implementation identifies several transformations that are important for performance: keeping short arrays in registers rather than shared memory, staging global-memory accesses to preserve coalescing, specializing operations by effective precision, and avoiding unnecessary materialization of intermediate arrays. Division is a useful stress test for such compiler work because it combines scans, shifts, comparisons, variable-size multiplications, and iterative refinement in a single exact arithmetic computation.

Finally, larger precisions raise the question of asymptotically faster multiplication. The shifted-inverse division algorithm is parameterized by multiplication, and can therefore in principle benefit from Karatsuba, FFT, or NTT-based methods. On GPUs, however, these methods introduce their own memory-management constraints. Near the fast-memory limit, CRT or NTT-based products may require several modular products for the largest multiplication, while shorter products arising during the Newton refinement may be able to use fewer moduli or longer transform vectors. This suggests that fast multiplication should not be treated as a black-box replacement only; it should be integrated with the changing precision requirements of the division algorithm.

In summary, the results show that exact multi-precision integer division based on whole shifted inverses can be implemented efficiently on GPUs for midsize operands. The implementation reaches near-model performance at the largest tested precisions, covers sizes beyond those supported by cgbn division in our experiments, and clarifies the low-level operations and compiler transformations needed for future high-level implementations.

References↩︎