Arbitrary-precision integer arithmetic (C)

From LiteratePrograms
Jump to: navigation, search

In many applications in number theory and other areas, it's convenient to be able to perform simple arithmetic operations (addition, subtraction, multiplication, division, remainder) on arbitrarily large integers. This article describes some very simple algorithms for computing these in C.

We will represent our datatype with a structure integer that encodes large integers as arrays of unsigned 16-bit integers, with the least significant component first:

<<type definitions>>=
typedef unsigned short component_t;
typedef unsigned long double_component_t;

<<integer type>>=
typedef struct {
    component_t* c;    /* least-significant word first */
    int num_components;
} integer;

<<prototypes>>=
integer create_integer(int components);
<<integer constructor>>=
integer create_integer(int components) {
    integer result;
    result.num_components = components;
    result.c = (component_t*)malloc(sizeof(component_t)*components);
    return result;
}

<<prototypes>>=
void free_integer(integer i);
<<integer destructor>>=
void free_integer(integer i) {
    free(i.c);
}

<<header files>>=
#include <stdlib.h>   /* for malloc() */

The precise definition of a component_t may vary from platform to platform, but a short is guaranteed to be at least 16 bits. Why do we use 16-bit components? Because to perform our operations efficiently we need a datatype double_component_t twice as wide as the component datatype, and C does not guarantee that any 64-bit integer datatype is supplied. On platforms supporting these, it is only necessary to modify the typedefs.

We'll find this small function helpful, which zeros out an integer using memset():

<<prototypes>>=
void set_zero_integer(integer i);
<<integer set to zero>>=
void set_zero_integer(integer i) {
    memset(i.c, 0, sizeof(component_t)*i.num_components);
}

<<header files>>=
#include <string.h>   /* for memset(), memmove() */

Similar is the function for copying an integer, which has the slight complication that it must either truncate or zero-extend if the lengths differ:

<<prototypes>>=
void copy_integer(integer source, integer target);
<<integer copy>>=
void copy_integer(integer source, integer target) {
    memmove(target.c, source.c,
            sizeof(component_t)*MIN(source.num_components, target.num_components));
    if (target.num_components > source.num_components) {
        memset(target.c + source.num_components, 0,
               sizeof(component_t)*(target.num_components - source.num_components));
    }
}

Contents

[edit] Addition

We can now define an addition function; this is the same addition used in grade-school methods, except that we're working with digits in base 216 (the components of the number). Keep in mind that the least significant component is at index zero. To avoid unnecessary costly allocations, the structure for the result is passed in to the method, and must have at least one more component than the larger of the two inputs to prevent overflow:

<<prototypes>>=
void add_integer(integer left, integer right, integer result);
<<integer addition>>=
void add_integer(integer left, integer right, integer result) {
    double_component_t carry = 0;
    int i;
    for(i=0; i<left.num_components || i<right.num_components || carry != 0; i++) {
        double_component_t partial_sum = carry;
        carry = 0;
        if (i < left.num_components)  partial_sum += left.c[i];
        if (i < right.num_components) partial_sum += right.c[i];
        if (partial_sum > MAX_COMPONENT) {
            partial_sum &= MAX_COMPONENT;
            carry = 1;
        }
        result.c[i] = (component_t)partial_sum;
    }
    for ( ; i < result.num_components; i++) { result.c[i] = 0; }
}

We must be careful to account for the case where either number is shorter, or where there is a carry beyond the end of both input numbers. MAX_COMPONENT is the maximum value of a component, which we define using a trick with −1 that works for unsigned types:

<<constants>>=
#define MAX_COMPONENT  ((component_t)(-1))

The above implementation also has the useful property that, because we never read index i of either input after writing index i of the result, we can specify the same integer for the result as for one of the inputs. Subtraction and multiplication by a small number will also have this property, but full multiplication will not.

[edit] Subtraction

For subtraction we again use the grade school method, but this is a bit trickier because borrowing is involved. If we find ourselves trying to subtract a larger component from a smaller component, we have to borrow:

