## Contents

## Introduction

The implementation of asymmetrical cryptographic schemes often requires
the use of numbers that are many times larger than the integer data types that are
supported natively by the compiler. In this article, we give an introduction to the
implementation of arithmetic operations involving large integers. We do not attempt
to give a full coverage of this topic since it is both complex and lengthy. For a
more detailed treatment, the reader is referred to the listed references.

The source code that accompanies this article implements the `BigInteger`

class supporting large integer arithmetic operations. Overloaded operators
includes +, -, *, /, %, >>, <<, ==, !=, >, <, >=, <=, &, |, ^, ++, -- and ~.
Other additional features such as modular exponential, modular inverse,
pseudoprime generation and probabilistic primality testing are also supported.

## Features

- Arithmetic operations involving large signed integers in 2's complement
representation.
- Prime number tests using Fermat's Little Theorem, Rabin Miller's method
and Solovay Strassen's method.
- Modular exponential with Barrett reduction.
- Modular inverse.
- Random Pseudoprime generation.
- Random Coprime generation.
- Greatest common divisor.

## Implementation Details

The most common way of representing numbers is by using the positional
notation system. Numbers are written using digits to represent multiples
of powers of the specified base. The base that we are most familiar with
and use everyday, is base 10. When we write the number 12345 in base 10,
it actually means

12345_{10}= 1x10^{4} + 2x10^{3} +
3x10^{2} + 4x10^{1} + 5x10^{0}

Similarly, if the same number is specified in base 16, then

12345_{16}= 1x16^{4} + 2x16^{3} +
3x16^{2} + 4x16^{1} + 5x16^{0}

Hence, it is important to know the base that the number is specified in,
since the same representation could represent different values when interpreted
in different bases. In general, a positive number can be written as a
polynomial of order k, where k is non-negative and a_{k} != 0.

n = a_{k}b^{k} + a_{k-1}b^{k-1} + ... +
a_{1}b^{1} + a_{0}

b represents the base with 0 <= a <= b-1.

For base 10 representation, b = 10 and the digits allowed at each
position is 0 <= a <= 9.

### Addition of Big Integers

The addition of two numbers can be represented as the addition of two
polynomials as shown below.

Let

n = a_{k}b^{k} + a_{k-1}b^{k-1} + ... +
a_{1}b^{1} + a_{0}

m = c_{k}b^{k} + c_{k-1}b^{k-1} + ... +
c_{1}b^{1} + c_{0}

Then

n + m = (a_{k} + c_{k})b^{k} + (a_{k-1} +
c_{k-1})b^{k-1} + ... + (a_{1} +
c_{1})b^{1} + (a_{0} + c_{0})

However, it is often easier to visualize the addition of two large
numbers as the digit-by-digit addition at each position. When we add
the digits at a particular position, the largest resulting value is 2b - 2.

#### Proof

Since the maximum value
of each digit is b - 1, we have

(b - 1) + (b - 1) = 2b - 2.

Since 2b - 2 < 2b, the maximum value that we have to carry over to the
next higher position is 1. For example, if we add 9 and 9 in base 10, we get
18, where 8 remains in the current position and the 1 gets *carried out*
to the next higher position. When adding the two digits in the next position, we must
remember to add the carry as well.

In our `BigInteger`

implementation, we store the number as an
array of consecutive bytes. This is equivalent to the representation of the
number in base 256 with each digit having values 0 to FF. A numerical
example that shows the addition of two numbers in this base is given
below.

1 1
FE A6 FF
+ FF FF
-----------
FF A6 FE
-----------

In this example, FF + FF = 1FE. This exceeds the maximum value of each
digit in base 256. Hence, we retain FE in the current position and carry 1
to the next higher position. This code segment that shows how this is done
is given below.

int carry = 0;
for(int i = 0; i < len; i++)
{
int sum = bi1.data[i] + bi2.data[i] + carry;
carry = sum >> 8;
result.data[i] = (byte)(sum & 0xFF);
}

### Subtraction of Big Integers

The subtraction of two numbers is very similar to addition and can be
easily implemented as byte-by-byte subtraction. However, when subtracting
two digits, we often encounter the situation where the first
digit is smaller than the second. In this case, we have to *borrow* 1
from the digit in the next position. This is actually done by subtracting
1 from the digit in the next higher position and adding
the value of the base to the current digit. For example,

