FAST MODULAR MULTIPLICATION USING
2-POWER RADIX
Colin D. Walter
Department of Computation
UMIST
PO Box 88
Manchester M60 1QD, UK
www.co.umist.ac.uk

Abstract. Details are given here of how to generalise Brickell’s fast
modular multiplication algorithm to when the number representations
have a general 2-power radix. Correct action depends upon the satisfaction of a complicated inequality and speed upon the use of a redundant
number system to enable parallel digit operations. The effect of varying
the radix on the efficiency of hardware implementations is considered.
Improved efficiency has repercussions in public key cryptography where
the RSA encryption scheme may use this type of algorithm for its modular exponentiations. However, it is shown that there is no advantage intaking a very large radix, although a small increase above 2 is beneficial.
Key words: Modular Multiplication, Fast Computer Arithmetic, Digital Arithmetic, Integer Division, Parallel Bit Operations, Redundant
Number Systems, RSA Algorithm, Cryptography
C.R.Categories F.2.1, B.7.1, E.3.

0Introduction

We investigate here a parallel digit algorithm for the calculation of A×B mod M .
A modular product R of B with part of A is created iteratively by using the
successive digits of A to add a multiple of B and predicting the right multiple of
M to subtract at the next cycle. So the whole process takes a time proportional
to the number of digits in A. There is an apparent advantage in throughput
speed from increasing the radix of the representation in order to decrease thenumber of digits in A. Here we prove a very general algorithm and look closely at
the major terms describing hardware area and clock speed in an implementationin order to clarify the best choice of radix.
The history of the subject includes the use of redundant number systems(e.g. [2]) applied to division ([1]). Arising from this, modular multiplication using

2

International J. of Computer Mathematics, vol. 3, 1991, pp. 21-28

these techniques is discussed by Simmons and Norris [7], and described for base
2 by Brickell [3], with proofs and further detail provided by Gibson [4] andWalter and Eldridge [8]. The recent rash of publications in this area also includes
notable work by Selby and Mitchell [6]. The main application of the results here
is in hardware implementations of the RSA algorithm [5] where the modular
exponentiation of encryption and decryption is performed by repeated modular
multiplications.

1

Notation

Numbers here are represented with 2-power radix r = 2ρ , ρ ≥ 1, making translation to or from a binary form fairly trivial. The digits of a numberPX are written
xi and lie in the range 0..xmax , say, so that X has the value i≥0 xi ri . The
i+1st digit of X is xi , the coefficient of ri , because the numbering starts at 0.
In an irredundant system xmax = r−1, which makes the choice of digits unique,and in a redundant system xmax ≥ r, which allows alternative representations. A
typical choice might be xmax = 2r−1 because full use is then made of the extra
bits needed for a minimal, properly redundant system. Since hardware area iscrucial xmax is kept below 2r so that a digit requires ρ bits if the representation
is not redundant, and ρ+1 if it is redundant.
Digit ranges are chosen to suit applications. Here, we are interested in repeated modular multiplication to perform the modular exponentiation of the
RSA two-key cryptosystem. Thus the modulus M is often fixed for many successive modular multiplications. Since M is usually in binary form, we obtain
M in a non-redundant form to base r by taking its bits in groups of ρ. So we
assume its digits are at most mmax = r−1.
The algorithm presented here calculates a residue R and an integer quotient Q
satisfying A×B = M ×Q+R where R is either the smallest non-negative residueof A×B mod M or differs by at most M from it. This occasional extra M can be
subtracted separately if required. Suppose M has N digits. If the output R has
more than N digits it is larger than M and can be reduced at virtually no cost.
So we may assume R has at most N digits. Since A and B are initially in therange 0..M in the RSA algorithm, and the output is used as subsequent input,
we can always assume A and B have at most N digits. Normally the hardware
fixes N , whilst M is shifted to put it in the right range, resulting in appropriate
shifts for the other quantities initially and finally. So we assume M lies in therange rN /2..rN −1 and also that each of A, B, R, Q and M is non-negative.
The multiplication is performed by repeated addition, but to perform addition in a single clock cycle requires each output digit to be dependent on only asmall number of input digits, not on all lower-indexed digits as in the case of theusual (serial) hand calculation. So carries must propagate at most only a small
fixed distance. This can only be done using a redundant number system for the
output. Since in our application, output from one modular multiplication becomes input for the next, we expect A, B, and any intermediate partial productR to have a common redundant number system, say with digits 0..rmax . Later

