Fast Addition of Large Integers

Barry Fagin

Thayer School of Engineering

Dartmouth College

Hanover, NH 03755

barry.fagin@dartmouth.edu

ABSTRACT

Carry-lookahead techniques for binary addition, formerly restricted to small word sizes due to hardware limitations, can now be extended to much larger inputs using massively parallel computers. We present three massively parallel algorithms for the addition of large integers, and analyze their performance.

1.0 Introduction

Unfortunately, carry lookahead circuits implemented in hardware suffer from fan-in and fan-out limitations. These limitations depend on the logic family employed and the size of the input, but in general restrict n to be quite small. TTL circuits, for example, are limited to n=12 [Rh84]. The reader is referred to [Hw79] and [Rh84] for a more detailed discussion of the effects of fan-in and fan-out restrictions on carry lookahead adders. For our purposes, it suffices to note that carry lookahead techniques have been restricted to no more than a few bits of input. This means that large integers must be added in linear time as a series of w-bit adds, where w is the word size of a given machine.

We believe that the recent development of a new computer architecture has made these restrictions obsolete. By permitting the efficient implementation of parallel prefix techniques for large input sizes, existing fan-in and fan-out restrictions no longer apply. This architecture, the massively parallel processor, has thousands of processors, each connected to a local memory and a communication network. By distributing the integers to be added among the processors, parallel prefix techniques can be employed to rapidly add large numbers faster than if conventional machines were employed.

We first discuss the basic computational model of a massively parallel processor. We then present the basic addition algorithm, in which each processor contains one bit of each of the two the integers to be added. (We assume, for now, that all integers are unsigned). Our work extends that of [Bu46] to non-Von Neumann architectures, in which multiple groups of bits may be added simultaneously, followed by parallel carry propagation. We show how performance can be improved by exploiting the average case behavior of large n-bit additions and the asymmetry of the computation time of

two particular operations. We then show even better performance by grouping multiple bits per processor. Finally, performance measurements of all these algorithms are presented and discussed.

2.0 The Computational Model

Figure 1

3.0 The Parallel Prefix Problem

and an associative operation ø, find all partial products

4.0 Addition: The Basic Algorithm

We present here our basic algorithm for the parallel addition of large integers on massively parallel processor. It differs from that of [Pe88] in that it is an iterative algorithm implemented on an existing multiprocessor, requiring log n + 1 communication steps. The algorithm described in [Pe88] is recursive, PRAM-based, and requires 2 log n steps.

The addition problem is a special case of the parallel prefix problem [La80]. Let A and B be two n-bit integers, and let A[i] and B[i] denote bits i of A and B. The carry out bits may be obtained by calculating all partial products

and ø is the following two-input, two-output Boolean function:

ø(xj,xk)[1] = xj[1] • xk[1]

ø(xj,xk)[0] = xj[0] | (xj[1] • xk[0])

The function implemented by ø is easier to see if depicted graphically, as shown below:

Figure 2

When the parallel prefix calculation is complete, the carry out bits are the xi[0]'s. Bit i of the sum, S[i], is calculated by shifting the carry bits over one position and performing a three way XOR on the carry, A[i], and B[i] bits:

The logic equations for the calculation of parallel prefix products may be unfolded and implemented in hardware using Boolean logic elements, but fan-in and

fan-out restrictions limit their application to values of n < 16. If, however, we employ a massively parallel processor, then we may use the algorithm described previously to calculate the products in O(log n) time for any n up to the number of available processors. The simplest way to implement the parallel prefix calculation for addition is to have each processor maintain two Boolean variables p and g. These are initialized to

(3) g <- A[i] • B[i]

These are simply the two bits of the input xi.

On communication step i, numbering from zero, all processors with neighbors 2i to the right perform the operation

(4) p <- p • neighbor[p]

(5) g <- g | (p • neighbor[g])

where neighbor[var] refers to the variable of the processor with which communication is occurring. Equations (4) and (5) are implemented in parallel; the p term on the right hand side of (5) refers to the value of p before it is updated by (4).

These equations implement the combining function ø. After log n communication steps, the carry outs are stored in the g's. These can then be shifted to the left to calculate the sum.

