Sputsoft

Introduction

This article is the first in a series dealing with algorithms for multiple-precision arithmetic. The goal is to present both a theoretical foundation with high-level algorithm descriptions (based on Section 4.3.1, The Classical Algorithms, of The Art of Computer Programming, Volume 2, by Donald E. Knuth) and a portable C++ implementation of the algorithms. The theory and high-level algorithms will be quite universal and generic, whereas the presented code will be just one way to implement the algorithms in a specific programming language.

We start out by considering only non-negative integers. A number u \geq 0 will be represented in radix b \geq 2 using the notation

We will call u an n-digit number and u_0, u_1, etc., its digits. Unless stated otherwise we will always have that the most-significant digit is non-zero, here u_{n-1} \neq 0, and we will represent zero with no digits, 0 = ()_b. We have b^{n-1} \leq u \leq b^n-1, implying n=1+\lfloor \log_b u \rfloor for u \geq 1.

Let the word size of the data type T used for each digit be b_T. For instance, if T is a 32 bit unsigned integer, we have b_T = 2^{32}. We will implement the algorithms using b = b_T and exploit the fact that C++ does arithmetic on unsigned integers modulo b_T (see paragraph 3.9.1 (4) of the C++ standard). This makes it possible to implement portable algorithms. They will not be optimal with respect to speed, however, and it will be noted when specialized operations, such as add-with-carry instructions, would lead to more efficient implementations.

Data Structures

A non-negative number will be represented in C++ as an instance of the class NonNegativeInteger:

template <typename T, typename V=detail::SimpleDigitVector<T> >
class NonNegativeInteger {
private:
  boost::shared_ptr<V> digitvec;
  // ...
public:
  NonNegativeInteger();
  NonNegativeInteger(T value);
  // ...
};

The type argument T is used to represent each digit. It must be integer and unsigned, so unsigned char, unsigned short, unsigned int, unsigned long, and unsigned long long can all be used (the type long long int is not standard C++, but is, e.g., supported by GCC). If a digit type with 8, 16, 32, or 64 bits is needed, the boost integer types uint8_t, uint16_t, uint32_t, uint64_t (from the namespace boost) can be used with portability ensured (uint64_t is not always available, but the macro BOOST_NO_INT64_T will tell you if it is not).

The type argument V is used as the container type for the digits. The default container is SimpleDigitVector (which at this time is also the only container supported). This default container simply wraps an array of size (at least) the number of digits.

Note that the digit container digitvec of NonNegativeInteger is wrapped in a boost shared pointer. Consider the small code excerpt:

  NonNegativeInteger<unsigned> a=value1(), b=value2(), c;
  c = a;
  b += a;
  c += b;

Because of the shared pointer, the first assignment, c = a, is very cheap and will not result in copying the digit container. Instead, the digitvec of both a and c will refer to the same instance of the digit container (which gets a reference count of 2). This also makes it cheap to pass arguments by value and returning numbers by value. The next statement, b += a, can add a to b in place (assuming there is enough space in b‘s digitvec to contain the result) because it is possible to check if b is the only instance referring to its digit container. Contrarywise, the last statement, c += a, cannot add a to c in place because a refers to the same container as c does.

Addition

Adding radix b-numbers is quite straightforward and is easily done using the familiar pencil-and-paper method.

To formalize, we first consider adding two n-digit numbers, n \geq 1, u=(u_{n-1} \ldots u_1 u_0)_b and v = (v_{n-1} \ldots v_1 v_0)_b, obtaining an (n+1)-digit sum w=(w_n \ldots w_1 w_0)_b. Note that we here may end up with w_n = 0 (in which case w_{n-1} \neq 0). We have

Note that \lfloor x/b \rfloor = [x \geq b] for x < 2 b, where [P] is equal to 1 if P is true and equal to 0 if P is false. This means k_i \in \{0,1\} and furthermore u_i + v_i + k_i \leq b-1 + b-1 + 1 = 2 b – 1. Using that 0 \leq u_i, v_i \leq b-1 and (x \;\mbox{mod}\; b) + \lfloor x/b \rfloor = x it is quite easy to show that, in fact, w = u + v.

The cases where u or v is zero makes addition trivial. Similarly, if the number of digits in u and v are different, is it quite easy to adjust the algorithm above.

Let us now look at some implementation details. How do we compute z = (x + y) \;\mbox{mod}\; b and k = \lfloor (x + y)/b \rfloor for 0 \leq x, y < b? We have several options, some of which are

1. Use 2 b \leq b_T. Then z and k can be computed directly.
2. Use b = b_T and the CPU’s add and add-with-carry instructions.
3. Use b = b_T, but the computations must be done in some portable C++ way.