C. D. Walter

3

on, sensible values for rmax will be considered. The integer quotient Q also turns
out to need a redundant representation, say with digits in the range 0..qmax .
One of the ideas tacit in Brickell’s paper [3], which we make explicit here, isto calculate A×B mod M as (AS×B mod M S)/S for a fixed r-power S = rE .
So a shift up by E places is performed. Here the number E of extra digits neededat the top of the registers is large enough to make the contribution of B to thepartial product R insignificant in a certain precise way. Henceforth, we therefore
work with RS, M S, etc., which we consider as individual numbers, not products,
with digits denoted rsi , msi , etc.
To predict the right multiple of M S to subtract from intermediate partial
products RS in order to reduce them mod M S we look only at the top few digitsof the numbers involved. The notation is that X denotes the number obtained by
setting the lowest bits of a number X to 0 and EX is the greatest number withthe same representation as X which is zero when is applied. The error involved
in using X as an approximation for X is expressed in X ≤ X ≤ X + EX . If thelowest δ digits and β bits from the next digit are set to zero in X, then EX =xmax (rδ −1)/(r−1) + rδ (2β −1). So, in the non-redundant case, EX = rδ 2β −1.

2Software Algorithm

Next is a pseudo-Pascal version of the modular multiplication algorithm understudy here. Missing detail and restrictions are derived later. In particular, L is
the main constant yet to be determined, suggestions are needed for defining thefunction ApproxQuot, and the register sizes for R and Q are outstanding. It
is an easy step subsequently to obtain (A×B) mod M and (A×B) div M precisely, because we will find that L < 2M , so that the property R < L for the
output residue makes it either the least non-negative one, or the second least.
Const
S = rE ;Type
Mdigits = 0..r-1; Qdigits = 0..qmax ; Rdigits = 0..rmax ;Procedure ModMult(A, B : Array[0..N-1] of Rdigits ;M : Array[0..N-1] of Mdigits ;Var R : Array[0..N+?] of Rdigits ;Var Q : Array[0..N+?] of Qdigits ) ;{ Pre-Condition:rN /2 ≤ M < rN , L = ?
}{ Post-Conditions: A*B = Q*M + R, R ∈ 0..L }
Var

J
MS
AS
RS

: 0..N+E-1 ;: Array[0..N+E-1] of Mdigits ;: Array[0..N+E-1] of Rdigits ;: Array[0..N+E+?] of Rdigits ;Function ApproxQuot(rRS,MS : Array...) : Qdigits ;{ Post-Condition : ApproxQuot ≈ rRS div MS }Begin ... End ;4International J. of Computer Mathematics, vol. 3, 1991, pp. 21-28
Begin { ModMult }MS := Shift(M,E); AS := Shift(A,E); RS := 0;
For J := N+E-1 DownTo 0 do
Begin
{ Invariant: RS ∈ 0..LS }Q[J] := ApproxQuot(rRS,MS) ;RS:= Shift(RS,1) + AS[J]*B - Q[J]*MS ;End ;R := Shift(RS,-E)
End ; { ModMult }Q := 0;Using the subscripts of our notation so far rather than the indexing [ ] ofP
PN +E−1
PN +E−1 i−J
Pascal, if AS = i asi ri then RSJ ≡ ( i=J
asi ri−J )B−( i=J
qi r )M S
satisfies RSN +E = 0 and RSJ ≡ rRSJ+1 + asJ B − qJ M S. Hence RSJ is thevalue of RS at the end of the iteration with J as the value of the control variable. Therefore the final iterate RS = RS0 satisfies RS ≡ AS×B − Q×M S,
and the output satisfies A×B = Q×M + R, with R ≤ L if the claimed loop
invariant holds. When a suitable L is determined, the outstanding register sizes
can be found easily, assuming representations of the numbers which require themaximum number of digits. In particular, we will find that R can be squashedinto N digits.