We refer to this algorithm for n-bit addition as Algorithm I.

4.1 An Improved Algorithm

We may improve Algorithm I by examining equations (4) and (5) more closely. We see that once a processor's p value is set to zero, no further computation by that processor is necessary. A p value, once zero, will remain zero, while a g value will remain unchanged. Thus the computation may terminate when all p values are zero. (We may ignore the rightmost p, as its value will not affect the computation). We detect termination by performing a logical OR of all p values. If the communication network is a Boolean n-cube, this can be done in O(log n) time.

It may seem that further improvements can be obtained by shutting down processors when g=1, regardless of the value of p, since once a carry out has been generated from a bit position no further work is necessary. In fact, the case of g=1 is subsumed by the case of p=0. A simple inductive argument will show that at the end of 0 or more communication steps, there are no processors for which g = 1 and p = 1.

4.2 Average Case Analysis

Proof: The proof proceeds by induction on i. If the input integers A and B are randomly chosen, any bit in their bitwise XOR will be 1 with probability 1/2. Thus the theorem is true for i=0.

The following theorem follows almost immediately:

Proof: For 0 <= j < n, let pj denote the value of p on processor Pj after i <= log n communication steps. Letting E denote expected value, the average number of active processors is

j=n-1 j=n-1

j=0 j=0

by Theorem 1. QED Theorem 2.

Since the average number of active processors is < 1 after | log log n | + 1 steps, we may improve the performance of Algorithm I by checking if all p values are 0 after step | log log n |, terminating once this occurs. We refer to this algorithm as Algorithm II.

4.3 A More Precise Analysis

To obtain a more detailed analysis of the number of communication steps required, we proceed as follows. We first define some terms. We modify our definition of an active processor, defining a processor to be active after step i if it has p=1 and a neighbor 2i to the right. We define a block of size m to be a group of m consecutive 1's in an n-bit input. We define a proper block to be a block with no 1's in the least significant bit. We say that a processor heads a block of size m if its initial input bit is 1 and the m-1 processors on its right, if they exist, have initial inputs of 1.

We will need the following theorem:

A processor is active after i communication steps if and only if it heads a proper block of size >= 2i.

We may now prove the following theorem:

Proof: Clearly the algorithm requires no communication steps if the size of the largest proper block is zero, as all p values will be zero before communication begins. We prove the non-trivial case below.

Note that the worst case occurs when the input is all 1's. Since m cannot exceed n-1, an input of all 1's requires log n communication steps. Notice as well that the addition algorithm as we have presented it requires one additional communication step to shift the g's after the p's are driven to zero.

Thus we have

where S(n) equals the average number of communication steps required to add two n-bit numbers, and D(i,n) denotes the number of n-bit digits that will cause Algorithm II to terminate in exactly i steps. The +1 term refers to the left shift of the carries at the end of the addition algorithm.

(8) D(0,n) = 2

noting that N(j,n-1)*2 is just the number of n-bit integers with a proper block of size j.

It remains, then, to calculate N(m,k), the number of k-bit integers with at least m consecutive 1's. Clearly if k < m or k = 0 then N(m,k) = 0, and if k=m and k > 0 then N(m,k) = 1. Suppose we know N(m,k-1). How do we calculate N(m,k)? To obtain all k-bit sequences with at least m consecutive 1's, we may first take all k-1-bit sequences with at least m consecutive 1's and add either a one or a zero. Thus

(10) N(m,k) >= 2N(m,k-1)

(11) N(m,k) = 0 k < m or k = 0

(12) N(m,k) = 1 k = m and k > 0

No closed form solution of (11)-(13) is known. However, [Us37] gives the following approximation for N(m,k):

where

_

i=2

Equations (14) through (17) give a way to estimate N(m,k), and hence S(n), the expected number of communication steps required to add two n-bit numbers using parallel prefix techniques. A proof that S(n) < log log n + 2 for n >= 4 is provided in [Fa91].

Some values of T(n) are shown in Table 1:

log n T(n)

1 1

2 1

3 1

4 1.0175

5 1.0001

6 1.0003

7 1.0008

8 1.0018