Option 1 is actually not an option because we insist on using b = b_T. Option 2 leads to the most efficient code, regarding both space and speed. The problem is that these special instructions are not directly accessible via the C++ standard. Some compilers, though, make it possible to use inline assembly. For instance, GCC has such capabilities.

Option 3 is the way to go. As mentioned earlier, C++ does calculations modulo b_T, so z \leftarrow (x + y) \;\mbox{mod}\; b comes ‘for free’ as simply z = x + y in C++. Left is how to detect whether a carry occurs during an addition. One way to do that is the following. Consider z = (x + y) \;\mbox{mod}\; b for which there are two possibilities. Either z = x + y (k = 0) which implies z \geq x and z \geq y, or we have z + b = x + y (k = 1) which implies z = x – (b – y) = y – (b – x), leading to z < x and z < y. So k = [z < x] = [z < y]. Another way to detect whether a carry occurs is to split x and y into a low and high part, and then adding the low and high parts seperately—keeping track of a possible intermediate carry, of course.

Subtraction

An algorithm for multiple-precision subtraction is similar to addition, as just considered. Again we set u=(u_{n-1} \ldots u_1 u_0)_b and v = (v_{n-1} \ldots v_1 v_0)_b with n \geq 1. To ensure that the result w = u – v is a non-negative integer we furthermore require that u \geq v. We now have the algorithm:

Note that any digit of the result may end up being zero—we only know that 0 \geq w \geq b^n – b^{n-1} – 1. Note furthermore that k_n = [u < v] = 0. Verification of the algorithm is easily done using the fact that x – y = ((x-y) \;\mbox{mod}\; b) – b [x < y] for 0 \leq x, y < b.

The options when implementing the algorithm is, again, much like for addition. Using the CPU’s subtract-with-borrow instruction would be ideal here, but it cannot be done portably in C++. The borrow when calculating x – y can be computed as simply k \leftarrow [x < y], or, to avoid branching instructions, we can split x and y into two parts and do the subtraction for each part seperately.

Multiplication

We seek an algorithm to compute w = u v where

We first, however, consider the simpler operation z \leftarrow y + \alpha u where

The following algorithm suggests itself:

If \alpha = 0 then obviously z_i = 0. If 1 \leq \alpha \leq b-1 then z_m may be zero, in which case z_ {m-1} \neq 0, since

Note that k_i < b since y_i + \alpha u_i + k_i \leq b-1 + (b-1)(b-1) + b-1 = b^2 – 1. The algorithm above is easily verified, using that z_i + b k_{i+1} = y_i + \alpha u_i + k_i.

We now turn to w = u v and get

From this we see that we can compute w by doing a number of (z \leftarrow y + \alpha u)-type operations, if we start with w=0 and then do in-place updates for each j. The b^j-factor simply determines the ‘digit offset’ on which the updates should be done:

for j = 0, 1, \ldots, n-1. Note that w_{m+n-1} may be zero, in which case w_{m+n-2} \neq 0, since b^{m+n-2} \leq w < b^{m+n}.

Now for some implementation details. We note that the only non-trivial computation is z \leftarrow y + \alpha x + k, where 0 \leq \alpha, k, x, y < b, followed by computing z \;\mbox{mod}\; b and \lfloor z/b \rfloor. Most CPUs have instructions that can multiply two word-size numbers and produce a double-word answer. As was the case for addition and subtraction, we don’t have access to these instructions from standard C++. We have 0 \leq z < b^2 but every multiplication result in our portable C++ implementation must be smaller than b. We can do this by using a new base number h where h^2 = b (we assume that b is chosen appropriately so h is integer) and setting

We can now use the multiplication algorithm above on a ‘smaller scale’ to compute the product \alpha x. We can furthermore expand the methods slightly and incorporate the addition of y and k in an elegant way:

We now have z \;\mbox{mod}\; b = (z_1 z_0)_h and \lfloor z/b \rfloor = (z_3 z_2)_h.

Concluding Remarks

We have now covered addition, subtraction, and multiplication of non-negative integers of arbitrary magnitude. Left is how to do division, which will be the subject of the next article.

I am currently working on a project that implements multiple-precision arithmetic following the principles of this, and the following, articles. You can follow the project as it develops at SourceForge. The project is still in its infancy and no releases have yet been made. You can, however, browse the source or check out the source code (using a Subversion client) as

svn co https://sputsoftnumbers.svn.sourceforge.net/svnroot/sputsoftnumbers/trunk sputsoftnumbers

Bug reports, suggestions for improvements, patches, and other comments are very welcome.

Update 2010-07-03: The code has undergone some changes since this article was written. See the project page for more information and the latest articles.