3

Convergence

If one partial product RS is in the range 0..LS, then we require that the next
partial product, say rRS+aB−qM S, be in the same range so that the process
converges − otherwise the output R may become too large. In this section weestablish the property required of L for this to be the case.
Ideally, the choice of q should make the expression for the next value ofRS minimally non-negative, allowing L to be chosen equal to M . However, thiswould require checking every digit of RS, aB and M S. Time is saved by using
only the topmost bits of these inputs to generate an approximation for q. This
forces a wider range 0..LS for the next RS, and possible extra final subtractionsof M to achieve R in the interval 0..M .
Since S will be chosen large enough to make the term aB relatively small, qwill satisfy rRS−qM S ≈ 0, yielding q ≈ rRS div M S, an approximate quotient.
As the value of q is predicted by looking at only the top few digits of M S and RS,
we define q = ApproxQuot(rRS, M S) for a suitable function. Here ApproxQuot
may be obtained from a pre-computed look-up table, parametrised only by RS,
which is re-set each time M is changed (cf Selby & Mitchell [6]). The table can becomputed easily from the values of qM S as q varies, these values being uniquelydetermined because the representation of M S is not redundant.
Although the choice of ApproxQuot is in the hands of the user, one possibility
is q = rRS div (M S + EM S +1). This is slightly less accurate than the tabulated
values of q suggested above, but easier to prove properties about. The inequalities

C. D. Walter5established from this definition of q will work even better for a more accurate
approximation to rRS div M S. Here EM S +1, rRS and M S are all multiples ofa very large power of 2, so that shifting them down leaves a simple definition ofq. Because M S ≤ M S + EM S +1, we knowrRS + aB − qM S
≥ rRS − {rRS div (M S + EM S +1)}×M S≥ rRS − {rRS div M S}×M S > 0
which is the desired lower bound. Similarly, an upper bound is obtained fromrRS + aB − qM S
= rRS + aB − q(M S+EM S +1) + q(EM S +1)= rRS + aB − {rRS div (M S + EM S +1)}×(M S+EM S +1)+{rRS div (M S + EM S +1)}×(EM S +1)≤ rRS + aB − {(rRS − rERS ) div (M S+EM S +1)}×(M S+EM S +1)+(rRS div M S)×(EM S +1)= X + rERS + aB − (X div Y )×Y + (rRS div M S)×(EM S +1)
where X = rRS − rERS and Y = M S+EM S +1
= X mod Y + rERS + aB + (rRS div M S)×(EM S +1)< Y + rERS + aB + (rRS div M S)×(EM S +1)≤ M S + EM S + 1 + rERS + rmax 2 (rN −1)/(r−1) + (rLS div M S)×(EM S +1).
This last expression is larger than M S, but only by a little if S is large enoughand sufficiently many top bits are used. It is almost the value to choose for LS,
which we need to choose slightly larger.
The other constraint required for the convergence of the algorithm is that
ApproxQuot produces a result in the range 0..qmax . Certainly q ≥ 0 in thiscase. Because here, as in any definition of q, rRS − qM S ≥ 0, we know q ≤rRS div M S ≤ rLS div M S if we pick LS as intended. Hence q ≤ qmax if
rLS div M S ≤ qmax . Now, with this assumption, define
LS = M S + EM S + 1 + rERS + rmax 2 (rN −1)/(r−1) + qmax (EM S +1).
Then the working of the last paragraph proves rRS + aB − qM S < LS, whichis the required upper bound to ensure that the loop invariant is satisfied. So theprocess will converge as long as rLS div M S ≤ qmax .
Inserting the value for LS in this assumption yields the sufficient conditionr{M S + rERS + rmax 2 (rN −1)/(r−1) + (qmax +1)(EM S +1)} div M S ≤ qmax
for everything to work, or, equivalently,
rERS + rmax 2 (rN −1)/(r−1) + (qmax +1)(EM S +1) < (qmax −r+1)M S/r.

6International J. of Computer Mathematics, vol. 3, 1991, pp. 21-28