<<prototypes>>=
void subtract_integer(integer left, integer right, integer result);
<<integer subtraction>>=
void subtract_integer(integer left, integer right, integer result) {
    int borrow = 0;
    int i;
    for(i=0; i<left.num_components; i++) {
        double_component_t lhs = left.c[i];
        double_component_t rhs = (i < right.num_components) ? right.c[i] : 0;
        deal with borrowing from previous step
        if (lhs < rhs) {
            borrow = 1;
            lhs += MAX_COMPONENT + 1;
        }
        result.c[i] = lhs - rhs;
    }
    for ( ; i < result.num_components; i++) { result.c[i] = 0; }
}

When we borrow, we just set a flag. When we reach the next step, we try to subtract one from partial_diff for the borrow. If left.c[i] ≤ right.c[i] however, this will require another borrow from the next step up:

<<deal with borrowing from previous step>>=
if (borrow) {
    if (lhs <= rhs) {
        /* leave borrow set to 1 */
        lhs += (MAX_COMPONENT + 1) - 1;
    } else {
        borrow = 0;
        lhs--;
    }
}

[edit] Multiplication by a small number

Multiplying by a small number, no larger than one component, is straightforward. We multiply by each component, set the result to the low COMPONENT_BITS bits, and add the higher bits to the next component as carry:

<<prototypes>>=
void multiply_small_integer(integer left, component_t right, integer result);
<<integer multiplication by a small number>>=
void multiply_small_integer(integer left, component_t right, integer result) {
    double_component_t carry = 0;
    int i;
    for(i=0; i<left.num_components || carry != 0; i++) {
        double_component_t partial_sum = carry;
        carry = 0;
        if (i < left.num_components)  partial_sum += left.c[i]*right;
        carry = partial_sum >> COMPONENT_BITS;
        result.c[i] = (component_t)(partial_sum & MAX_COMPONENT);
    }
    for ( ; i < result.num_components; i++) { result.c[i] = 0; }
}

[edit] Multiplication

The multiplication algorithm is again based on grade-school algorithms, and is similar to addition. As in addition, we compute the columns of the result right-to-left, and any carry from a column is added to the next column to its left. However, the numbers we're adding are different: to get these numbers, we must multiply components of the input number whose indexes add to a fixed value. For example, the component result.c[2] is given by the low 16 bits of this expression:

left.c[0]*right.c[2] + left.c[1]*right.c[1] + left.c[2]*right.c[0] + the carry from result.c[1]

Notice that in all these products, the indexes into left and right add to 2.

The code works by first fixing the result column index and then iterating over all possible pairs of indexes into left and right that sum to that index - the bounds are a bit tricky and depend on the sizes of left and right. Also, because the multiplications above produce double-size values, we must separate out the upper half often to prevent overflow.

<<prototypes>>=
void multiply_integer(integer left, integer right, integer result);
<<integer multiplication>>=
void multiply_integer(integer left, integer right, integer result) {
    int i, lidx, ridx;
    double_component_t carry = 0;
    int max_size_no_carry;
    int left_max_component  = left.num_components - 1;
    int right_max_component = right.num_components - 1;
    minimize max component indexes
    max_size_no_carry = left_max_component + right_max_component;
    for(i=0; i <= max_size_no_carry || carry != 0; i++) {
        double_component_t partial_sum = carry;
        carry = 0;
        lidx = MIN(i, left_max_component);
        ridx = i - lidx;
        while(lidx >= 0 && ridx <= right_max_component) {
            partial_sum += ((double_component_t)left.c[lidx])*right.c[ridx];
            carry += partial_sum >> COMPONENT_BITS;
            partial_sum &= MAX_COMPONENT;
            lidx--; ridx++;
        }
        result.c[i] = partial_sum;
    }
    for ( ; i < result.num_components; i++) { result.c[i] = 0; }
}

Above, COMPONENT_BITS is the number of bits in a component, which can be computed at compile time using sizeof and CHAR_BIT, since C guarantees sizeof(char) == 1:

<<constants>>=
#define COMPONENT_BITS  (sizeof(component_t)*CHAR_BIT)

<<type definition headers>>=
#include <limits.h>  /* for CHAR_BIT */

