BigInteger Library

User 6918454

4.77/5 (18 votes)

Feb 22, 2010

CPOL

5 min read

98598

5603

A .NET 2.0 and Mono library for the 64 bit optimized handling of very large integers, up to 10240 binary digits or approximately (safe to use) 3000 decimal digits

Introduction

BigInteger is a .NET 2.0 and Mono library for handling very large integers, up to 10240 binary digits or approximately (safe to use) 3000 decimal digits.

Background

I required a custom type (class) to manage variable-length integers for modular arithmetic and cryptographic operations for my B.Sc. thesis in cryptography. In 2007, I began developing the type and it functioned very well for its' intended purpose. At the beginning of 2010, I started a thorough refactoring of the project, which resulted in this article.

Numeration Base

The first decision relates to the numeration base (radix) to be used. The naive solution involves employing the decimal base, due to familiarity. In the context of computer architecture, this solution is inappropriate because:

it wastes the 32 bit (or 64 bit) integer processor registers by using only 3-4 of the available bits
there is no formulaic correspondence between the number of binary digits and the number of decimal digits for the same number. One only knows the approximation 2¹⁰ = 10³. Therefore, converting binary data into the decimal format will provide inconsistent data sizes.

The decision to use the base 2 radix is almost as unsuitable, since, in this case, although data can be seamlessly imported as numbers, the waste in terms of processor registers' bits is tremendous. The proper option is to select a poly-binary radix, meaning a numeration base that is a power of 2.

Considering the potential overflow of arithmetic operations, it becomes clear that one cannot hope to use a full processor register for storing a digit, without twisting the results of the calculations. The aforementioned reason is coupled with the necessity to utilize a digit size that is divisible by a byte unit (for I/O).

After studying the overflow of the most demanding arithmetic operation algorithms (division), as well as the requirement to accommodate large size integers, made me settle on the 2¹⁶ = 65536 radix, representable on 2 bytes and stored in an 8-byte type to manage in-place the heavy overflow of division on very big numbers.

Implementation-wise, the class' fields look like this:

/// <summary>
/// The number's sign, where Positive also stands for the number zero.
/// </summary>
private enum Sign { Positive, Negative };

/// <summary>
/// 2^16 numeration base for internal computations, in order to benefit the most from the
/// 32 bit (or 64 bit) integer processor registers.
/// </summary>
private const long NumberBase = 65536;

/// <summary>
/// Maximum size for numbers is up to 10240 binary digits or 
/// approximately (safe to use) 3000 decimal digits.
/// The maximum size is, in fact, double the previously specified amount, 
/// in order to accommodate operations'
/// overflow.
/// </summary>
internal const int MaxSize = 2 * 640;

/// <summary>
/// Ratio for the conversion of a BigInteger's size to a binary digits size.
/// </summary>
private const int RatioToBinaryDigits = 16;

/// Integer constants
private static readonly BigInteger Zero = new BigInteger();
private static readonly BigInteger One = new BigInteger(1);
private static readonly BigInteger Two = new BigInteger(2);
private static readonly BigInteger Ten = new BigInteger(10);

/// <summary>
/// The array of digits of the number.
/// </summary>
private long[] digits;

/// <summary>
/// The actual number of digits of the number.
/// </summary>
private int size;

/// <summary>
/// The number sign.
/// </summary>
private Sign sign;

Constructors

The following types of BigInteger constructors have been implemented:

Default constructor - sets the number to 0
Constructor taking a long (64 bit) base-10 number as parameter
Constructor accepting a variable-length string, containing base-10 digits
Constructor using a byte array parameter, for converting raw binary data into a BigInteger, at the ratio 2 raw bytes to 1 BigInteger digit

Addition, Subtraction and Multiplication

The basic arithmetic operations of addition, subtraction and multiplication of two large integers can be implemented easily as a natural extension to the methods used in primary school.

Addition Helper Method Code

/// <summary>
/// Adds two BigNumbers a and b, where a >= b, a, b non-negative.
/// </summary>
private static BigInteger Add(BigInteger a, BigInteger b)
{
     BigInteger res = new BigInteger(a);
     long trans = 0, temp;
     int i;

     for (i = 0; i < b.size; i++)
     {
          temp = res.digits[i] + b.digits[i] + trans;
          res.digits[i] = temp % NumberBase;
          trans = temp / NumberBase;
     }

     for (i = b.size; ((i < a.size) && (trans > 0)); i++)
     {
          temp = res.digits[i] + trans;
          res.digits[i] = temp % NumberBase;
          trans = temp / NumberBase;
     }

     if (trans > 0)
     {
          res.digits[res.size] = trans % NumberBase;
          res.size++;
          trans /= NumberBase;
     }

     return res;
}

