Cpp bignum arithmetic: part I

Adding two C preprocessor bignums is an easy task, we have to merely reproduce the basic algorithm taught at elementary school; if the bignums to add are bn1 and bn2 and length(bn1) ≥ length(bn2) (otherwise we just swap them), we can fold over reverse(bn1), similarly as how we did for bignum comparison, using the following recursive state definition:

state_i = (carry_i,seq_i,res_i);
state₀ = (0, reverse(bn2),empty sequence),
carry_i+1 = addmsd(reverse(bn1)_i,head(seq_i),carry_i),
seq_i+1 = tail(seq_i),
res_i+1 = append(res_i, sublsd(reverse(bn1)_i,head(seq_i),carry_i)),

where addmsd and addlsd (whose suffixes stand for "most significant digit" and "least significant digit", respectively) are defined as:

addmsd(x,y,carry) = floor((x + y + carry) /10),
addlsd(x,y,carry) = (x + y + carry) mod 10.

The bignum corresponding to bn1 + bn2 is obtained from the last res_i, to which the last carry_i must be appended if it is not null. The whole algorithm is implemented with Boost.Preprocessor as follows:

‎#define BN_ADD(bn1,bn2) \
‎BOOST_PP_IF( \
‎  BOOST_PP_GREATER_EQUAL( \
‎    BOOST_PP_SEQ_SIZE(bn1), \
‎    BOOST_PP_SEQ_SIZE(bn2)), \
‎  BN_ADD_IMPL(bn1,bn2), \
‎  BN_ADD_IMPL(bn2,bn1))
‎
‎#define BN_ADD_IMPL(bn1,bn2) \
‎BN_ADD_POSTPROCESS( \
‎  BOOST_PP_SEQ_FOLD_LEFT( \
‎    BN_ADD_OP, \
‎    (0)(BOOST_PP_SEQ_REVERSE(bn2))((~)), \
‎    BOOST_PP_SEQ_REVERSE(bn1)))
‎
‎#define BN_ADD_POSTPROCESS(state) \
‎BOOST_PP_IF( \
‎  BOOST_PP_NOT_EQUAL(BN_ADD_STATE_CARRY(state),0), \
‎  BOOST_PP_IDENTITY((BN_ADD_STATE_CARRY(state))), \
‎  BOOST_PP_EMPTY)() \
‎BOOST_PP_SEQ_POP_BACK(BN_ADD_STATE_RESULT(state))
‎
‎#define BN_ADD_OP(s,state,x) \
‎BN_ADD_OP_IMPL( \
‎  x, \
‎  BN_ADD_STATE_Y(state), \
‎  BN_ADD_STATE_CARRY(state), \
‎  BN_ADD_STATE_RBN2(state), \
‎  BN_ADD_STATE_RESULT(state))
‎
‎#define BN_ADD_OP_IMPL(x,y,carry,rbn2,result) \
‎(BN_ADD_MSD(x,y,carry)) \
‎(BOOST_PP_SEQ_TAIL(rbn2(0))) \
‎((BN_ADD_LSD(x,y,carry))result)
‎
‎#define BN_ADD_STATE_CARRY(state) BOOST_PP_SEQ_ELEM(0,state)
‎#define BN_ADD_STATE_RBN2(state) BOOST_PP_SEQ_ELEM(1,state)
‎#define BN_ADD_STATE_Y(state) \
‎BOOST_PP_SEQ_HEAD(BN_ADD_STATE_RBN2(state))
‎#define BN_ADD_STATE_RESULT(state) BOOST_PP_SEQ_ELEM(2,state)
‎
‎#define BN_ADD_MSD(x,y,carry) \
‎BOOST_PP_DIV(BOOST_PP_ADD(BOOST_PP_ADD(x,y),carry),10)
‎
‎#define BN_ADD_LSD(x,y,carry) \
‎BOOST_PP_MOD(BOOST_PP_ADD(BOOST_PP_ADD(x,y),carry),10)

Subtraction is similar. As our bignums only represent positive integers, we define BN_SUB(bn1,bn2) to be 0 when bn2 > bn1, just as BOOST_PP_SUB does. Again, we assume that length(bn1) ≥ length(bn2) and fold over reverse(bn1); the associated recursive state definition is this:

state_i = (borrow_i,prefix_i,seq_i,res_i);
state₀ = (0,0, reverse(bn2),empty sequence),
borrow_i+1 = submsd(reverse(bn1)_i,head(seq_i),borrow_i),
prefix_i+1 =
‎ · if sublsd(reverse(bn1)_i,head(seq_i),borrow_i) ≠ 0, 1,
‎‎ · else prefix_i + 1.
seq_i+1 = tail(seq_i),
res_i+1 = append(res_i, sublsd(reverse(bn1)_i,head(seq_i),borrow_i)),
submsd(x,y,borrow) = floor((10 + x − y − borrow) / 10) − 1,
sublsd(x,y,borrow) = (10 + x − y − borrow) mod 10.

prefix_i keeps track of the leading zeros in the res_i sequence; these zeros must be deleted on postprocessing time, since the normal representation for bignums always begins with a non-zero digit (except for the bignum corresponding to 0). The observant reader might have noticed that prefix_i is actually the number of leading zeros plus one; this is merely an implementation artifact that helps us distinguish the special case when the result of the subtraction is 0. The algorithm allows us to easily detect when bn2 > bn1: this is the case if and only if the final borrow_i is not zero. We will not show the whole coding of BN_SUB as a Boost.Preprocessor-based macro since it basically follows the same implementation patterns as BN_ADD.

Multiplication is a little harder. Although it is by no means the fastest algorithm possible, we will implement the O(n²) elementary school method, which consists of multiplying the first factor by every digit of the second factor, shifting appropriately and adding the intermediate results. So, we need first to devise a macro BN_DIGIT_MUL(bn,y) to multiply a bignum by a digit. No surprise by now, we can implement such operation by folding over reverse(bn) with this state structure:

state_i = (y,carry_i,res_i);
state₀ = (y,0,empty sequence),
carry_i+1 = mulmsd(reverse(bn1)_i,y,carry_i),
res_i+1 = append(res_i, mullsd(reverse(bn1)_i,y,carry_i)),
mulmsd(x,y,carry) = floor((xy + carry) /10),
mullsd(x,y,carry) = (xy + carry) mod 10.

Using BN_DIGIT_MUL, we can implement BN_MUL(bn1,bn2) by folding over reverse(bn1) as follows:

state_i = (shift_i,bn2,res_i);
state₀ = (empty sequence,bn2,0),
shift_i+1 = append(shift_i,0),
res_i+1 = resi + append(BN_DIGIT_MUL(bn2,reverse(bn1)_i), shift_i).

shift_i is a sequence of increasingly many zero digits used to perform fast multiplication of the intermediate results by powers of 10 --multiplying a bignum bn by 10ⁿ reduces to appending a sequence of n zero digits to the representation of bn.

To keep this entry short, we will defer to later occasions the study of bignum division, when we will also provide complete code for all the arithmetic operations and make some performance measurements.