MIN() is just the obvious macro finding the minimum of two side-effect free expressions, which we use to ensure that lidx is in range at the beginning of the loop. We similarly define MAX, which we'll use later.

<<macros>>=
#define MIN(x,y)  ((x)<(y) ? (x) : (y))
#define MAX(x,y)  ((x)>(y) ? (x) : (y))

Since the time required increases rapidly as we add more components, we want to ensure the max left and right component indexes are as small as possible. We do this by ignoring leading zero components, pointing the indexes at the most significant nonzero component:

<<minimize max component indexes>>=
while(left.c[left_max_component]   == 0) left_max_component--;
while(right.c[right_max_component] == 0) right_max_component--;

[edit] Comparison operator

It's helpful in many algorithms to be able to determine whether an integer is less than, equal to, or greater than another integer. Much like comparing strings, we simply compare the integers lexicographically by their components, starting with the most significant component. Any components beyond the end of an integer are treated as zero.

<<prototypes>>=
int compare_integers(integer left, integer right);
<<integer comparison>>=
int compare_integers(integer left, integer right) {
    int i = MAX(left.num_components - 1, right.num_components - 1);
    for ( ; i >= 0; i--) {
        component_t left_comp =
            (i < left.num_components) ? left.c[i] : 0;
        component_t right_comp =
            (i < right.num_components) ? right.c[i] : 0;
        if (left_comp < right_comp)
            return -1;
        else if (left_comp > right_comp)
            return  1;
    }
    return 0;
}

Like the comparators passed to library routines like qsort(), we return a negative, zero, or positive value depending on whether left is less than, equal to, or greater than right, respectively.

[edit] Shift left/right

Performing a bit shift on a larger integer is handy for rapidly multiplying or dividing it by a power of 2. The easiest operation is shifting left by 1. We shift each component left by one, but also make sure to propagate the high bit of each component into the low bit of the next component up. In this implementation we discard the highest bit.

<<prototypes>>=
void shift_left_one_integer(integer arg);
<<integer shift left one>>=
void shift_left_one_integer(integer arg) {
    int i;
    arg.c[arg.num_components - 1] <<= 1;
    for (i = arg.num_components - 2; i >= 0; i--) {
        arg.c[i + 1] |= arg.c[i] >> (COMPONENT_BITS - 1);
        arg.c[i] <<= 1;
    }
}

Shifting right by one is similar, but propagating from the lowest to the highest bit and moving right to left instead of left to right:

<<prototypes>>=
void shift_right_one_integer(integer arg);
<<integer shift right one>>=
void shift_right_one_integer(integer arg) {
    int i;
    arg.c[0] >>= 1;
    for (i = 1; i < arg.num_components; i++) {
        arg.c[i - 1] |= (arg.c[i] & 1) << (COMPONENT_BITS - 1);
        arg.c[i] >>= 1;
    }
}

[edit] Remainder with a small divisor

Next, we need a function to find the remainder of an integer after division by a small divisor, no larger than one component. The simplest way to do this is to look at the binary representation of the dividend. This expresses it as a sum of powers of 2. For example, for 13 we have:

13 = 11012 = 23 + 22 + 20
13 mod n = ((23 mod n) + (22 mod n) + (20 mod n)) mod n

Here we're using "mod" in the sense of "remainder". As it turns out, computing the sequence 2i mod n is relatively easy:

  • Start with k=1. Multiply k by 2 repeatedly (we can effect this with shift_left_integer).
  • If ever kn, subtract n from k (using subtract_integer).

This works because if k < n, it follows that 2k < 2n, and so we cannot subtract out more than one n. This doesn't work for larger bases. As we produce the sequence, we add together the sum of powers of 2 above. Since each value we add is also less than n, by subtracting n whenever the value exceeds n the result will be less than n.

Here's some C code for this algorithm, relying on operations we've already written.

<<prototypes>>=
component_t mod_small_integer(integer left, component_t right);
<<integer remainder with a small divisor>>=
component_t mod_small_integer(integer left, component_t right) {
    double_component_t mod_two_power = 1;
    double_component_t result = 0;
    int i, bit;
    for(i=0; i<left.num_components; i++) {
        for(bit=0; bit<COMPONENT_BITS; bit++) {
            if ((left.c[i] & (1 << bit)) != 0) {
                result += mod_two_power;
                if (result >= right) {
                    result -= right;
                }
            }
            mod_two_power <<= 1;
            if (mod_two_power >= right) {
                mod_two_power -= right;
            }
        }
    }
    return (component_t)result;
}

