The bignum Library

This library is an adaptation of the “miracl” bignum library published by Mike Scott in 1989. Since then, Mike has developed his library further. The latest version can be downloaded from his web site

ftp.compapp.dcu.ie/pub/crypto/miracl.zip

I had worked a lot in the 1989 version, rewriting the essential parts of it in pure assembly. It runs quite fast now, and I am convinced it is a very useful package.

The interface this library presents is described in “bignums.h”. Basically, most operators (addition, subtraction, etc.) are overloaded there, so that you can use the library using a practical infix notation (you write c = a + b instead of quadsum(a,b,c) to make an addition).

Practical considerations:

Note that the bignum library is designed to work with the garbage collector. Most operators that need to, for instance “+”, “-“ etc., return a new bignum, that will be garbage collected automatically if it is no longer needed. This can be costly, and if you prefer trading speed for notational convenience, you can use the normal C interface and instead of writing:

c = (a + b)/(a - b);//creates temporary objects for a+b and a-b

you write:

     c = newBignum(0);  //allocate result

     quadsum(a,b,tmp);  // tmp = a + b

     quadsub(a,b,c);    // c = a - b

     quaddiv(tmp,c,c);  // c = tmp/c = (a+b)/(a-b)

In this latter case, you create only a single temporary that can later be reused. In the former case, you create two. This is absolutely unnoticeable in MOST cases. The convenience of the infix notation largely compensates the occasional stops to collect unneeded structures. Of course, you can trade execution speed for development speed.

Note that the public interface of the bignum library is an opaque data structure:

typedef struct _mp {

        void *mp;

} pBignum;

This allows for changes in the bignum representation to be made transparently to your code. You could plug in a completely different bignum library, or get the latest version of Mike Scott and plug it in, without any changes to the bignum client code. Note that this is just a glorified pointer. Before using it, you must initialize it to a new bignum, either as a result of an operation or as the result of newBignum().

The library is presented as a DLL. When loaded, the DLL will automatically initialize the bignum system to a precision of 300 bits. You can change this using the BignumPrecision() API.   

Example:

This example prints the powers of two from 127 to zero.

#include <bignums.h>

int main(void)

{

        pBignum a,b,c;

        char buffer[4096];

        c = 2; // Note the overloaded assignment

        for (int i = 127;i> 0;i--) {

                b = newBignum(i);

                a = quadpow(c,b);

                quadformat(a,buffer);

                printf("[%3d] %s\n",i,buffer);

        }

        return 0;

}

Output:

[127] 170141183460469231731687303715884105728

[126] 85070591730234615865843651857942052864

[125] 42535295865117307932921825928971026432

[124] 21267647932558653966460912964485513216

[123] 10633823966279326983230456482242756608

[122] 5316911983139663491615228241121378304

[121] 2658455991569831745807614120560689152

[120] 1329227995784915872903807060280344576

[119] 664613997892457936451903530140172288

…

Interface overview

Debugger support.

Bignums are recognized in the latest versions of wedit and the debugger shows them as numbers.

Important considerations when writing bignum code

·        A pBignum is actually a pointer to a number. This is more efficient, but can lead to mistakes. When you write:

pBignum a;

// Some code here

pBignum b = a;

 After this assignment both a and b point to the same number! If you modify b, a will get modified too. To avoid this problem you should write

     pBignum b = quadcopy(a);

This forces b to point to a different area of memory, i.e. a copy of a.  This semantics are obvious for double precision numbers or for integers, but here they must be explicitly enforced by the programmer.

·        Avoid floating point numbers when working with bignums if possible. For instance, do not write:

pBignum a = 0.5; // WRONG

You should write:

pBignum a = quadratio(1,2); //OK

The later expression stores exactly 1/2. The former stores an approximation to 1/2 that will be good in the first 16 digits only.

Bignum implementation

 Bignums are stored in two ways:

·        As a large integer: The first word (32 bits) contains the length and the sign, followed by the data, 30 bits for each 32 bit word.

·        As a ratio of two integers, represented as a numerator and denominator. This allows for a floating point (or maybe floating slash) representation.

In both cases, the library returns a pointer to these opaque data. No direct access is provided, so the underlying library or DLL can be changed as needed, without affecting user code, as long as the interface of the dll is the same.

Here is an excerpt from the "Miracl" reference manual that explains the basic implementation of the library.

Floating-Slash numbers

The straightforward way to represent rational numbers is as reduced fractions, as a numerator and denominator with all common factors cancelled out. These numbers can then be added, subtracted, multiplied and divided in the obvious way and the result reduced by dividing both numerator and denominator by their Greatest Common Divisor. An efficient GCD subroutine, using Lehmers modification of the classical euclidean algorithm for multiprecision numbers [Knuth81], is included in the MIRACL package.

An alternative way to represent rationals would be as a finite continued fraction [Knuth81]. Every rational number p/q can be written as