2 AF
- FF
-----
1 B0
-----

In this case, since AF < FF, we subtract 1 from the next higher digit 2
and compute 1AF - FF = B0. The following code segment shows how byte-by-byte
subtraction is performed.

int carryIn = 0;
for(int i = 0; i < len; i++)
{
int diff, totalSub;
totalSub = bi2.data[i] + carryIn;
if(bi1.data[i] < totalSub)
{
diff = (bi1.data[i] + 0x100) - totalSub;
carryIn = 1;
}
else
{
diff = bi1.data[i] - totalSub;
carryIn = 0;
}
result.data[i] = (byte)diff;
}

### Multiplication of Big Integers

The simplest way to implement multiplication is by using repeated
addition. For simplicity of analysis, we assume that both the multiplier
and the multiplicand has the same number of digits.
To compute *a * b*, we set the result to 0, then we repeatedly add
*a* to the result for *b* number of times.
Although this approach is simple, its complexity is O(b) large integer
additions and is equivalent to O(bn)
byte-by-byte additions, where *n* is the number of digits in the large
integer. For large *b*, O(bn) > O(n^{2}). Hence, this approach does
not scale to large multipliers.

The next easiest approach is by using techniques that we have learnt in
elementary schools. More precisely, we multiply *a* by the each digit
of *b* and sum the partial results. In the following example, we illustrate
in base 256 the computation of 02AF x 01FF.

We first compute 02AF x FF

AE
02 AF
x FF
-------
2 AC 51
-------

Note that when we multiply AF by FF, we get AE51. Since this exceeds
base 256, we retain 51 in the current position and carry AF over to the next
higher position. We then proceed to the next position and
multiply 02 by FF, remembering to add the carry to get 02AC.

Next, we compute 02AF x 01

02 AF
x 01
-------
02 AF
-------

Finally, we sum the partial results.

1
02 AC 51 (result from 02AF x FF)
+ 02 AF 00 (result from 02AF x 01 with 00 appended)
----------
05 5B 51
----------

From this illustration, it is clear that the time complexity of
multiplying two large integers is O(n^{2}) byte multiplications where n is
the number of digits in the number. The process of adding the partial results has
a complexity of O(n) large integer additions, equivalent to O(n^{2})
byte additions. Thus, the total complexity is O(n^{2}) byte multiplications
and O(n^{2}) byte additions. We have implemented this approach in our
`BigInteger`

class as shown in the following code segment.

for(int i = 0; i < bi1.byteLength; i++)
{
int mcarry = 0;
for(int j = 0; j < bi2.byteLength; j++)
{
int val = (bi1.data[i] * bi2.data[j]) +
tempResult[i+j] + mcarry;
tempResult[i+j] = (byte)(val & 0xFF);
mcarry = (val >> 8);
}
tempResult[i+bi2.byteLength] = (byte)mcarry;
}

### Division of Big Integers

Please refer to the listed references.

### Primality Testing using Fermat's Little Theorem

The purpose of primality testing is to find out whether a given number is
prime. A prime number is a number that is only divisible by 1 and itself. It is
easy to determine the primality of a small number *n*
by testing whether it is divisible by a prime number that is less than
*sqrt(n)*. For a large integer, this method is not practical. Instead,
we can only determine whether the large integer is *probably prime*.
This is known as *probabilistic primality testing* and numbers that
passes the test are known as *pseudoprimes*. There are many
*probabilistic primality testing* methods, but in this article,
we will discuss only the method based on Fermat's Little Theorem.

From Fermat's Little Theorem, we know that

a^{n} mod n = a

for any integer *a* and any prime *n*. This means that if
*n* is prime then *a*^{n} mod n = a for
any integer *a*. The contrapositive of this statement says that if
*a*^{n} mod n != a for some integer *a* then *n* is
NOT prime.

Unfortunately, the converse of the theorem is not true. This means that
if *a*^{n} mod n = a, we cannot be sure whether n is prime.
However, we can make use of the contrapositive to determine if *n* is
composite. In other words, if there exists an integer *a* such that
*a*^{n} mod n != a, then *n* is not prime.
By testing out sufficient numbers of *a's*, we can exclude numbers that
are composite and retain numbers that
are probably prime. The following source code illustrates this.

for(int round = 0; round < confidence; round++)
{
a.genRandomBits(testBits, rand);
BigInteger expResult = a.modPow(thisVal, thisVal);
if(expResult != a)
return false;
}
return true;