Subtraction Helper Method Code

/// <summary>
/// Subtracts the BigInteger b from the BigInteger a, where a >= b, a, b non-negative.
/// </summary>
private static BigInteger Subtract(BigInteger a, BigInteger b)
{
     BigInteger res = new BigInteger(a);
     int i;
     long temp, trans = 0;
     bool reducible = true;

     for (i = 0; i < b.size; i++)
     {
          temp = res.digits[i] - b.digits[i] - trans;
          if (temp < 0)
          {
               trans = 1;
               temp += NumberBase;
          }
          else trans = 0;
          res.digits[i] = temp;
     }

     for (i = b.size; ((i < a.size) && (trans > 0)); i++)
     {
          temp = res.digits[i] - trans;
          if (temp < 0)
          {
               trans = 1;
               temp += NumberBase;
          }
          else trans = 0;
          res.digits[i] = temp;
     }

     while ((res.size - 1 > 0) && (reducible == true))
     {
          if (res.digits[res.size - 1] == 0)
          res.size--;
          else reducible = false;
     }

     return res;
}

Multiplication Helper Method Code

/// <summary>
/// Multiplies two BigIntegers.
/// </summary>
private static BigInteger Multiply(BigInteger a, BigInteger b)
{
     int i, j;
     long temp, trans = 0;

     BigInteger res = new BigInteger();
     res.size = a.size + b.size - 1;
     for (i = 0; i < res.size + 1; i++)
          res.digits[i] = 0;

     for (i = 0; i < a.size; i++)
          if (a.digits[i] != 0)
               for (j = 0; j < b.size; j++)
                    if (b.digits[j] != 0)
                          res.digits[i + j] += a.digits[i] * b.digits[j];

     for (i = 0; i < res.size; i++)
     {
          temp = res.digits[i] + trans;
          res.digits[i] = temp % NumberBase;
          trans = temp / NumberBase;
     }

     if (trans > 0)
     {
          res.digits[res.size] = trans % NumberBase;
          res.size++;
          trans /= NumberBase;
     }

     return res;
}

Division

In opposition to the simplicity in implementing addition, subtraction and multiplication, the division algorithm stands out not only as difficult to implement, but also as one for which good documentation is very hard to procure. After dabbling into Knuth's Seminumerical Algorithms Vol. 2 lacunar explanations, I stumbled upon this wonderful paper by Per Brinch Hansen [1].

There are 5 division helper methods that perform the operation, each of them corresponding to a method in the paper mentioned beforehand.

Classes of Implemented Algorithms

Simple arithmetic operations are not sufficient for a large integer manipulation library. The following algorithms, modular arithmetic functions and primality tests have been implemented, all of which are necessary for modern public-key cryptography algorithms such as the RSA:

Power through fast exponentiation

/// <summary>
/// Returns the power of a BigInteger base to a non-negative exponent by using the
/// fast exponentiation algorithm (right to left binary exponentiation).
/// </summary>
/// <param name="number">The BigInteger base</param>
/// <param name="exponent">The non-negative exponent</param>
/// <returns>The power of the BigInteger base to the non-negative exponent</returns>
/// <exception cref="BigIntegerException">
/// Cannot raise a BigInteger to a negative power exception</exception>
public static BigInteger Power(BigInteger number, int exponent)
{
    if (exponent < 0)
         throw new BigIntegerException
		("Cannot raise an BigInteger to a negative power.", null);

    if (a == Zero)
         return new BigInteger();

    BigInteger res = new BigInteger(1);
    if (exponent == 0)
         return res;

    BigInteger factor = new BigInteger(number);

    while (exponent > 0)
    {
         if (exponent % 2 == 1)
              res *= factor;

         exponent /= 2;

         if (exponent > 0)
              factor *= factor;
    }

    return res;
}