We may now estimate the average case running time of the modified addition algorithm. Let c0 denote the interprocessor communication time, and c1log n the time required to perform a global OR. The average number of communication steps is S(n), done in constant time. We conjecture that the average number of termination steps is

5.0 The Effect of Multiple Bits Per Processor

This method differs from the 1-bit per processor algorithm only in the initial and final stages. In addition to each processor keeping a p and g bit, it also maintains a d-bit digit D. Let A and B be the two n-bit integers to be summed, and let Ai and Bi denote the d-bit digit of A and B assigned to processor i. Let Di[j] denote bit j of the digit maintained by processor i. The initial equations (2) and (3) become

(19) Di <- Ai plus Bi

(20) p <- logical AND of all bits of Di

(21) g <- carry out of (19)

Each processor adds its own two digits locally, storing g and p as the carry out of the sum and the logical AND of all the sum bits, respectively. The p's and g's are then updated using the parallel prefix algorithm. The carries are shifted to the left, as in the last stage of Algorithms I and II, but are then added to each processor's local sum to compute the final result.

To solve this problem, we recalculate the recurrence in (11)-(13) to include a third variable b, the total of number of 1 bits. Let N(m,k,b) denote the number of k-bit integers with exactly b 1 bits and at least m consecutive 1 bits. Clearly, we must have N(m,k,b) = 0 unless k >= b >= m. Likewise, N(m,k,b) = 0 if k = b = m = 0, or if k > 0, b = 1, and m = 0, and N(m,k,b) = 1 for m = k = b. Let us call a k-bit integer with m consecutive and b total 1 bits an (m,k,b) integer Following the original derivation of N(m,k), we may obtain an (m,k,b) integer by a) appending a 1 to an (m,k,b-1) integer, b) appending a 0 to an (m,k-1,b) integer, or c) appending a 1 to a (m-1,k-1,b-1) integer that ends in m-1 consecutive 1's that does not already contain m consecutive 1's. This gives:

(22) N(m,k,b) = 0 unless (k >= b >= m)

N(m,k,b) = 0 k = b = m = 0

N(m,k,b) = 1 k > 0, b = 0, m = 0

N(m,k,b) = 0 k > 0, b = 1, m = 0

N(m,k,b) = 1 k = b = m, k > 0

N(m,k,b) = N(m,k-1,b-1) + N(m,k-1,b) + bi(k-m-1,b-m)

- N(m,k-m-1,b-m) k >= b >= m, k > m, (b > 1 or m = 1)

where bi(x,y) denotes the binomial coefficient of x and y.

This gives D(i,n,d), the number of input pairs of n*d-bit integers that cause Algorithm II to terminate in i steps, as

D(i,n,d) = 0 < i <= log n

D(i,n,d) = 0 i > log n

We may now calculate the expected number of communication steps as

i=log n

i=0

Since the probability of any input p-value being 1 declines exponentially with d, we should expect Sd(n) to be significantly smaller than S1(n) for d > 1, with Td(n) drawing closer to 1 (recall that Td(n) is the expected number of termination checks). Values of Sd(n) for small n are shown in Table 2. Values less than 1.01 are indicated by a '-'.

Table 2

log n log d=0 log d=1 log d=2 log d=3 log d=4 log d=5

1 1.5 1.25 1.06 - - -

2 2.25 1.69 1.18 1.01 - -

3 2.88 2.16 1.39 1.03 - -

4 3.34 2.55 1.67 1.06 - -

The average case running time of the algorithm is

As d increases, the middle two terms decrease exponentially while the first term increases linearly. The optimal value of d is machine dependent, occurring when the time saved by using large d to reduce the probability of interprocessor communication matches the time required to add d-bit digits.

We refer to Algorithm II modified to handle multiple bits per processor as Algorithm III.

6.0 Performance Results

Algorithms I,II 3 + 2 log n

Algorithm III 3d + log n + 132

n(3d + log n + 132) <= 231

The ratio of virtual processors to physical processors is the VP ratio. When the VP ratio is greater than 1, the system is saturated. For saturated systems, the

execution time must grow linearly, as the algorithms described previously require as

many processors as inputs for log n or log log n growth.

For more information on the Connection Machine, the reader is referred to [TMC89a] and [TMC89b].

6.1 Addition on a VAX