There is a class of composite integer known as Carmichael number that
cannot be differentiated from real primes using this test. Such integers satisfy
*a*^{n} mod n = a for all positive integers *a*
with *gcd(a,n) = 1*.

## Sample Code

The generation of a 512-bits random pseudoprime number.

public static void Main(string[] args)
{
Console.Write(
"\nGenerating 512-bits random pseudoprime. . .");
Random rand = new Random();
BigInteger prime = BigInteger.genPseudoPrime(
512, 5, rand);
Console.WriteLine("\n" + prime);
}

Arithmetic operations involving `BigInteger`

.

public static void Main(string[] args)
{
BigInteger bn1 = new BigInteger("12345678901234567890",
10);
BigInteger bn2 = new BigInteger("4321", 10);
BigInteger bn3 = bn1 / bn2;
BigInteger bn4 = bn1 % bn2;
BigInteger bn5 = (bn3 * bn2) + bn4;
Console.WriteLine("bn1 = " + bn1);
Console.WriteLine("bn5 = " + bn5);
}

Asymmetrical encryption and decryption.

public static void Main(string[] args)
{
BigInteger bi_e = new BigInteger(
"a932b948feed4fb2b692609bd22164fc9edb" +
"59fae7880cc1eaff7b3c9626b7e5b241c27a" +
"974833b2622ebe09beb451917663d4723248" +
"8f23a117fc97720f1e7", 16);
BigInteger bi_d = new BigInteger(
"4adf2f7a89da93248509347d2ae506d683dd" +
"3a16357e859a980c4f77a4e2f7a01fae289f" +
"13a851df6e9db5adaa60bfd2b162bbbe31f7" +
"c8f828261a6839311929d2cef4f864dde65e" +
"556ce43c89bbbf9f1ac5511315847ce9cc8d" +
"c92470a747b8792d6a83b0092d2e5ebaf852" +
"c85cacf34278efa99160f2f8aa7ee7214de07b7", 16);
BigInteger bi_n = new BigInteger(
"e8e77781f36a7b3188d711c2190b560f205a" +
"52391b3479cdb99fa010745cbeba5f2adc08" +
"e1de6bf38398a0487c4a73610d94ec36f17f" +
"3f46ad75e17bc1adfec99839589f45f95ccc" +
"94cb2a5c500b477eb3323d8cfab0c8458c96" +
"f0147a45d27e45a4d11d54d77684f65d48f1" +
"5fafcc1ba208e71e921b9bd9017c16a5231af7f", 16);
Console.WriteLine("e =\n" + bi_e.ToString(10));
Console.WriteLine("\nd =\n" + bi_d.ToString(10));
Console.WriteLine("\nn =\n" + bi_n.ToString(10) + "\n");
BigInteger bi_data = new BigInteger(
"12345678901234567890", 10);
BigInteger bi_encrypted = bi_data.modPow(bi_e, bi_n);
BigInteger bi_decrypted = bi_encrypted.modPow(
bi_d, bi_n);
Console.WriteLine("bi_data = " + bi_data);
Console.WriteLine("\nbi_encrypted =\n" + bi_encrypted);
Console.WriteLine("\nbi_decrypted = " + bi_decrypted);
}

## Performance Evaluation

The performance of basic arithmetic operators were evaluated and presented in this section.
We observe significant performance improvements by moving from byte-level to 32-bit word-level operations.
The use of Barrett reduction in modular exponentiation, and the use of a faster algorithm for the
computation of the Jacobi symbol has contributed significantly to improved speed in RSA and
probabilistic primality testing.

All evaluation code was compiled using

`csc BigInteger.cs /o`

and executed on a Pentium III
800Mhz laptop. For tests that involve the generation of random numbers, we use the built-in

`Random`

class and the same initial seed value of

*1* for each test. Timing of the
operations were done using the following code.

int dwStart = System.Environment.TickCount;
Console.WriteLine(System.Environment.TickCount - dwStart);

### Basic Arithmetic Operators

To test the performance of the addition operator, we generated two random `BigIntegers`

of
1024-bits each and added them. This test was repeated 100000 times and the total time taken to complete
the task was recorded. The subtraction and multiplication operators were tested in the similar manner.
For the division operator, we generated two random `BigIntegers`