[edit] Remainder

Generalizing the above to a remainder by a large integer is straightforward, just replacing some of the component-sized operations with previously implemented arbitrary-precision operations. Note that we assume for convenience that result has one more component than right, the modulus.

<<prototypes>>=
void mod_integer(integer left, integer right, integer result);
<<integer remainder>>=
void mod_integer(integer left, integer right, integer result) {
    integer mod_two_power = create_integer(right.num_components + 1);
    int i, bit;
    set_zero_integer(result);
    set_zero_integer(mod_two_power);
    mod_two_power.c[0] = 1;
    for(i=0; i<left.num_components; i++) {
        for(bit=0; bit<COMPONENT_BITS; bit++) {
            if ((left.c[i] & (1 << bit)) != 0) {
                add_integer(result, mod_two_power, result);
                if (compare_integers(result, right) >= 0) {
                    subtract_integer(result, right, result);
                }
            }
            shift_left_one_integer(mod_two_power);
            if (compare_integers(mod_two_power, right) >= 0) {
                subtract_integer(mod_two_power, right, mod_two_power);
            }
        }
    }
    free_integer(mod_two_power);
}

[edit] Division by a small number

Dividing by a small number, the size of one component, is simple using the grade-school method of division. We simply proceed left-to-right, dividing each component by the divisor. The quotient becomes the corresponding component of the result, and the remainder becomes the high 16 bits of the dividend in the next component's division:

<<prototypes>>=
void divide_small_integer(integer left, component_t right, integer result);
<<integer division by a small number>>=
void divide_small_integer(integer left, component_t right, integer result) {
    double_component_t dividend = 0;
    int i;
    for (i = left.num_components - 1; i >= 0; i--) {
        dividend |= left.c[i];
        result.c[i] = dividend/right;
        dividend = (dividend % right) << COMPONENT_BITS;
    }
}

[edit] Conversion to/from string

Producing an integer from a decimal string is a simple matter of repeatedly alternating between multiplying by 10 and adding the next digit. To determine how many components we need, we use this formula, where b is COMPONENT_BITS:

\log_{2^b} n = (\log_2 10 \cdot \log_{10} n)/b

This follows from basic laws of exponents. Noting that the string length is the ceiling of log10n, this allows us to find a fairly tight upper bound on \log_{2^b} n:

<<prototypes>>=
integer string_to_integer(char* s);
<<string to integer conversion>>=
integer string_to_integer(char* s) {
    integer result = create_integer((int)ceil(LOG_2_10*strlen(s)/COMPONENT_BITS));
    set_zero_integer(result);
    integer digit = create_integer(1);
    int i;
    for (i = 0; s[i] != '\0'; i++) {
        multiply_small_integer(result, 10, result);
        digit.c[0] = s[i] - '0';
        add_integer(result, digit, result);
    }
    free_integer(digit);
    return result;
}

<<header files>>=
#include <math.h>  /* for ceil */

The constant LOG_2_10 is precomputed:

<<constants>>=
#define LOG_2_10    3.3219280948873623478703194294894

To convert in the other direction, we alternate between dividing by 10 and using modulus to extract the next digit, making a special case for zero. To determine how large a string to allocate we rearrange the formula above to obtain:

\log_{10} n = (b\log_{2^b} n)/(\log_2 10)
<<prototypes>>=
char* integer_to_string(integer x);
<<integer to string conversion>>=
char* integer_to_string(integer x) {
    int i, result_len;
    char* result =
        (char*)malloc((int)ceil(COMPONENT_BITS*x.num_components/LOG_2_10) + 2);
    integer ten = create_integer(1);
    ten.c[0] = 10;

    if (is_zero_integer(x)) {
        strcpy(result, "0");
    } else {
        for (i = 0; !is_zero_integer(x); i++) {
            result[i] = (char)mod_small_integer(x, 10) + '0';
            divide_small_integer(x, 10, x);
        }
        result[i] = '\0';
    }

    reverse string

    free_integer(ten);
    return result;
}

