1 Introduction↩︎

An algorithm is robust if it does not fail for small perturbations of degenerate inputs [1]. Software can be made more robust by the application of multiprecision arithmetic. The algorithms [2] to extend 32-bit floating-point arithmetic using originated in the late sixties [3] and available software packages are QDlib [4] and CAMPARY [5].

A multiple double is an unevaluated sum of nonoverlapping doubles. The multiple double arithmetic multiplies the accuracy of the results, from 16 decimal places to 32 and 64, respectively for double doubles and quad doubles. However, the cost overhead is significant, and parallel computing allows to compensate for this cost overhead. The overhead factors are summarized in Table ¿tbl:tabOverhead?.

This paper describes some recent additions to PHCpack [6], extending contributions described in [7] and [8]. One application is the computation of power series [9]. As the errors in the leading coefficients propagate to the trailing coefficients, the leading coefficients must be computed at higher accuracy than what can be computed in 64-bit double precision.

The code described in this paper serves as a computational preparation for use of tensor cores on a graphic processing unit, in a stepping stone to implement the Ozaki scheme [10].

2 An Error Free Summation↩︎

Assuming all 64-bit doubles have the same exponent, we work with 52-bit integers (fractions of the doubles). The idea for an error free summation is introduced in Figure 1.

If the number of additions does not exceed some threshold, then we have sufficiently many zero bits left at the end of the numbers to represent the result exactly, without any error.

The idea in Figure 1 will be applied to computing inner products with double double arithmetic. Given are vectors \({\boldsymbol{x}}\) and \({\boldsymbol{y}}\) both of length \(n\), of double double numbers, we compute \(\displaystyle \sum_{k=1}^n x_k \star y_k\), where \(\star\) is the double double multiplication.

The double double \(x_k\) is represented by \((x^{hi}_k, x^{lo}_k)\), where the high double \(x^{hi}_k\) and the low double \(x^{lo}_k\) of \(x_k\) are splitted in quarters:

\[(\overbracket[1.2pt][6pt]{x_{k,0},x_{k,1},x_{k,2},x_{k,3}}^{x^{hi}_k}, \overbracket[1.2pt][6pt]{x_{k,4},x_{k,5},x_{k,6},x_{k,7}}^{x^{lo}_k}).\] After splitting also \(y_k\), we compute in double arithmetic: \[s_0 = \sum_{k=1}^n x_{k,0} y_{k,0}, ~ s_1 = \sum_{k=1}^n x_{k,1} y_{k,0} + x_{k,0} y_{k,1},\] and \[s_i = \sum_{k=1}^n \sum_{j=0}^i x_{k,j} y_{k,i-j},\] for \(i=2,\ldots,7\). Then, add \(s_0 + s_1 + \cdots + s_7\) in double double arithmetic.

To examine the computational efficiency, random 64-bit doubles are generated with a fraction of 52 bits in following pattern:

\[1 \underbrace{b b \cdots b}_{12 bits} 1 \underbrace{b b \cdots b}_{12 bits} 1 \underbrace{b b \cdots b}_{12 bits} 1 \underbrace{b b \cdots b}_{12 bits}, \quad b \in \{0, 1\}.\]

Splitting such double into four leads to doubles with fractions

\[\begin{array}{c} 1 b \cdots b ~ 0 0 \cdots 0 ~ 0 0 \cdots 0 ~ 0 0 \cdots 0, \\ 0 0 \cdots 0 ~ 1 b \cdots b ~ 0 0 \cdots 0 ~ 0 0 \cdots 0, \\ 0 0 \cdots 0 ~ 0 0 \cdots 0 ~ 1 b \cdots b ~ 0 0 \cdots 0, \\ 0 0 \cdots 0 ~ 0 0 \cdots 0 ~ 0 0 \cdots 0 ~ 1 b \cdots b. \end{array}\]

By virtue of the placement of the ones in the random fractions, all quarters have fixed exponents, e.g.: 0, \(-13\), \(-26\), \(-39\). All doubles in a multiple double are generated according this pattern.

3 Computational Results↩︎

The results on computing 1,024 times \(\displaystyle \sum_{k=1}^{6144} a_k \star b_k\) in increasing precision is shown in Table ¿tbl:tabResults?.

Times in Table ¿tbl:tabResults? were obtained on an Intel Xeon 5318Y Ice Lake-SP, up to 3.40GHz, 256GB of internal memory at 3200MHz, GNU/Linux, Microway 2024, compiled with GNAT 12.2.0, flags -O3 -gnatp -gnatf.

It takes 9 seconds for 1,024 inner products in hexa double precision (16d). The wall clock time is 9s 308ms, with 85ms for generating the vectors.

In a high level multithread computation, every thread does one inner product. On two 24-core Intel Xeon 5318Y Ice Lake-SP, up to 3.40GHz, 256GB of internal memory at 3200MHz, GNU/Linux, Microway 2024, compiled with GNAT 12.2.0, flags -O3 -gnatp -gnatf, the wall clock time drops to 293 milliseconds, using 96 threads.

Comparing the 293 milliseconds to the 318 milliseconds with one thread in quad double precision, we can quadruple the precision and compute as fast as in quad double precision, using 96 threads, achieving quality up.

4 Conclusions↩︎

Postponing renormalizations of multiple doubles benefits the efficiency.

The convolutions \(\displaystyle \sum_{k=1}^n \sum_{j=0}^i x_{k,j} y_{k,i-j}\) allow to rewrite the inner products in multiple double arithmetic as matrix multiplications in double precision floating-point arithmetic, to prepare for better acceleration with graphics processing units, in particular tensor cores.

References↩︎