For comparison purposes, we present the time required for a VAX 11/785 to add large integers using the following code:

LABEL:

adwc (r0)+ , (r1)+

acbl $DIGITS-1, $1, r2, LABEL

where $DIGITS represents the number of 32-bit digits to be added. "adwc" is the VAX add-word-with-carry instruction, while "acbl" stands for "add, compare, and branch longword". We assume the two integers are stored in contiguous regions of memory, with the least significant digits addressed by r0 and r1. The sum is stored in the region addressed by r1, while r2 contains the number of digits.

The performance of Algorithms I, II, and III as compared with the VAX is shown below.

Figure 3

Each data point for the Connection Machine represents the average addition time of between 100 and 10000 pairs of integers.

We see that the massively parallel algorithms can grow logarithmically only if sufficient processors are available. Once the system is saturated, the running time becomes linear.

For small values of n, the advantages of carry lookahead are not enough to overcome the overhead of interprocessor communication and the fast addition time of the VAX. The breakeven point occurs near saturation, at which point all algorithms grow linearly. For n = 32, the VAX is about 500 times faster. For n = 32M, Algorithm I on the CM is about 4 times faster.

As expected, the running time of Algorithm II grows as log log n, until n equals the number of processors. The breakeven point occurs at approximately log n =13. For n = 32, the VAX is about 400 times faster. For n = 32M, Algorithm II on the CM is about 25 times faster.

Figure 4 shows the effect of digit size on execution time for Algorithm III.

Figure 4

We observe that for large d, the time required to add s=n*d bit numbers is independent of n. This is because for non-saturated systems with large d, no carries are generated, and the algorithm is equivalent to computing n d-bit sums in parallel.

7.0 Conclusions

These algorithms are made effective through the utilization of a massively parallel processor, a computer architecture that does not suffer from the fan-in/fan-out restrictions that have previously limited the size of fast adders. We have presented performance results of all three algorithms on a 32K Connection Machine, a massively parallel processor developed by Thinking Machines Incorporated. These results were compared with a VAX 11/785 running an extended precision addition routine.

For small n, the VAX remains superior. We suspect this will hold for any conventional uniprocessor architecture. At small values of n, the advantages of extended carry lookahead are outweighed by the overhead of interprocessor communication and bit-serial arithmetic. As n increases, however, the performance advantages of a log n algorithm begin to approach those of a linear one, until a breakeven point is reached. For the Connection Machine and the VAX 11/785, the breakeven points were at log n = 15, 13, and 12 bits for Algorithms I, II, and III.

We note that the superiority of linear algorithms for small integer sizes does not render the CM ineffective. A P-processor, 150-ns CM can perform P local N-bit additions in .15N microseconds, subject to memory and microcode limitations. For small values of N this is significantly faster than Algorithm III, and is much more competitive with uniprocessor architectures. Since a CM can perform P additions in the same amount of time that it takes to perform one, applications which require large numbers of small additions may run faster on the CM.

After the breakeven point, the lookahead algorithms on the CM significantly outperform the VAX. For the largest values of n, Algorithms I, II, and III run 4, 25, and 1200 times faster. These figures are generous to the VAX and other minicomputers in that their execution times ignore memory management limitations. It is not possible, for example, for an ordinary user to access 512Mbits of memory on the VAX. We note as well that execution times on both machines could be significantly improved through the use of special-purpose microcode.

Algorithms II and III remain effective even when cost is considered. A 32K CM costs an order of magnitude more than a VAX, while algorithms III and III provide 25 and 1200 times more performance. Thus these algorithms offer both improved performance and improved cost/performance. We note, however, that the CM is a highly specialized machine, with an architecture well suited to the problem under study. Other problems that do not map as closely to the CM will exhibit poorer performance.

We have limited our discussion to unsigned integer addition. If 2's complement representation is used, subtraction can be performed by permitting a carry-in C0 to the least significant processor P0. Let g0 be the g variable of processor 0, A0 and B0 the least significant input bits. If g0 is initialized to A0•B0 OR A0•C0 OR B0•C0, and if the least significant bit of the sum is computed as a three-way XOR of A0, B0, and C0, then A - B can be computed by adding A to the 1's complement of B with C0=1. For d-bit digits, we modify (19) for processor 0 to include C0 in its local sum.

