C# BigInteger Class

• Download source files - 26 Kb
• Download demo - 11 Kb

Contents

• Introduction
• Features
• Implementation Details
• Addition of Big Integers
• Subtraction of Big Integers
• Multiplication of Big Integers
• Division of Big Integers
• Primality Testing using Fermat's Little Theorem
• Sample Code
• Performance Evaluation
• Basic Arithmetic Operators
• Modular Exponentiation
• Jacobi Symbol Computation
• Speed of Primality Testing
• Limitations<
• Future Work
• Conclusion
• Change Log
• References

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₁₀= 1x10⁴ + 2x10³ + 3x10² + 4x10¹ + 5x10⁰

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

12345₁₆= 1x16⁴ + 2x16³ + 3x16² + 4x16¹ + 5x16⁰

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_kb^k + a_k-1b^k-1 + ... + a₁b¹ + a₀

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_kb^k + a_k-1b^k-1 + ... + a₁b¹ + a₀
m = c_kb^k + c_k-1b^k-1 + ... + c₁b¹ + c₀

Then
n + m = (a_k + c_k)b^k + (a_k-1 + c_k-1)b^k-1 + ... + (a₁ + c₁)b¹ + (a₀ + c₀)

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.

// code segment that performs byte-by-byte addition

int carry = 0;
for(int i = 0; i < len; i++)
{
        // bi1 and bi2 are the two large numbers to add.

        // data is the byte array holding the large number.


        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.

// code segment that performs byte-by-byte subtraction


int carryIn = 0;

for(int i = 0; i < len; i++)
{
    int diff, totalSub;

    // calculate the total value to be subtracted

    // from the current

    // digit.

    totalSub = bi2.data[i] + carryIn;

    // if total value to be subtracted is larger than the

    // currentdigit.


    if(bi1.data[i] < totalSub)
    {
        // we add the value of the base to the current

        // digit and remember to subtract 1 from the

        // next higher digit.


        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²). 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²) 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²) byte additions. Thus, the total complexity is O(n²) byte multiplications and O(n²) 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++)
    {
        // compute the byte multiplication and the addition

        // of the partial results in a single step.


        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ⁿ mod n = a

for any integer a and any prime n. This means that if n is prime then aⁿ mod n = a for any integer a. The contrapositive of this statement says that if aⁿ 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ⁿ 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ⁿ 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.

// confidence determines the number of times to test.


for(int round = 0; round < confidence; round++)
{
    // generate random a < n

    a.genRandomBits(testBits, rand);

    // calculate a^(n) mod n, where thisVal = n

    BigInteger expResult = a.modPow(thisVal, thisVal);

    // is NOT prime if a^(p) mod p != a

    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ⁿ 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)
{
    // define two BigIntegers in base 10

    BigInteger bn1 = new BigInteger("12345678901234567890",
        10);
    BigInteger bn2 = new BigInteger("4321", 10);

    // Determine the quotient and remainder by dividing

    // the first number by the second.

    BigInteger bn3 = bn1 / bn2;
    BigInteger bn4 = bn1 % bn2;

    // Recalculate the number

    BigInteger bn5 = (bn3 * bn2) + bn4;

    Console.WriteLine("bn1 = " + bn1);
    Console.WriteLine("bn5 = " + bn5);
}

Asymmetrical encryption and decryption.

public static void Main(string[] args)
{
    // private and public key

    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");

    // data

    BigInteger bi_data = new BigInteger(
        "12345678901234567890", 10);

    // encrypt and decrypt data

    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;

    // Perform the test.


    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.

1. Test the number for divisibility by prime numbers below 2000.
2. 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.
3. 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.

1. Test the number for divisibility by prime numbers below 2000.
2. Test whether the number is a strong pseudoprime to base 2. If the number passes this test, proceed on to the next step.
3. Test whether the number is a strong Lucas pseudoprime.
4. 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).

Limitations

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.

References

[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.