This is clearly sharpest when M S is at a minimum and illustrates that Q musthave a redundant representation, i.e. there is only a solution if at least qmax ≥ rholds. It is certainly solved by taking S large enough and using sufficiently many
digits to approximate M S and RS for any fixed such qmax . Since the general
upper bound qmax < 2r was an initial hypothesis, rLS div M S ≤ qmax < 2r sothat L < 2M , proving the contention that the output is at most M larger thanthe minimum non-negative residue.

4

A Solution

Assume that the lowest α+E digits and the lowest β bits of the next digit are
set to 0 in M S by , and that σ+E digits and τ bits are annihilated similarlyin the case of RS. Moreover, suppose as usual that M ≥ Mmin = rN /2. The
final inequality of the last section may now be re-written asr(rmax (rσ+E −1)/(r−1) + rσ+E (2τ −1)) + rmax 2 (rN −1)/(r−1) + (qmax +1)rα+E 2β< (qmax −r+1)rN +E /2r
This is what we need to satisfy for the algorithm to behave.
If q is obtained from a look-up table pre-computed from M S and indexed
by RS, then we want to minimise the number of significant bits in RS, that is,
maximise σ(ρ+1) + τ . To do this, take qmax as large as possible, viz 2r−1. Then,from choosing just the two most significant terms on the left in the above, at
least
rσ+E+1 2τ + rα+E+1 2β+1 < rN +E /2
must hold, giving σ = N −2, τ = ρ−2 and α = N −2, β = ρ−4 as the most
optimistic solution possible.
How many bits of RS are needed? Using RS < LS < 2M S < 2rN +E , there
may be 1 significant bit in the N +E+1st digit of RS. Also ρ+1 bits come fromthe N +Eth digit, and 3 from the N +E−1st, since ρ+1 are needed for each digit.
This is a total of ρ+5 significant bits in RS. However, it is easy to OR the only
bit of the N +E+1st digit with the ρ+1st bit of the N +Eth digit to reduce thisby 1 (so the N +E+1st digit of RS is not needed). For MS only ρ+4 bits areneeded if the above solution works because M S has a non-redundant form with
only N +E digits. In this case, however, M S will have been shifted up initiallyso that its most significant bit is the top one of its N +Eth digit. So only ρ+3
bits of M S affect the computation of ApproxQuot. Thus, direct calculation ofthat function requires a total input of 2ρ+8 bits if the above choices satisfy ourinequality.
When does the proposed solution above not work? Substituting the values
into the inequality yieldsr(rmax (rN −2+E −1)/(r−1) + rN −2+E (r/4−1)) + rmax 2 (rN −1)/(r−1) + rN +E /8< rN +E /2

C. D. Walter

7

Suppose we insist, reasonably, that rmax ≤ 2(r−1). Then this holds if
2r(rN −2+E −1) + rN −1+E (r/4 − 1) + 4(rN −1)(r−1) < 3rN +E /8
i.e. ifrN −1+E − 2r + 4(rN −1)(r−1) < rN +E /8
For large enough E, this is therefore satisfied if r > 8, i.e. r ≥ 16, or ρ ≥ 4,
which, of course, is the only range of values for which the choice of β(≥ 0) makes
sense.
The largest value needed for E is when r is least, i.e. r = 16. The lastinequality is satisfied if
rE−1 + 4(r−1) ≤ rE /8
i.e. if 4(r−1) ≤ rE−1 (r/8−1), and so if 4r ≤ rE−1 . Thus E = 3 suffices forr ≥ 16, with E = 2 easily if r ≥ 64.

5

Summary

The general radix version of Brickell’s modular multiplication algorithm given
in section 2 works when L is given by
LS = M S + EM S + 1 + rERS + rmax 2 (rN −1)/(r−1) + qmax (EM S +1)
and qmax = 2r−1, rmax < 2r−1, Mmin = rN /2, M S is obtained by setting thelowest N +E−2 digits and the lowest ρ−4 bits of the next digit to 0, RS is
obtained similarly by setting the lowest N +E−2 digits and next ρ−2 bits to 0,
E = 3 for r ≥ 16, E = 2 for r ≥ 64, and r ≥ 16.
We do not need to work out LS explicitly, but we know it is less than 2M S,
so that the output R is bounded by 2M . The least value Mmin of M can bejustified by shifting M up by a power of 2 to make full use of the N th digit.
Solutions exist to the inequality for the values 2, 4 and 8 of r and can beobtained in exactly the same way as above. This is also true for different choicesof rmax and qmax provided they each allow redundancy.