, each having random length
between 2 to 1024 bits. The larger of the two is then divided by the smaller value. This test was
repeated 100000 times and the total time taken to complete the test was recorded. The results are
shown in the following graph.

### Modular Exponentiation

The speed of modular exponentiation was evaluated by generating three random `BigIntegers`

with
*a = 512 bits*, *e = 512 bits* and *n = 1024 bits*. The exponentiation is then computed as
*a*^{e} mod n. The test was repeated with 100 different values of *a*, *e*, *n*,
and the average time required for each exponentiation was calculated. The following graph shows the
result.

### Jacobi Symbol Computation

The speed of computing the Jacobi symbol was evaluated by generating two random `BigIntegers`

of 1024 bits each. The value of the Jacobi symbol was computed for the two numbers, and the test was
repeated 100 times using different pairs of numbers. The following graph shows the average time required
for each computation.

### Speed of Primality Testing

In our `BigInteger`

class, we provided two functions that performs probabilitic primality
testing of an integer. The first method `bool isProbablePrime(int confidence)`

, uses Rabin-Miller's
strong pseudoprime test to determine whether the integer is probably prime. The algorithm is described
as follows.

- Test the number for divisibility by prime numbers below 2000.
- Perform strong pseudoprime test on the number using
*n* random bases. The number of random bases
used in the test depends on the *confidence* parameter. - The function returns
`true`

if the number passes both step 1 and 2.

The second primality test method `bool isProbablePrime()`

does not require the specification of
a confidence parameter. A brief description of the test is given as follows.

- Test the number for divisibility by prime numbers below 2000.
- Test whether the number is a strong pseudoprime to base 2. If the number passes this test, proceed on to
the next step.
- Test whether the number is a strong Lucas pseudoprime.
- The function returns
`true`

if the number passes steps 1, 2 and 3.

This test works based on the assumption that it is extremely rare for a composite number to be both
a base 2 strong pseudoprime and a strong Lucas pseudoprime. For a detailed discussion of this method,
the reader is referred to [6],[7] and [8].

We compare the performance of the two algorithms by measuring the time it takes for each algorithm
to declare a known prime number as probably prime. The result for 512 and 1024-bit numbers is
given in the following graph.

From the results, it is evident that for a 512-bit prime number,
`isProbablePrime()`

is marginally slower than
`isProbablePrime(10)`

, which uses 10 rounds of Rabin-Miller test.
However, for a 1024-bit prime number, `isProbablePrime()`

is significantly faster
than `isProbablePrime(10)`

.

The following positive pseudoprime and several others passes our implementation of primality test but
failed in JDK's `isProbablePrime`

test. Any advice on why this is so would be greatly
appreciated.

byte[] pseudoPrime1 = { (byte)0x00,
(byte)0x85, (byte)0x84, (byte)0x64, (byte)0xFD,
(byte)0x70, (byte)0x6A, (byte)0x9F, (byte)0xF0,
(byte)0x94, (byte)0x0C, (byte)0x3E, (byte)0x2C,
(byte)0x74, (byte)0x34, (byte)0x05, (byte)0xC9,
(byte)0x55, (byte)0xB3, (byte)0x85, (byte)0x32,
(byte)0x98, (byte)0x71, (byte)0xF9, (byte)0x41,
(byte)0x21, (byte)0x5F, (byte)0x02, (byte)0x9E,
(byte)0xEA, (byte)0x56, (byte)0x8D, (byte)0x8C,
(byte)0x44, (byte)0xCC, (byte)0xEE, (byte)0xEE,
(byte)0x3D, (byte)0x2C, (byte)0x9D, (byte)0x2C,
(byte)0x12, (byte)0x41, (byte)0x1E, (byte)0xF1,
(byte)0xC5, (byte)0x32, (byte)0xC3, (byte)0xAA,
(byte)0x31, (byte)0x4A, (byte)0x52, (byte)0xD8,
(byte)0xE8, (byte)0xAF, (byte)0x42, (byte)0xF4,
(byte)0x72, (byte)0xA1, (byte)0x2A, (byte)0x0D,
(byte)0x97, (byte)0xB1, (byte)0x31, (byte)0xB3,
};
byte[] pseudoPrime2 = { (byte)0x00,
(byte)0x99, (byte)0x98, (byte)0xCA, (byte)0xB8,
(byte)0x5E, (byte)0xD7, (byte)0xE5, (byte)0xDC,
(byte)0x28, (byte)0x5C, (byte)0x6F, (byte)0x0E,
(byte)0x15, (byte)0x09, (byte)0x59, (byte)0x6E,
(byte)0x84, (byte)0xF3, (byte)0x81, (byte)0xCD,
(byte)0xDE, (byte)0x42, (byte)0xDC, (byte)0x93,
(byte)0xC2, (byte)0x7A, (byte)0x62, (byte)0xAC,
(byte)0x6C, (byte)0xAF, (byte)0xDE, (byte)0x74,
(byte)0xE3, (byte)0xCB, (byte)0x60, (byte)0x20,
(byte)0x38, (byte)0x9C, (byte)0x21, (byte)0xC3,
(byte)0xDC, (byte)0xC8, (byte)0xA2, (byte)0x4D,
(byte)0xC6, (byte)0x2A, (byte)0x35, (byte)0x7F,
(byte)0xF3, (byte)0xA9, (byte)0xE8, (byte)0x1D,
(byte)0x7B, (byte)0x2C, (byte)0x78, (byte)0xFA,
(byte)0xB8, (byte)0x02, (byte)0x55, (byte)0x80,
(byte)0x9B, (byte)0xC2, (byte)0xA5, (byte)0xCB,
};