Our analysis of the average-case behavior of Algorithms I-III has assumed random data. For many practical applications, this may not be the case, particularly if 2's complement notation is used For example, small positive numbers contain many leading zeros, while small negative numbers contain many leading 1's. The resulting vector of p-values will contain many consecutive 1's, requiring more communication steps. For these types of inputs, the algorithm will require more steps than predicted by our analysis. Under such circumstances, it may be appropriate to delay the termination check of p-values, or to simply execute log n + 1 communication steps by default.

For future work, we are attempting to implement other arithmetic functions on the CM, including the multiplication of large integers. We are currently investigating the use of Polynomial, Fermat, and Fast Fourier transforms to compute large integer products. The interested reader is referred to [Fa89], [Fa90a], and [Fa90b].

8.0 Acknowledgements

9.0 References

[Ba80] Baer, J. Computer Systems Architecture, Computer Science Press, Potomac Md., 1980.

[Bu46] Burks, A. W. et. al. Preliminary Discussion of a Logical Computing Instrument, Report to the U. S. Army Ordinance Department, 1946.

[Fa89] Fagin, B. Negacyclic Convolution Using Polynomial Transforms on Hypercubes, to appear in IEEE Transactions on Signal Processing .

[Fa90a] Fagin, B. Multiplication of Large Integers on Hypercubes, to appear in the Journal of Parallel and Distributed Computing (Mar 1990).

[Fa90b] Fagin, B. Large Integer Multiplication on Massively Parallel Processors, Frontiers of Massively Parallel Computation (1990), pp 38-42.

[Fa91] Fagin, B. Fast Addition of Large Integers, Technical Report, Thayer School of Engineering, Dartmouth College, Hanover NH 03755.

[Fl75] Floyd, R. The Exact Time Required to Perform Generalized Addition, Sixteenth Annual Symposium on Foundations of Computer Science, (1975), pp 3-5.

[Hi86] Hillis, W. and Steele, G. Data Parallel Algorithms, Communications of the ACM 29,12 (Dec 1986), pp 1170-1183.

[Hw79] Hwang, K. Computer Arithmetic: Principles, Architecture, and Design, Wiley, New York, NY, 1979.

[Hy82] Hyman, A. Charles Babbage: A Pioneer of the Computer. Princeton University Press, Princeton NJ, 1982.

[Kn81] Knuth, D. The Art of Computer Programming, Vol. II, Seminumerical Algorithms, Addison-Wesley, Reading Mass., 1981.

[La80] Ladner, R. and Fischer, M. Parallel Prefix Computation, JACM 27,4 (Oct 1980), pp 831-838.

[Pe88] Peng, S. and Hudson, T. Parallel Algorithms For Multiplying Very Large Integers, 18th International Conference on Parallel Processing (1988), vol. 3 pp 173-177.

[Rh84] Rhyne, T. Limitations on Carry Lookahead Networks, IEEE Transactions on Computers (Apr 84), C-33,4 pp 373-374.

[Sp69] Spira, P. The Time Required for Group Multiplication, JACM 16,2 (Apr 1969), pp 235-243.

[TMC89a] Connection Machine Technical Summary, Thinking Machines Corporation, Cambridge, Mass.,1989.

[TMC89b] C/Paris 5.1 User's Manual, Thinking Machines Corporation, Cambridge, Mass., 1989.

[Us37] Uspensky, J. "Introduction to Mathematical Probability", McGraw-Hill, New York, NY, 1937, pp 77-84.

[Wi65] Winograd, S. On the Time Required to Perform Addition, JACM 12,2 (Apr 1965), pp 277-285.

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

FOOTNOTES:

Affiliation of author: Thayer School of Engineering, Dartmouth College, Hanover NH 03755. e-mail: barry.fagin@dartmouth.edu.

4 We were unable to implement some combinations of n and d that met the constraints. These combinations came very close to the upper memory limit; some portion of memory must be left free for system use.

FIGURE CAPTIONS

Figure 1: A Massively Parallel Processor

Figure 3: Comparative Performance

Figure 4: Effect of Digit Size