This particular implementation destroys the integer x in the process, to avoid copying. Since the digits are extracted in reverse order, we must reverse the string at the end to get the usual printing order:

<<reverse string>>=
result_len = strlen(result);
for(i=0; i < result_len/2; i++) {
    char temp = result[i];
    result[i] = result[result_len - i - 1];
    result[result_len - i - 1] = temp;
}

The helper function is_zero_integer tests whether an integer is zero by comparing each component to zero:

<<prototypes>>=
int is_zero_integer(integer x);
<<test if integer is zero>>=
int is_zero_integer(integer x) {
    int i;
    for(i=0; i < x.num_components; i++) {
        if (x.c[i] != 0) return 0;
    }
    return 1;
}

[edit] Files

<<integer.h>>=
#ifndef _INTEGER_H_
#define _INTEGER_H_

type definition headers
type definitions
constants

macros

integer type

prototypes

#endif /* #ifndef _INTEGER_H_ */
<<integer.c>>=
header files
#include "integer.h"

integer constructor

integer destructor

integer set to zero

test if integer is zero

integer copy

integer addition

integer subtraction

integer multiplication by a small number

integer multiplication

integer comparison

integer shift left one

integer shift right one

integer remainder with a small divisor

integer remainder

integer division by a small number

string to integer conversion

integer to string conversion

[edit] Sample code

This sample code allows a user to interactively use several of the functions above, as well as some new ones built using them. It also incidentally stresses a number of other functions through string conversion, which depends on mod, which in turn depends on several others.

<<integer_sample.c>>=
#include <stdio.h>
#include <string.h>
#include <stdlib.h>  /* for atoi */
#include "integer.h"

int main(int argc, char* argv[]) {
    integer a1, a2, result;
    a1.num_components = a2.num_components = 0; /* quelch compiler warnings */
    if (argc > 2) a1 = string_to_integer(argv[2]);
    if (argc > 3) a2 = string_to_integer(argv[3]);
    if (strcmp(argv[1], "add") == 0) {
        result = create_integer(a1.num_components + a2.num_components + 1);
        add_integer(a1, a2, result);
    } else if (strcmp(argv[1], "subtract") == 0) {
        result = create_integer(a1.num_components + a2.num_components + 1);
        subtract_integer(a1, a2, result);
    } else if (strcmp(argv[1], "multiply") == 0) {
        result = create_integer(a1.num_components + a2.num_components + 1);
        multiply_integer(a1, a2, result);
    } else if (strcmp(argv[1], "mod") == 0) {
        result = create_integer(a2.num_components + 1);
        mod_integer(a1, a2, result);
    } else if (strcmp(argv[1], "factorial") == 0) {
        int i, n = atoi(argv[2]);
        result = create_integer(n*n + 1);
        set_zero_integer(result);
        result.c[0] = 1;
        for(i=2; i<=n; i++) {
            multiply_small_integer(result, i, result);
        }
    }
    puts(integer_to_string(result));
    return 0;
}

Here's some sample input and output (some lines have been split using \ for readability):

$ ./integer_sample add 5000000000000 4000000000000
9000000000000
$ ./integer_sample subtract 5000000000000 4000000000000
1000000000000
$ ./integer_sample multiply 123456789 987654321
121932631112635269
$ ./integer_sample factorial 100
9332621544394415268169923885626670049071596826438162146859296389521759999322991\
5608941463976156518286253697920827223758251185210916864000000000000000000000000

[edit] Potential improvements

Although simple, the algorithms described in this article are too inefficient for very large numbers with millions of digits. There are a number of better multiplication algorithms, the most well known of which runs in O(n log n) time and is based on the Fast Fourier Transform. The addition and subtraction algorithms are hard to beat on a single processor, but these implementations parallelize poorly. Finally, the modulus depends on many integer-size computations, and could be written more efficiently in terms of component-sized operations like the others.

Download code
hijacker
hijacker
hijacker
hijacker