Integer square root (mostly for trying to factor an RSA modulus into its' primes, assuming that the difference between the primes is small enough)

/// <summary>
/// Integer square root of the given BigInteger using Newton's numeric method.
/// </summary>
/// <param name="n">The BigInteger whose integer square root is to be computed
/// </param>
/// <returns>The integer square root of the given BigInteger</returns>
/// <exception cref="BigIntegerException">Cannot compute the 
/// integer square root of a negative number exception</exception>
public static BigInteger IntegerSqrt(BigInteger n)
{
     if (n.sign == Sign.Negative)
          throw new BigIntegerException
	   ("Cannot compute the integer square root of a negative number.", null);

     BigInteger oldValue = new BigInteger(0);
     BigInteger newValue = new BigInteger(n);

     while (BigInteger.Abs(newValue - oldValue) >= One)
     {
          oldValue = newValue;
          newValue = (oldValue + (n / oldValue)) / Two;
     }

     return newValue;
}

Euclid's greatest common divisor algorithm and the generalized version of the GCD (useful for the key generation part of the RSA algorithm)

/// <summary>
/// Euclidean algorithm for computing the 
/// greatest common divisor of two BigIntegers.
/// </summary>
/// <returns>The greatest common divisor of the two given BigIntegers</returns>
/// <exception cref="BigIntegerException">
/// Cannot compute the Gcd of negative BigIntegers exception</exception>
public static BigInteger Gcd(BigInteger a, BigInteger b)
{
     if ((a.sign == Sign.Negative) || (b.sign == Sign.Negative))
          throw new BigIntegerException
		("Cannot compute the Gcd of negative BigIntegers.", null);

     BigInteger r;

     while (b > Zero)
     {
          r = a % b;
          a = b;
          b = r;
     }

     return a;
}

/// <summary>
/// Extended Euclidian Gcd algorithm, returning the 
/// greatest common divisor of two BigIntegers,
/// while also providing u and v, where: a*u + b*v = gcd(a,b).
/// </summary>
/// <param name="u">Output BigInteger parameter, where a*u + b*v = gcd(a,b)
/// </param>
/// <param name="v">Output BigInteger parameter, where a*u + b*v = gcd(a,b)
/// </param>
/// <returns>The greatest common divisor of the two given BigIntegers</returns>
/// <exception cref="BigIntegerException">
/// Cannot compute the Gcd of negative BigIntegers exception</exception>
public static BigInteger ExtendedEuclidGcd(BigInteger a, BigInteger b,
                                           out BigInteger u, out BigInteger v)
{
     if ((a.sign == Sign.Negative) || (b.sign == Sign.Negative))
          throw new BigIntegerException
		("Cannot compute the Gcd of negative BigIntegers.", null);

     BigInteger u1 = new BigInteger();
     BigInteger u2 = new BigInteger(1);
     BigInteger v1 = new BigInteger(1);
     BigInteger v2 = new BigInteger();
     BigInteger q, r;

     u = new BigInteger();
     v = new BigInteger();

     while (b > Zero)
     {
          q = a / b;
          r = a - q * b;
          u = u2 - q * u1;
          v = v2 - q * v1;

          a = new BigInteger(b);
          b = new BigInteger(r);
          u2 = new BigInteger(u1);
          u1 = new BigInteger(u);
          v2 = new BigInteger(v1);
          v1 = new BigInteger(v);
          u = new BigInteger(u2);
          v = new BigInteger(v2);
     }

     return a;
}

Modular arithmetic functions: modular inverse, modular exponentiation by repeated squaring (for the Miller-Rabin primality test and the RSA key generation, encryption and decryption processes)

/// <summary>
/// Computes the modular inverse of a given BigInteger.
/// </summary>
/// <param name="a">The non-zero BigInteger whose inverse is to be computed</param>
/// <param name="n">The BigInteger modulus, 
/// which must be greater than or equal to 2
/// </param>
/// <returns>The BigInteger equal to a^(-1) mod n</returns>
/// <exception cref="BigIntegerException">Invalid number or modulus exception
/// </exception>
public static BigInteger ModularInverse(BigInteger a, BigInteger n)
{
     if (n < Two)
          throw new BigIntegerException("Cannot perform a 
		modulo operation against a BigInteger less than 2.", null);

     if (Abs(a) >= n)
          a %= n;
     if (a.sign == Sign.Negative)
          a += n;

     if (a == Zero)
          throw new BigIntegerException("Cannot obtain the 
				modular inverse of 0.", null);

     if (Gcd(a, n) != One)
          throw new BigIntegerException("Cannot obtain the 
	  modular inverse of a number that is not 
		coprime with the modulus.", null);

     BigInteger u, v;
     ExtendedEuclid(n, a, out u, out v);
     if (v.sign == Sign.Negative)
          v += n;

     return v;
}

/// <summary>
/// Returns the power of a BigInteger to a non-negative exponent modulo n, 
/// by using the fast exponentiation algorithm 
/// (right to left binary exponentiation) and modulo optimizations.
/// </summary>
/// <param name="a">The BigInteger base</param>
/// <param name="exponent">The non-negative exponent</param>
/// <param name="n">The modulus, which must be greater than or equal to 2</param>
/// <returns>The power of the BigInteger to the non-negative exponent</returns>
/// <exception cref="BigIntegerException">Invalid exponent or modulus exception
/// </exception>
public static BigInteger ModularExponentiation
	(BigInteger a, BigInteger exponent, BigInteger n)
{
     if (exponent < 0)
          throw new BigIntegerException
		("Cannot raise a BigInteger to a negative power.", null);

     if (n < Two)
          throw new BigIntegerException("Cannot perform a 
		modulo operation against a BigInteger less than 2.", null);

     if (Abs(a) >= n)
          a %= n;
     if (a.sign == Sign.Negative)
          a += n;

     if (a == Zero)
          return new BigInteger();

     BigInteger res = new BigInteger(1);
     BigInteger factor = new BigInteger(a);

     while (exponent > Zero)
     {
          if (exponent % Two == One)
               res = (res * factor) % n;
          exponent /= Two;
          factor = (factor * factor) % n;
     }

     return res;
}

Composite primality testing method (involving trial division against the primes less than 1000 and then the Miller-Rabin test)

/// <summary>
/// Determines whether the given BigInteger number is probably prime, 
/// with a probability of at least 1/(4^MaxMillerRabinIterations), 
/// using the Miller-Rabin primality test.
/// </summary>
private static bool IsPrimeMillerRabinTest(BigInteger number)
{
     BigInteger s = new BigInteger();
     BigInteger t = number - One;
     BigInteger b = new BigInteger(2);
     BigInteger nmin1 = number - One;
     BigInteger r, j, smin1;

     if (number == One)
          return false;
     else if (number == Two)
          return true;
     else if (number == Three)
          return true;

     while (t % Two == Zero)
     {
          t /= Two;
          s++;
     }

     smin1 = s - One;

     for (int i = 0; i < MaxMillerRabinIterations; i++)
     {
          r = BigInteger.ModularExponentiation(b, t, number);

          if ((r != One) && (r != nmin1))
          {
               j = new BigInteger(One);
               while ((j <= smin1) && (r != nmin1))
               {
                    r = (r * r) % number;
                    if (r == One)
                         return false;
                    j++;
               }
               if (r != nmin1)
                    return false;
          }

          if (b == Two)
               b++;
          else
               b += Two;
      }

      return true;
}

The theory and algorithms for most of the previously mentioned operations can be found in my B.Sc. thesis [2], chapters 2 and 3.

Additional Mentions

To display a BigInteger in the base 10 (decimal system), the library uses a toned-down version of the BigInteger class called Base10BigInteger, where the numeration base is set to 10 and the number converted on-the-fly from the base 65536 of the BigInteger class. The availability in C# of a C++-style generics abstraction allowing parameterization by value (radix) would have removed the necessity of this helper class.

The BigInteger library is made to be simple and straightforward to use. It overloads every arithmetic operator, so a client programmer can work with BigIntegers as with ordinary built-in numeric types.

Also, the BigInteger class implements the interfaces ISerializable, IEquatable<T>, IComparable and IComparable<T>.

Sample Code

Here is a small sample of interacting with the library:

BigInteger m = new BigInteger("982451159");
BigInteger m2 = m + 31567;
Console.WriteLine(BigInteger.IsPrime(m));
BigInteger n = BigInteger.Power(m, 10);
BigInteger p = BigInteger.Power(m, 4);
Console.WriteLine(BigInteger.IntegerSqrt(n) / p == m);
BigInteger q = BigInteger.ModularExponentiation(m, n, p);
Console.WriteLine(q.ToString());

Library and Cryptography Applications

The complete source code and binaries of the BigInteger and the BigInteger-Mono libraries (belonging to the same Google Code project) can be accessed here.

I also made available a set of cryptography applications implementing the RSA and the Rabin cryptosystems, respectively, that rely on an older version of the BigInteger library. These can be found here.

I would gladly welcome contributions and feedback to this BigInteger open-source (LGPL) project.

References

[1] P. B. Hansen, Multiple-length Division Revisited: a Tour of the Minefield, Software - Practice and Experience, vol. 24(6), 579-601, 1994

[2] The Prime Pages, http://primes.utm.edu/

[3] The Microsoft Developer Network (MSDN) pages

History

Version 0.1 - Initial submission - 20/02/2010
Version 0.2 - Added the project content to CodeProject - 21/02/2010
Version 0.5 - Project content and explanations update - 23/02/2010
Version 1.0 - Source code revisions - 24/02/2010
Version 1.1 - Added download links at the head of the article - 28/02/2010
Version 1.2 - Added full XML comments to the public methods inside the code - 13/03/2010
Version 1.3 - Added binaries for Mono - 04/06/2010
Version 1.4 - Updated content, sources and binaries - 21/09/2011