Burt is chief scientist of RSA Laboratories, a division of RSA Data Security. He received a PhD in computer science from MIT in 1988 and is interested in cryptography and fast arithmetic techniques. He can be contacted at [email protected] and at 10 Twin Dolphin Drive, Redwood City, CA 94065.
Computer arithmetic gets better and faster all the time. Once you could only add 8-bit numbers, then 16 bits became the standard, now 32 bits, soon 64....
But what if you want to add, subtract, multiply, or divide 512-bit numbers? Few computers have machine-language instructions for big-number arithmetic and, for obvious reasons, even fewer programming languages have built-in operations supporting it.
The ability to manipulate 512-bit (or larger) numbers is necessary, for instance, when you're implementing mathematically oriented encryption schemes; both the RSA and DSS schemes involve such numbers. (For more-detailed information on RSA and DSS, see the sidebar entitled, "Public-Key Cryptography Meets the Real World" on page 22 of the May 1992 DDJ.) So the question arises, how can you operate on 512-bit numbers in a language such as C that only goes as far as 32 bits?
For the answer, we can look at one simple implementation of what is called "multiple-precision" (MP) arithmetic. MP arithmetic has been implemented in RSA Laboratories' cryptographic toolkit, RSAREF, and the code has been ported to many machines without modification. (RSAREF is available to U.S. and Canadian citizens at no charge for non-commercial use; contact RSA Laboratories, 10 Twin Dolphin Drive, Redwood City, CA 94065 or at rsaref @rsa.com.)
Representing MPs
MPs are represented as arrays of type NN_DIGIT, where NN_DIGIT depends on the machine. For instance, NN_DIGIT might be defined as an unsigned long, which has 32 bits on most machines. For technical reasons, the number of bits b must be even. Each element of the array is a digit in the base B representation of the MP, where B= 2{b}. The minimum digit is 0 and the maximum digit is B-1.
Lower-indexed elements of the array are less significant than higher-indexed elements. We call a[0] the 1s digit of array a, a[1] the Bs digit (similar to 10s digit), a[2] the B{2}s digit (similar to 100s digit), and so on.
For example, the ninth Fermat number, 2{512}+1, would be represented as an array of 17 32-bit NN_DIGITS:
a[0]=1 a[1]=a[2]=...=a[14]=a[15]=0 a[16]=1
Listing One (page 116) defines NN_DIGIT and some other items helpful in handling MPs, including function prototypes.
MP Tools
We start with four tools: C's built-in addition, subtraction, multiplication, and division operators. The tools let us:
- Add two NN_DIGITs, and get the 1s digit of the sum (but not the carry-out).
- Subtract an NN_DIGIT from an NN_DIGIT, and get the 1s digit of the remainder (but not the borrow-out).
- Multiply two NN_DIGITs, and get the 1s digit of the product (but not the Bs digit).
- Divide an NN_DIGIT by an NN_DIGIT, and get the quotient, also an NN_DIGIT.
First Step: MP Addition
In grade school you were taught that to add two numbers, you wrote them down, one above the other, then added columns of digits from right to left. You wrote down the 1s digit of the sum at the bottom of the column. If the sum had a 10s digit, you wrote it down at the top of the column to the left. This was the "carry," and the last carry became the leftmost sum digit.
Adding two MPs is much the same: We have a carry-in that's either 0 or 1 and two addend digits; we want a sum digit and a carry-out. But it's harder to get the carry-out than in the grade-school method because our tools only let us see the 1s digit of a sum. So let's take it a step at a time.
We first add the carry-in to the first addend digit, and look at the 1s digit of the sum. Let's call this the "subsum". Here's where we do a "twist" on the grade-school addition. We don't know directly whether there is a carry. But we do know that if there is, then the subsum must be less than the carry-in, because the real sum has wrapped past the maximum digit--but it can't get as far as the carry-in. (We detect carries this way throughout our MP implementation.)
So if the subsum is less than the carry-in, we have to carry out. We also know that the carry-in is 1, not 0 (otherwise we couldn't have carried out); and the subsum is 0 (since it's less than the carry-in). Since the subsum is 0, we write down the second addend digit as the sum digit, and go on to the next digit.
If we don't carry out right away, we next add the subsum to the second addend digit and look at the 1s digit of the sum, which we write down as the sum digit. If the sum digit is less than the second addend digit, we have to carry out. Otherwise, we don't. Figure 1 gives an example of this twist on grade-school addition. The first row (1s and 0s) are the carries; the third is the subsum. In the 1s column, the subsum is more than the carry, and the sum digit is more than the second addend digit, so there's no carry out. In the 10s and 100s column, the subsum is more than the carry, and the sum digit is less than the second addend digit, so there is a carry out. In the 1000s column, the subsum is less than the carry, so there is a carry out. The twist is that we can do everything with an addition operation that gives only the 1s digit of its result.
Figure 1: Grade-school addition with a twist.
11100 <--carry 9876 <--first addend 0976 <--subsum 5432 <--second addend 15308 <--sum
The procedure NN_Add (Listing Two, page 116) adds MPs with our limited tools. The addends are b and c, and the sum is a. They all have length digits. Index variable i moves through the digits from right to left, with carry as the carry-in and -out. The procedure returns the final carry-out.
NN_Add adds MPs with a loop around a simple three-part conditional. The first part computes the subsum. If there's a carry out, the first part sets the sum digit to the second addend digit, and stops. The second part computes the sum digit. If there's a carry-out, the second part sets carry to 1. The third part sets carry to 0 when there's no carry-out.
The sum can share memory with the addends, since NN_Add stores intermediate results in the variable ai. (The sum must share exactly the same memory; if it overlaps only partially, the result is undefined.) All procedures described here have the shared-memory feature.
MP Subtraction
Subtracting two MPs is just like adding two MPs, except that we borrow instead of carry. Here, we have a borrow-in that's either 0 or 1, a subtrahend digit, and a minuend digit (how's that for terminology!); we want a remainder digit and a borrow-out.
We first subtract the borrow-in from the subtrahend digit and look at the 1s digit of the remainder. We call this the "subremainder." If the subremainder is more than the maximum digit minus the borrow-in, we have to borrow out. We also know that the borrow-in is 1, and the subsum is the maximum digit. (This is just like addition, except everything's reversed: "Less" is "more," "plus" is "minus," "0" is "maximum digit.")
Since the subremainder is the maximum digit, we write down the maximum digit minus the minuend digit as the remainder digit and go on to the next digit.
If we don't borrow out right away, we next subtract the minuend digit from the subremainder and look at the 1s digit of the remainder, which we write down as the remainder digit. If the remainder digit is more than the maximum digit minus the minuend digit, we have to borrow out. Otherwise, we don't.
The procedure NN_Sub, shown in Listing Three (page 117), subtracts MPs. The subtrahend is b, the minuend is c, and the remainder is a. As in NN_Add, they all have length digits, and the index variable i moves through the digits. The variable borrow is the borrow-in and -out. The procedure returns the final borrow-out, which is 1 when the subtrahend is less than the minuend.
More Tools
Before going on to MP multiplication and addition, we need two additional tools. These add to those built into the C language, and let us:
- Multiply two NN_DIGITs and get a two-NN-DIGIT product with a 1s digit and a Bs digit.
- Divide a two -NN_DIGIT dividend with a 1s digit and a Bs digit by an NN_DIGIT and get the quotient, also an NN_DIGIT.
(These extra tools are in the instruction sets of most machines, even when NN_DIGIT is 32 bits, but they're not in Standard C. So if you're working in assembly language, your work is done. Read on to see how to do it in C.)
Digit Multiplication
We multiply two NN_DIGITs by combining four NN_HALF_DIGIT x NN_HALF_DIGIT multiplications. We can multiply two NN_HALF_DIGITs with the tools already built into the C language; the product is an NN_DIGIT. (The C language's type-conversion rules require that we cast the NN_HALF_DIGITs to an NN_DIGIT to get an NN_DIGIT product. A clever compiler will notice the casting and generate an NN_HALF_DIGITxNN_HALF_DIGIT machine-language multiplication.)
The procedure NN_DigitMult, shown in Listing Four (page xx), multiplies NN_DIGITs. The multiplicand and multiplier are b and c, and the product is a (a two-NN_DIGIT array).
A familiar algebraic formula describes the mapping to four multiplications. Let bHigh, bLow, cHigh, and cLow denote the high-order and low-order halves of b and c, so that
b=bHighx\/B+bLow ;
c=cHighx\/B+cLow
Then we have:
bxc=bHighxcHighxB+((bHighxcLow)
+(bLowxcHigh)squareroot of B+bLowxcLow
NN_DigitMult computes the four multiplications and stores the product of low-order halves in a[0], the product of high-order halves in a[1], and the two cross products in t and u.
The low-order halves of the cross products need to be added into the high-order half of a[0]. Similarly, the high-order halves need to be added into the low-order half of a[1].
NN_DigitMult adds the intermediate values t and u. If the sum is less than u--there is a carry out -- the procedure adds 1 to the high-order half of a[1]. NN_DigitMult then adds the low-order half of t+u to a[0], carrying into a[1]. It finally adds the high-order half of t+u to a[1].
Digit Division
We divide a two-NN_DIGIT dividend by an NN_DIGIT divisor by combining two NN_DIGITxNN_HALF_DIGIT divisions, four NN_HALF_DIGITxNN_HALF_DIGIT multiplications, and some other arithmetic. We can divide an NN_DIGIT dividend by an NNHALF_DIGIT divisor with the C-language tools; the quotient is an NN_HALF_DIGIT, assuming the high half of the dividend is less than the divisor. (It's worth mentioning that we're interested in the "floor" of the quotient--the greatest integer not more than the quotient.)
The procedure NN_DigitDiv, shown in Listing Five (page 117), divides a two-NN_DIGIT dividend by an NN_DIGIT divisor. The dividend is b (a two-NN_DIGIT array), the divisor is c, and the quotient is a. The two-NN_DIGIT variable t holds a copy of the dividend. NN_DigitDiv assumes that the quotient is not more than the maximum digit. We guarantee this later when we call NN_DigitDiv.
NN_DigitDiv follows the estimate-and-correct method. First, it estimates the high-order half of the quotient with an NN_DIGITxNN_HALF_DIGIT division. It then subtracts the product of the estimate and the divisor from the dividend. It finally corrects the estimate by subtracting the divisor (shifted left a half digit) from the dividend, until the dividend is less than the divisor (shifted left a half digit). Next, NN-DigitDiv estimates the low-order half of the quotient, given what remains of the dividend; it subtracts the product of the estimate and the divisor from the dividend, and it finally corrects the estimate. Usually there's only one correction for each half, so NN_DigitDiv gets the correct quotient pretty fast.
NN_DigitDiv estimates the high-order half of the quotient by dividing t[1] by the high-order half of c, plus 1. (If the high-order half of c were the maximum half digit, division would be by square root of B, which is too large, so the estimate is just the high-order half of t[1].)
NN_DigitDiv subtracts the product of the estimate aHigh and the divisor, shifted left a half digit, from t. Then it subtracts the divisor, shifted left a half digit, and increments the estimate until what remains of t is less than the divisor shifted left a half digit.
The assumption that the quotient is not more than the maximum half digit ensures that the estimate doesn't exceed the maximum half digit.
NN_DigitDiv then estimates the low-order half of the quotient by similar means. It divides the low-order half of t[1], shifted left a half digit, plus the high-order half of t[0], by the high-order half of c, plus 1. (If the high-order half of c is the maximum half digit, the estimate is just the low-order half of t[1].) Again, the estimate aLow is an underestimate, and the procedure corrects it.
It's a good idea to normalize the divisor by shifting as far left as possible before calling NN_DigitDiv, because otherwise it can take a long time to correct the estimates. If the divisor is normalized, then the initial estimate can't be off by more than 2. NN_DigitDiv thus makes at most two and usually one correction. We guarantee that the divisor is normalized, too, when we call NN_DigitDiv later.
MP Multiplication
Returning to grade school, remember that to multiply two numbers (a multiplier and a multiplicand), you wrote them down, one above the other. Then you took the 1s digit of the multiplier, multiplied the multiplicand by it, and wrote down the product. You took the 10s digit of the multiplier, multiplied the multiplicand by it, and wrote down the product, below and one digit to the left of the first product. And so on; after you wrote down all the intermediate products, you added them up, and the sum became the final product.
Now let's consider grade-school multiplication with a twist. Instead of waiting until the end to add, suppose we accumulate the sum as we go along. Then we don't have to store the intermediate products, just the accumulator, and the accumulator becomes the final product.
Multiplying two MPs thus involves two procedures: one that multiplies an MP by a digit and adds the product to an accumulator, and one that moves through the digits of the multiplier.
The procedure NN_AddDigitMult, shown in Listing Six (page 118), does the first part. The digit is c, the multiplicand is d, the input accumulator is b, and the output is a. The multiplicand and accumulators all have length digits. Index variable i moves through the digits of the multiplicand from right to left, with carry as the carry-in and -out. NN_AddDigitMult returns the final carry-out.
NN_AddDigitMult is similar to NN_Add, except for two things. It calls the tool NN_DigitMult to multiply digits, rather than adding them with a built-in NN_DIGIT add (since we're multiplying, not adding); and carry is a digit, not just 0 or 1. Here there are two conditionals. The first conditional adds the carry. If there's a carry-out, the first conditional sets the carry to 1; otherwise, it sets the carry to 0. The second conditional adds the low-order digit of the NN_DIGITxNN_DIGIT product. If there's a carry out, the second conditional increments the carry. NN_AddDigitMult then unconditionally adds the high-order digit of the NN_DIGITxNN_DIGIT product to the carry and continues.
The procedure NN_Mult, shown in Listing Seven (page 118), multiplies MPs. The multiplier is b, the multiplicand is c, and the product is a. The inputs have length digits and the output has length 2xdigits. Index variable i moves through the digits of the multiplier from right to left. In fact it doesn't go all the way to the left, only as far as the leftmost non-zero digit. (The procedure NN_Digits, in Listing Ten, page 119, determines how far to go.) NN_Mult similarly reduces the length of the multiplicand to avoid unnecessary multiplications. It's easy to see how NN_Mult works with NN_AddDigitMult to accumulate the product. Each call has a new multiplier digit, b[i], and operates on a new part of the accumulator, &t[i]. NN_Mult stores NN_AddDigitMult's carry as a new high-order digit of the accumulator. This is all quite similar to the illustration of grade-school multiplication that Figure 2 describes, in which the first row is the multiplicand, the second is the multiplier, and the last is the product. Intermediate rows give accumulator values and intermediate products. Intermediate products are between a multiplier digit and the multiplicand. They are shifted left as many digits as the multiplier digit. The accumulator value is the sum of the previous accumulator value and an intermediate product. Intermediate products can be added directly into the accumulator. (NN_Mult also calls two other procedures we haven't discussed: NN_Assign, which copies an MP, and NN_AssignZero, which sets an MP to 0. See Listing Ten.)
Figure 2: Grade-school multiplication with a twist.
5432 <-- multiplicand
9876 <-- multiplier
________
0000 <-- accumulator value
32592 <-- intermediate product
________
32592 <-- accumulator value
38024 <-- intermediate product
________
412832
43456
________
4758432
48880
________
53646432 <-- final product
Last Step: MP Division
We're now ready to put together all our tools to do what's generally considered the most difficult arithmetic operation: division.
We've already seen how to divide a two-NN_DIGIT dividend by an NN_DIGIT divisor. The estimate-and-correct method has served us well. We extend that method to MP division.
For each digit of the quotient, we first estimate the digit by dividing the high-order two digits of the dividend (or what remains of it) by the high-order digit of the divisor, plus 1. (If the high-order digit of the divisor is the maximum digit, the estimate is just the high-order digit of the dividend.) We then subtract the product of the estimate and the divisor, shifted left some number of digits, from what remains of the divisor.
The estimate is an underestimate, so we may have to increase it, just as in NN_DigitDiv. We do so by subtracting the divisor, shifted left some number of digits, from the dividend, until what remains of the dividend is less than the shifted divisor. When we reach the 1s digit of the quotient, we have a remainder that's less than the divisor, and we're done. Figure 3 illustrates grade-school long division with the latest twist. The first row is the quotient. The first part of the second row is the divisor, and the second part is the dividend. The last row is the remainder. Quotient digits are estimated by dividing the two high-order digits of what remains of the dividend by the high-order digit of the divisor, plus 1. If what remains after subtracting a shifted multiple of the divisor is too large, the estimate is corrected. For instance, the digits 53 are divided by (5 + 1) to get an estimate 8 for the 1000s digit of the quotient. What remains after subtracting 8000x5432 is too large, so the estimate is incremented, and another 1000x5432 is subtracted.
Figure 3: Grade-school division with a twist.
09876
_________
5432)053646432
0000 <-- estimate: 5/6 = 0
_________
53646432
43456 <-- estimate: 53/6 = 8
_________
10190432 (too large)
5432 <-- corrected estimate: 9
_________
4758432
38024 <-- estimate: 47/6 = 7
_________
956032 (too large)
5432 <-- corrected estimate: 8
_________
412832
32592 <-- estimate: 41/6 = 6
_________
86912 (too large)
5432 <-- corrected estimate: 7
_________
32592
27160 <-- estimate: 32/6 = 5
_________
5432 (too large)
5432 <-- corrected estimate: 6
_________
0 <-- remainder
Dividing MPs, like multiplying MPs, involves two new procedures: one that multiplies an MP by a digit and subtracts the product from an accumulator, and one that moves through the digits of the quotient.
The procedure NN_SubDigitMult, shown in Listing Eight (page 118), does the first part. It is just like NN_AddDigitMult. NN_SubDigitMult returns the final borrow-out.
The procedure NN_Div, shown in Listing Nine (page 118), divides MPs. The dividend is c, the divisor is d, the quotient is a, and the remainder is b. The dividend and quotient have length cDigits, while the divisor and remainder have length dDigits. (Notice that NN-Div returns both quotient and remainder, whereas NN_DigitDiv returns only the quotient.)
The index variable i moves through the digits of the quotient from left to right. (This is the first time we've moved in that direction, and we do so because division reveals the most significant digit first.) NN_Div estimates the quotient digit with NN_DigitDiv, subtracts the product of the estimate and the divisor with NN_SubDigitMult, and corrects the estimate with NN_Cmp (MP comparison, see Listing Ten) and NN_Sub.
NN_Div is a little more complicated than the other procedures; this is due to NN_DigitDiv's two requirements. For efficiency, NN_Div normalizes the divisor. It does this with NN_LShift (MP left shift, Listing Ten). The number of bits to shift left is determined with NN_DigitBits. NN_Div shifts both the divisor and the dividend, and at the end, shifts back the remainder with NN_RShift (MP right shift, Listing Ten). The quotient is unaffected by the normalization.
For correctness, NN_Div extends the dividend with a leading 0 digit. This guarantees that the leading digit of the dividend is less than the leading digit of the divisor, as NN_DigitDiv requires. The condition continues throughout the main loop.
Conclusions
In this article, we've gone from built-in C tools to multiple-precision addition, subtraction, multiplication, and division. One possible next step is modular arithmetic, where all results are divided by a predetermined quantity called the "modulus," and only the remainder is kept. Modular arithmetic, and modular exponentiation in particular, are essential to the RSA and DSS schemes. You can also implement some number-theoretic operations such as greatest common divisor with the tools presented here. Other next steps include converting an MP to base 10--useful if you want to print your results--and even computing the digits of pi.
It's not difficult to implement the arithmetic operations more efficiently than we've done here, since we've tried to keep things simple. Assembly-language speedups would help a lot. Some technical refinements would help, too, and indeed are available in various commercial and public-domain implementations. Better and faster implementations are encouraged.
Editor's Note
Public-key cryptosystems such as RSA and DSS are subject to U.S. patents 4,218,582 and 4,405,829, as well as other patents issued to Stanford University and MIT. Public Key Partners of Sunnyvale, California holds exclusive licensing rights. Licensed software embodying the cryptosystems is available from several vendors, including RSA Data Security.
Recommended Reading
Baker, Henry G. "Computing A*B (mod N) Efficiently in ANSI C." ACM SIGPLAN Notices 27 (January, 1992).
Buell, Duncan A. and Robert L. Ward. "A Multiprecise Integer Arithmetic Package." The Journal of Supercomputing (vol. 3, 1989).
Dusse, Stephen R. and Burton S. Kaliski, Jr. "A Cryptographic Library for the Motorola DSP56000." Advances in Cryptology--EUROCRYPT '90 Proceedings, volume 473 of Lecture Notes in Computer Science, New York: Springer-Verlag, 1991.
Knuth, Donald E. "Seminumerical Algorithms." The Art of Computer Programming, Volume 2, Second edition. Reading, MA: Addison-Wesley, 1983.
Montgomery, Peter L. "Modular Multiplication Without Trial Division" Mathematics of Computation (vol. 44, 1985).
_MULTIPLE-PRECISION ARITHMETIC IN C_ by Burton S. Kaliski, Jr.[LISTING ONE]
<a name="01b6_0014">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* PROTOTYPES should be set to one if and only if the compiler supports
function argument prototyping. The following makes PROTOTYPES default to 0
if it has not already been defined with C compiler flags. */
#ifndef PROTOTYPES
#define PROTOTYPES 0
#endif
/* UINT2 defines a two byte word */
typedef unsigned short int UINT2;
/* UINT4 defines a four byte word */
typedef unsigned long int UINT4;
/* PROTO_LIST is defined depending on how PROTOTYPES is defined above. If using
PROTOTYPES, PROTO_LIST returns the list, otherwise it returns empty list. */
#if PROTOTYPES
#define PROTO_LIST(list) list
#else
#define PROTO_LIST(list) ()
#endif
/* RSA key lengths. */
#define MAX_RSA_MODULUS_BITS 1024
#define MAX_RSA_MODULUS_LEN ((MAX_RSA_MODULUS_BITS + 7) / 8)
typedef UINT4 NN_DIGIT;
typedef UINT2 NN_HALF_DIGIT;
/* Length of digit in bits */
#define NN_DIGIT_BITS 32
#define NN_HALF_DIGIT_BITS 16
/* Length of digit in bytes */
#define NN_DIGIT_LEN (NN_DIGIT_BITS / 8)
/* Maximum length in digits */
#define MAX_NN_DIGITS \
((MAX_RSA_MODULUS_LEN + NN_DIGIT_LEN - 1) / NN_DIGIT_LEN + 1)
/* Maximum digits */
#define MAX_NN_DIGIT 0xffffffff
#define MAX_NN_HALF_DIGIT 0xffff
/* Macros. */
#define LOW_HALF(x) (NN_HALF_DIGIT)((x) & MAX_NN_HALF_DIGIT)
#define HIGH_HALF(x) \
(NN_HALF_DIGIT)(((x) >> NN_HALF_DIGIT_BITS) & MAX_NN_HALF_DIGIT)
#define TO_HIGH_HALF(x) (((NN_DIGIT)(x)) << NN_HALF_DIGIT_BITS)
#define DIGIT_MSB(x) (unsigned int)(((x) >> (NN_DIGIT_BITS - 1)) & 1)
#define DIGIT_2MSB(x) \
(unsigned int)(((x) >> (NN_DIGIT_BITS - 2)) & 3)
/* NOTE: A bug in the MPW 3.2 C compiler causes an incorrect sign extension in
the routine NN_DigitDiv. To overcome this bug, change the definition of the
macro HIGH_HALF to:
#define HIGH_HALF(x) (((x) >> NN_HALF_DIGIT_BITS) & MAX_NN_HALF_DIGIT) */
void NN_Assign PROTO_LIST ((NN_DIGIT *, NN_DIGIT *, unsigned int));
void NN_AssignZero PROTO_LIST ((NN_DIGIT *, unsigned int));
NN_DIGIT NN_Add PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT *, unsigned int));
NN_DIGIT NN_Sub PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT *, unsigned int));
void NN_Mult PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT *, unsigned int));
int NN_Cmp PROTO_LIST ((NN_DIGIT *, NN_DIGIT *, unsigned int));
unsigned int NN_Bits PROTO_LIST ((NN_DIGIT *, unsigned int));
unsigned int NN_Digits PROTO_LIST ((NN_DIGIT *, unsigned int));
NN_DIGIT NN_LShift PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, unsigned int, unsigned int));
NN_DIGIT NN_RShift PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, unsigned int, unsigned int));
void NN_Div PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT *, unsigned int, NN_DIGIT *,
unsigned int));
NN_DIGIT NN_AddDigitMult PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT, NN_DIGIT *, unsigned int));
NN_DIGIT NN_SubDigitMult PROTO_LIST
((NN_DIGIT *, NN_DIGIT *, NN_DIGIT, NN_DIGIT *, unsigned int));
unsigned int NN_DigitBits PROTO_LIST ((NN_DIGIT));
void NN_DigitMult PROTO_LIST ((NN_DIGIT [2], NN_DIGIT, NN_DIGIT));
void NN_DigitDiv PROTO_LIST ((NN_DIGIT *, NN_DIGIT [2], NN_DIGIT));
<a name="01b6_0015">
<a name="01b6_0016">[LISTING TWO]<a name="01b6_0016">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b+c. Returns carry. Lengths: a[digits], b[digits], c[digits]. */
NN_DIGIT NN_Add (a, b, c, digits)
NN_DIGIT *a, *b, *c;
unsigned int digits;
{
NN_DIGIT ai, carry;
unsigned int i;
carry = 0;
for (i = 0; i < digits; i++) {
if ((ai = b[i] + carry) < carry)
ai = c[i];
else if ((ai += c[i]) < c[i])
carry = 1;
else
carry = 0;
a[i] = ai;
}
return (carry);
}
<a name="01b6_0017">
<a name="01b6_0018">[LISTING THREE]<a name="01b6_0018">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b-c. Returns borrow. Lengths: a[digits], b[digits], c[digits]. */
NN_DIGIT NN_Sub (a, b, c, digits)
NN_DIGIT *a, *b, *c;
unsigned int digits;
{
NN_DIGIT ai, borrow;
unsigned int i;
borrow = 0;
for (i = 0; i < digits; i++) {
if ((ai = b[i] - borrow) > (MAX_NN_DIGIT - borrow))
ai = MAX_NN_DIGIT - c[i];
else if ((ai -= c[i]) > (MAX_NN_DIGIT - c[i]))
borrow = 1;
else
borrow = 0;
a[i] = ai;
}
return (borrow);
}
<a name="01b6_0019">
<a name="01b6_001a">[LISTING FOUR]<a name="01b6_001a">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b*c, where b and c are digits. Lengths: a[2]. */
void NN_DigitMult (a, b, c)
NN_DIGIT a[2], b, c;
{
NN_DIGIT t, u;
NN_HALF_DIGIT bHigh, bLow, cHigh, cLow;
bHigh = HIGH_HALF (b);
bLow = LOW_HALF (b);
cHigh = HIGH_HALF (c);
cLow = LOW_HALF (c);
a[0] = (NN_DIGIT)bLow * (NN_DIGIT)cLow;
t = (NN_DIGIT)bLow * (NN_DIGIT)cHigh;
u = (NN_DIGIT)bHigh * (NN_DIGIT)cLow;
a[1] = (NN_DIGIT)bHigh * (NN_DIGIT)cHigh;
if ((t += u) < u)
a[1] += TO_HIGH_HALF (1);
u = TO_HIGH_HALF (t);
if ((a[0] += u) < u)
a[1]++;
a[1] += HIGH_HALF (t);
}
<a name="01b6_001b">
<a name="01b6_001c">[LISTING FIVE]<a name="01b6_001c">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Sets a=b/c, where a and c are digits. Lengths: b[2]. Assumes b[1] < c
and HIGH_HALF (c) > 0. For efficiency, c should be normalized. */
void NN_DigitDiv (a, b, c)
NN_DIGIT *a, b[2], c;
{
NN_DIGIT t[2], u, v;
NN_HALF_DIGIT aHigh, aLow, cHigh, cLow;
cHigh = HIGH_HALF (c);
cLow = LOW_HALF (c);
t[0] = b[0];
t[1] = b[1];
/* Underestimate high half of quotient and subtract product
of estimate and divisor from dividend. */
if (cHigh == MAX_NN_HALF_DIGIT)
aHigh = HIGH_HALF (t[1]);
else
aHigh = (NN_HALF_DIGIT)(t[1] / (cHigh + 1));
u = (NN_DIGIT)aHigh * (NN_DIGIT)cLow;
v = (NN_DIGIT)aHigh * (NN_DIGIT)cHigh;
if ((t[0] -= TO_HIGH_HALF (u)) > (MAX_NN_DIGIT - TO_HIGH_HALF (u)))
t[1]--;
t[1] -= HIGH_HALF (u);
t[1] -= v;
/* Correct estimate. */
while ((t[1] > cHigh) ||
((t[1] == cHigh) && (t[0] >= TO_HIGH_HALF (cLow)))) {
if ((t[0] -= TO_HIGH_HALF (cLow))
> MAX_NN_DIGIT - TO_HIGH_HALF (cLow))
t[1]--;
t[1] -= cHigh;
aHigh++;
}
/* Underestimate low half of quotient and subtract product of
estimate and divisor from what remains of dividend. */
if (cHigh == MAX_NN_HALF_DIGIT)
aLow = LOW_HALF (t[1]);
else
aLow =
(NN_HALF_DIGIT)
((NN_DIGIT)(TO_HIGH_HALF (t[1]) + HIGH_HALF (t[0]))
/ (cHigh + 1));
u = (NN_DIGIT)aLow * (NN_DIGIT)cLow;
v = (NN_DIGIT)aLow * (NN_DIGIT)cHigh;
if ((t[0] -= u) > (MAX_NN_DIGIT - u))
t[1]--;
if ((t[0] -= TO_HIGH_HALF (v)) > (MAX_NN_DIGIT - TO_HIGH_HALF (v)))
t[1]--;
t[1] -= HIGH_HALF (v);
/* Correct estimate. */
while ((t[1] > 0) || ((t[1] == 0) && t[0] >= c)) {
if ((t[0] -= c) > (MAX_NN_DIGIT - c))
t[1]--;
aLow++;
}
*a = TO_HIGH_HALF (aHigh) + aLow;
}
<a name="01b6_001d">
<a name="01b6_001e">[LISTING SIX]<a name="01b6_001e">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b+c*d, where c is a digit. Returns carry.
Lengths: a[digits], b[digits], d[digits]. */
NN_DIGIT NN_AddDigitMult (a, b, c, d, digits)
NN_DIGIT *a, *b, c, *d;
unsigned int digits;
{
NN_DIGIT carry, t[2];
unsigned int i;
carry = 0;
for (i = 0; i < digits; i++) {
NN_DigitMult (t, c, d[i]);
if ((a[i] = b[i] + carry) < carry)
carry = 1;
else
carry = 0;
if ((a[i] += t[0]) < t[0])
carry++;
carry += t[1];
}
return (carry);
}
<a name="01b6_001f">
<a name="01b6_0020">[LISTING SEVEN]<a name="01b6_0020">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b*c. Lengths: a[2*digits], b[digits], c[digits].
Assumes digits < MAX_NN_DIGITS. */
void NN_Mult (a, b, c, digits)
NN_DIGIT *a, *b, *c;
unsigned int digits;
{
NN_DIGIT t[2*MAX_NN_DIGITS];
unsigned int bDigits, cDigits, i;
NN_AssignZero (t, 2 * digits);
bDigits = NN_Digits (b, digits);
cDigits = NN_Digits (c, digits);
for (i = 0; i < bDigits; i++)
t[i+cDigits] += NN_AddDigitMult (&t[i], &t[i], b[i], c, cDigits);
NN_Assign (a, t, 2 * digits);
}
<a name="01b6_0021">
<a name="01b6_0022">[LISTING EIGHT]<a name="01b6_0022">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=b-c*d, where c is a digit. Returns borrow.
Lengths: a[digits], b[digits], d[digits]. */
NN_DIGIT NN_SubDigitMult (a, b, c, d, digits)
NN_DIGIT *a, *b, c, *d;
unsigned int digits;
{
NN_DIGIT borrow, t[2];
unsigned int i;
borrow = 0;
for (i = 0; i < digits; i++) {
NN_DigitMult (t, c, d[i]);
if ((a[i] = b[i] - borrow) > (MAX_NN_DIGIT - borrow))
borrow = 1;
else
borrow = 0;
if ((a[i] -= t[0]) > (MAX_NN_DIGIT - t[0]))
borrow++;
borrow += t[1];
}
return (borrow);
}
<a name="01b6_0023">
<a name="01b6_0024">[LISTING NINE]<a name="01b6_0024">
/* Copyright (C) 1991-2 RSA Laboratories, a division of RSA Data
Security, Inc. All rights reserved. */
/* Computes a=c div d and b=c mod d. Lengths: a[cDigits], b[dDigits],
c[cDigits], d[dDigits]. Assumes d > 0, cDigits < 2 * MAX_NN_DIGITS,
dDigits < MAX_NN_DIGITS. */
void NN_Div (a, b, c, cDigits, d, dDigits)
NN_DIGIT *a, *b, *c, *d;
unsigned int cDigits, dDigits;
{
NN_DIGIT ai, cc[2*MAX_NN_DIGITS+1], dd[MAX_NN_DIGITS], t;
int i;
unsigned int ddDigits, shift;
ddDigits = NN_Digits (d, dDigits);
if (ddDigits == 0)
return;
<a name="01b6_0025"><a name="01b6_0026">[LISTING TEN]<a name="01b6_0026">
/* Normalize operands. */
shift = NN_DIGIT_BITS - NN_DigitBits (d[ddDigits-1]);
NN_AssignZero (cc, ddDigits);
cc[cDigits] = NN_LShift (cc, c, shift, cDigits);
NN_LShift (dd, d, shift, ddDigits);
t = dd[ddDigits-1];
NN_AssignZero (a, cDigits);
for (i = cDigits-ddDigits; i >= 0; i--) {
/* Underestimate quotient digit and subtract. */
if (t == MAX_NN_DIGIT)
ai = cc[i+ddDigits];
else
NN_DigitDiv (&ai, &cc[i+ddDigits-1], t + 1);
cc[i+ddDigits] -=
NN_SubDigitMult (&cc[i], &cc[i], ai, dd, ddDigits);
/* Correct estimate. */
while (cc[i+ddDigits] || (NN_Cmp (&cc[i], dd, ddDigits) >= 0)) {
ai++;
cc[i+ddDigits] -= NN_Sub (&cc[i], &cc[i], dd, ddDigits);
}
a[i] = ai;
}
/* Restore result. */
NN_AssignZero (b, dDigits);
NN_RShift (b, cc, shift, ddDigits);
}
Copyright © 1992, Dr. Dobb's Journal