or more elegantly as p/q = [a₀/a₁/a₂/..../a_n] where the a_i are positive integers, usually quite small.

For example

Note that the a_i elements of the above continued fraction representation are easily found as the quotients generated as a by-product when the euclidean GCD algorithm is applied to p and q.

As we are committed to fixed length representation of rationals, a problem arises when the result of some operation exceeds this fixed length. There is a necessity for some scheme of truncation, or rounding. While there is no obvious way to truncate a large fraction, it is a simple matter to truncate the continued fraction representation. The resulting, smaller, fraction is called a best rational approximation, or a convergent, to the original fraction.

Consider truncating 277/642 = [0/2/3/6/1/3/3]. Simply drop the last element from the CF representation, giving [0/2/3/6/1/3] = 85/197, which is a very close approximation to 277/642 (error = 0.0018%). Chopping more terms from the CF expansion gives the successive convergents as 22/51, 19/44, 3/7, 1/2, 0/1. As the fractions get smaller, the error increases. Obviously the truncation rule for a computer implementation should be to choose the biggest convergent that fits the computer representation.

The type of rounding described above is also called ‘Mediant rounding’. If p/q and r/s are two neighbouring representable slash numbers astride a gap, then their mediant is the unrepresentable (p+r)/(q+s). All larger fractions between p/q and the mediant will round to p/q, and those between r/s and the mediant will round to r/s. The mediant itself rounds to the ‘simpler’ of p/q and r/s.

 This is theoretically a very good way to round, much better than the rather arbitrary and base-dependent methods used in floating-point arithmetic, and is the method used here. The full theoretical basis of floating-slash arithmetic is described in detail by Matula & Kornerup [Matula85]. It should be noted that our flash representation is in fact a cross between the fixed- and floating-slash systems analysed by Matula & Kornerup, as our slash can only float between words, and not between bits. However the characteristics of the flash data-type will tend to those of floating-slash, as the precision is increased.

 The MIRACL routine round implements mediant rounding. If the result of an arithmetic operation is the fraction p/q, then the euclidean GCD algorithm is applied as before to p and q. However this time the objective is not to use the algorithm to calculate the GCD per se, but to use its quotients to build successive convergents to p/q. This process is stopped when the next convergent is too large to fit the flash representation. The complete algorithm is given below (Kornerup & Matula [Korn83])

            Given p³0 and q³1

        b_-2=p                x_-2=0                y_-2=1

        b_-1=q                x_-1=1                y_-1=0

Now for i=0,1,..... and for b_i-1>0, find the quotient a_i and remainder b_i when b_i-2 is divided by b_i-1, such that

        b_i = -a_i.b_i-1 + b_i-2

Then calculate

        x_i = a_i.x_i-1 + x_i-2

        y_i = a_i.y_i-1 + y_i-2

Stop when  is to big to fit the flash representation, and take  as the rounded result.

If applied to 277/642, this process will give the same sequence of convergents as stated earlier.

Since this rounding procedure must be applied to the result of each arithmetic operation, and since it is potentially rather slow, a lot of effort has been made to optimise its implementation. Lehmer's idea of operating only with the most significant piece of each number for as long as possible [Knuth81] is used, so that for most of the iterations only single-precision arithmetic is needed. Special care is taken to avoid the rounded result overshooting the limits of the flash representation [Scott89a]. The application of the basic arithmetic routines to the calculation of elementary functions such as log(x), exp(x), sin(x), cos(x), tan(x) etc., uses the fast algorithms described by Brent [Brent76].

A disadvantage of using a flash type of variable to approximate real arithmetic is the non-uniformity in gap-size between representable values (Matula & Kornerup [Matula85]).

To illustrate this consider a floating-slash system which is constrained to have the product of numerator and denominator less than 256. Observe that the first representable fraction less than 1/1 in such a system is 15/16, a gap of 1/16. The next fraction larger than 0/1 is 1/255, a gap of 1/255. In general, for a k-bit floating-slash system, the gap size varies from smaller than 2^-k to a worst case 2^-k/2. In practise this means that a real value that falls into one of the larger gaps, will be represented by a fraction which will be accurate to only half its usual precision. Fortunately such large gaps are rare, and increasingly so for higher precision, occurring only near simple fractions. However it does mean that real results can only be completely trusted to half the given decimal places. A partial solution to this problem would be to represent rationals directly as continued fractions. This gives a much better uniformity of gap-size (Kornerup & Matula [Korn85]), but would be very difficult to implement using a high level language.

Arithmetic on flash data-types is undoubtedly slower than on an equivalent sized multiprecision floating-point type (e.g. [Brent78]). The advantages of the flash approach are its ability to exactly represent rational numbers, and do exact arithmetic on them. Even when rounding is needed, the result often works out correctly, due to the tendency of mediant-rounding to prefer a simple fraction over a complex one. For example the roots program (Chapter 8) when asked to find the square root of 2 and then square the result, comes back with the exact answer of 2, despite much internal rounding.