## Future Work

- Faster implementation of arithmetic operations.
- More robust primality testing methods.
- Faster pseudoprime generation

## Conclusion

In this article, we have provided a short introduction to the topic of
large integer arithmetic. We have looked at how large integer addition, subtraction and
multiplication can be implemented. We have also examined the problem
of primality testing and introduced the concept of primality testing based
on Fermat's Little Theorem. Our implementation of `BigInteger`

class
can be downloaded from this page and provides the overloading of most
arithmetic operators. We have pointed out the limitations of our implementations of primality
testing and we are working towards more robust primality testing methods and faster
implementation of arithmetic operators.

## Change Log

- September 23, 2002
- Version 1.03
- Fixed
`operator-`

to give correct data length. - Added Lucas sequence generation.
- Added Strong Lucas Primality test.
- Added integer square root method.
- Added
`setBit`

/`unsetBit`

methods. - New
`isProbablePrime()`

method which do not require the
confidence parameter.

- August 29, 2002
- Version 1.02
- Fixed bug in the exponentiation of negative numbers.
- Faster modular exponentiation using Barrett reduction.
- Added
`getBytes()`

method. - Fixed bug in
`ToHexString`

method. - Added overloading of ^ operator.
- Faster computation of Jacobi symbol.

- August 19, 2002
- Version 1.01
- Big integer is stored and manipulated as unsigned integers (4 bytes)
instead of individual bytes. This gives significant performance improvement.
- Updated Fermat's Little Theorem test to use a
^{(p-1)} mod p =
1 - Added
`isProbablePrime`

method. - Updated documentation.

- August 9, 2002
- Version 1.0
- Initial Release.

[1] D. E. Knuth, "Seminumerical Algorithms", The Art of Computer
Programming Vol. 2,
3rd Edition, Addison-Wesley, 1998.

[2] K. H. Rosen, "Elementary Number Theory and Its Applications", 3rd Ed,
Addison-Wesley, 1993.

[3] B. Schneier, "Applied Cryptography", 2nd Ed, John Wiley & Sons,
1996.

[4] A. Menezes, P. van Oorschot, and S. Vanstone, "Handbook of Applied Cryptography",
CRC Press, 1996, www.cacr.math.uwaterloo.ca/hac

[5] A. Bosselaers, R. Govaerts, and J. Vandewalle, "Comparison of Three Modular
Reduction Functions", Proc. CRYPTO'93, pp.175-186.

[6] R. Baillie and S. S. Wagstaff Jr, "Lucas Pseudoprimes", Mathematics of Computation,
Vol. 35, No. 152, Oct 1980, pp. 1391-1417.

[7] H. C. Williams, "Édouard Lucas and Primality Testing", Canadian Mathematical
Society Series of Monographs and Advance Texts, vol. 22, John Wiley & Sons, New York,
NY, 1998.

[8] P. Ribenboim, "The new book of prime number records", 3rd edition, Springer-Verlag,
New York, NY, 1995.

[9] M. Joye and J.-J. Quisquater, "Efficient computation of full Lucas sequences",
Electronics Letters, 32(6), 1996, pp 537-538.