6

Hardware and Efficiency

In the loop each addition and the generation of the next q can be done together
in a single clock cycle, with the result that N +E cycles are required for thealgorithm to execute. Clock speed is bound by the number of gates on the longest
or critical path, which is found by adding the lengths of the critical paths needed
to compute q and to perform the addition rRS+aB−qM S. The addition, which
includes the digit multiplications, is like adding about 2ρ+6 binary numbers,
and so has a depth in the order of log2 (2ρ+6). However, binary carries mustthen propagate over the length of the digit, so that the critical path length has
an effective length of order ρ.

8International J. of Computer Mathematics, vol. 3, 1991, pp. 21-28

The number of bits for determining q has order ρ and so a look-up table
would require a critical path of the order log2 ρ, leaving the adder as the main
source of length in the critical path.
Thus the clock cycle time is asymptotically proportional to ρ. Also, the number of iterations in the loop is asymptotically proportional to the number ofdigits, and so to 1/ρ. This means that the throughput time is asymptoticallyconstant as r is increased, and no real benefit accrues from having a large base.
However, for small r, the various switching circuits etc. independent of r will
dominate critical path lengths, so that a modest increase in r beyond 2 should
initially reap a reward, but that return will decrease every time r is further
increased.
The area of the chip is dominated by that of the adder, and the registers
containing A, B, M and R. Again, the common hardware for any choice ofr dominates for very small r. The adder can be made from around 3(2ρ+6)gates for each bit position, giving an area 3N ρ(2ρ+6) for the adder. This is
proportional to ρ as r is increased. Hence chip efficiency measured as the inverseof the product of time by area actually decreases eventually as r is increased.
However, one immediate advantage of increasing the base is that the registers
contain fewer bits: for a given natural number, the number of bits required inour redundant systems is proportional to (ρ+1)/ρ, which makes the step frombase 2 to r = 4 or 8 quite attractive.

7

Conclusion

A very general modular multiplication algorithm has been presented and verified
and the efficiency of hardware implementations discussed. Although the order
estimates supply an accurate view of large-scale behaviour, lower order terms
dominate for small values of the base r of the number representations. It seems
best to choose the base r > 2, but not much larger.

References
1. Atkins, D.E., Higher-radix division using estimates of the divisor and partial remainders, IEEE Trans. Comp. C17 (1968) 925-934.
2. Avizienis, A., Signed Digit Number Representations for Fast Parallel Arithmetic,
IEEE Trans. Elec. Comp. EC10 (1961) 389-400.
3. Brickell, E.F., A Fast Modular Multiplication Algorithm with Application to Two
Key Cryptography, Advances in Cryptology (Proc. of CRYPTO 82), Chaum et al.
eds., Plenum, 1983, 51-60.
4. Gibson, J.K., A generalisation of Brickell’s algorithm for fast modular multiplication, Bit 28 (1988) 755-763.
5. Rivest, R.L., Shamir, A., & Adleman, L., A Method of obtaining Digital Signatures
and Public Key Cryptosystems, Comm. ACM 21 (1978) 120-126.
6. Selby, A., & Mitchell, C., Algorithms for software implementations of RSA, IEE
Proc. 136, part E, (1989) 166-170.

C. D. Walter

97. Simmons, G.J., & Norris, M.J., High speed arithmetic utilizing redundant number
systems, Nat. Telecom. Conf. 1980, Houston, Texas, IEEE, New York, 1980, 49.3/12.
8. Walter, C.D., & Eldridge, S.E., A Verification of Brickell’s Fast Modular Multiplication Algorithm, Intern. J. Computer Math. 33 (1990) 153-169.