Natural numbers as an abstract data type

You'd certainly want addition and subtraction, multiplication and division. In the case of division, you'd want to be able to determine both the (integer) quotient and the remainder. You'd probably also want exponentiation -- since the bases and exponents are natural numbers, you wouldn't have to look out for any problematical cases (except for zero to the zeroth power). Various special cases of these operations occur commonly enough to justify functions of their own: adding 1, subtracting 1, doubling, halving (discarding the remainder), squaring, and cubing.

The preceding operations are ``non-destructive'' -- the operands are intact after the operation. You might want destructive versions of a few of them: incrementing or decrementing a natural-number variable by a specified amount and scaling a natural number variable up or down by a specified factor (that is, multiplying it or dividing it by that factor). Again, there are special cases that merit separate procedures: incrementing and decrementing by one, doubling and halving.

You'd want to be able to compare natural numbers; probably, you'd want the full set of six comparison functions (<, =, >, >=, <>, <=), and perhaps also the special case of testing equality with 0. You'd want the max and min functions that select and return the larger or the smaller of two natural-number arguments. You might want a function that takes two natural numbers and determines whether one is a multiple of the other; the special case of determining whether a number is a multiple of 2 would also be handy.

You'd want to be able to read natural numbers from text files (and specifically from standard input) and to write them to text files (and specifically to standard output and standard error); also, independently, to convert natural numbers to and from strings of characters that express those numbers in decimal numeration.

Finally, there's one less obvious operation: You'd want to be able to assign any natural number to any natural-number variable. The assignment operation available in your implementation language might not be appropriate. (For instance, the implementer might choose to represent natural numbers as structures accessed through pointers, in which case the built-in assignment operation would copy the pointer rather than the structure to which it points.)

Language designers have occasionally provided more than this (the Scheme programming language provides functions for computing greatest common divisors and least common multiples, for instance), but this list of operations constitutes the basic repertoire. Here's how the procedure and function headers should look in HP Pascal:

function sum (augend, addend: natural): natural;
function difference (minuend, subtrahend: natural): natural;
function product (multiplicand, multiplier: natural): natural;
function quotient (dividend, divisor: natural): natural;
function remainder (dividend, divisor: natural): natural;
function power (base, exponent: natural): natural;

function successor (given: natural): natural;
function predecessor (given: natural): natural;
function twice (given: natural): natural;
function half (given: natural): natural;
function square (given: natural): natural;
function cube (given: natural): natural;

procedure increment (var given: natural; amount: natural);
procedure decrement (var given: natural; amount: natural);
procedure scale_up (var given: natural; scale: natural);
procedure scale_down (var given: natural; scale: natural);
procedure up1 (var given: natural);
procedure down1 (var given: natural);
procedure double (var given: natural);
procedure halve (var given: natural);

function less (first, second: natural): Boolean;
function equal (first, second: natural): Boolean;
function greater (first, second: natural): Boolean;
function not_less (first, second: natural): Boolean;
function not_equal (first, second: natural): Boolean;
function not_greater (first, second: natural): Boolean;
function is_zero (given: natural): Boolean;

function max (first, second: natural): natural;
function min (first, second: natural): natural;

function is_multiple (first, second: natural): Boolean;
function is_even (given: natural): Boolean;
function is_odd (given: natural): Boolean;

procedure read_natural (var source: text; var legendum: natural);
procedure input_natural (var legendum: natural);
procedure write_natural (var target: text; scribendum: natural);
procedure output_natural (scribendum: natural);
procedure output_natural_to_standard_error (scribendum: natural);

function string_to_natural (var s: string): natural;
procedure natural_to_string (n: natural; var s: string);

procedure assign (var target: natural; source: natural);

Designing and choosing the implementation

To come as close as possible to the ideal of being able to work with natural numbers of any size, we need a representation that uses a data structure with dynamically allocated components, so that we can take as many as we need to represent whatever natural numbers turn up in our computations. When we consider how to break an arbitrary natural number into components of fixed size, it's plausible to think of a sequence of digits. A strictly linear structure is okay, and since we'll hardly ever want to inspect the middle digits of a number without first looking at either the most or the least significant digits, sequential access to the components is a limitation that we can live with. However, it should be possible to run through the digits from least significant to most significant (for such operations as addition, which is easiest if one works from right to left) or from most significant to least significant (for such algorithms as halving, which is easiest if one works left to right). So a doubly-linked list of digits seems to be the best choice.

The digits need not be decimal digits. In fact, using base ten in the internal representation would mean that a large natural number would have a lot of components, each containing less than four bits of information. Considering that each component will also require two pointers (which on our system occupy thirty-two bits each), this seems a little extravagant. Using a much larger base of numeration would allow us to pack more information into each component.

On the other hand, if we make the base of numeration too large, we run up against the problem of being unable to perform arithmetic on the individual digits without causing overflows. Specifically, the algorithm for long division that will be presented below presupposes that we can directly compute numbers having as many as three digits in whatever base we settle on -- in other words, that the base should not exceed the cube root of MAXINT + 1. Furthermore, it will simplify the parity functions (even and odd), as well as the division function, if we ensure that the base of numeration is an even number.

There are actually two good choices for the base of numeration. To minimize the number of components, we could use the largest even integer not exceeding the square root of MAXINT, which on the HPs turns out to be 1290. To simplify the input and output routines, we could use 1000. In the naturals module that I've developed below, I chose the larger base, but since I've used a defined constant for the base it would be straightforward to switch.

In addition to the components making up the doubly-linked list, it is convenient to have a ``header node'' containing pointers to the most and least significant digits and a count of the number of digits in the number. The individual digits are accessible only through these header nodes; one of the module invariants is that no pointer to an individual component ever escapes from the module. The header nodes are allocated dynamically. The natural data type is simply equated with a pointer to a header node.

Here, then, are the type definitions for the implementation I propose:

const
  BASE = 1290;
     { = largest even integer less than the cube root of MAXINT + 1 }
type
  base_digits = 0 .. BASE - 1;
  link = ^component;
  component = record
                digit: base_digits;
                higher, lower: link
              end;
  header = record
             size: integer;
             highest, lowest: link
           end;
  natural = ^header;

function integer_to_natural (n: integer): natural;
function natural_to_integer (n: natural): integer;

To keep the doubly-linked lists as short as possible and to avoid some useless computation, we can establish as an invariant of the module that our representations of natural numbers will never contain leading zeroes. Any of the procedures or functions that generate new values of type natural will trim off any leading zeroes before returning. Moreover, we shall represent the natural number 0 with a header node in which both of the links are NIL (and the doubly-linked list therefore has no components).

Implementing the functions and procedures

Addition

Here's how this algorithm is expressed in HP Pascal:

type
  bit: 0 .. 1;
function sum (augend, addend: natural): natural;
var
  result: natural;
    { the sum of the arguments, as it is constructed }
  augend_traverser: link;
    { a pointer to successive components of the augend }
  addend_traverser: link;
    { a pointer to successive components of the addend }
  carry: bit;
    { a carry from one digit position to the next }
  augend_digit, addend_digit: base_digits;
    { digits from corresponding positions in the augend and addend }
  place_sum: integer;
    { the sum of the digits in the same position in the augend and addend,
      plus any carry from the next lower digit position }
begin
  allocate_zero (result);
  augend_traverser := augend^.lowest;
  addend_traverser := addend^.lowest;
  carry := 0;
  while (augend_traverser <> NIL) or (addend_traverser <> NIL) or
        (carry <> 0) do begin
    if augend_traverser = NIL then
      augend_digit := 0
    else begin
      augend_digit := augend_traverser^.digit;
      augend_traverser := augend_traverser^.higher
    end;

    if addend_traverser = NIL then
      addend_digit := 0
    else begin
      addend_digit := addend_traverser^.digit;
      addend_traverser := addend_traverser^.higher
    end;
    place_sum := augend_digit + addend_digit + carry;
    if place_sum < BASE then begin
      prepend (place_sum, result);
      carry := 0
    end
    else begin
      prepend (place_sum - BASE, result);
      carry := 1
    end
  end;
  sum := result
end;

Subtraction

Since the result of the subtraction may contain fewer significant digits than the minuend, whereas the algorithm generates exactly the same number of digits in the difference as in the minuend, the utility procedure strip_leading_zeroes is invoked to remove any zero digits from the result before it is returned.

The version presented here presupposes that the subtrahend is less than or equal to the minuend. In the actual naturals module, the corresponding procedure checks this precondition and signals an error if it is not met.

function difference (minuend, subtrahend: natural): natural;
var
  result: natural;
    { the difference of the arguments, as it is constructed }
  minuend_traverser: link;
    { a pointer to successive components of the minuend }
  subtrahend_traverser: link;
    { a pointer to successive components of the subtrahend }
  borrow: bit;
    { a borrow from one digit position to the next; a 0 indicates no
      borrow from the current position, a 1 indicates that there was
      a borrow }
  minuend_digit, subtrahend_digit: base_digits;
    { digits from corresponding positions in the minuend and subtrahend }
begin
  allocate_zero (result);
  minuend_traverser := minuend^.lowest;
  subtrahend_traverser := subtrahend^.lowest;
  borrow := 0;
  while minuend_traverser <> NIL do begin
    minuend_digit := minuend_traverser^.digit;
    minuend_traverser := minuend_traverser^.higher;
    if subtrahend_traverser = NIL then
      subtrahend_digit := 0
    else begin
      subtrahend_digit := subtrahend_traverser^.digit;
      subtrahend_traverser := subtrahend_traverser^.higher
    end;
    if minuend_digit - borrow < subtrahend_digit then begin
      prepend (BASE + minuend_digit - borrow - subtrahend_digit, result);
      borrow := 1
    end
    else begin
      prepend (minuend_digit - borrow - subtrahend_digit, result);
      borrow := 0
    end
  end;
  strip_leading_zeroes (result);
  difference := result
end;

Multiplication

function product (multiplicand, multiplier: natural): natural;
var
  result: natural;
    { the product, as it is built up from partial products }
  lower_positions: integer;
    { the number of digit positions below the current one }
  traverser: link;
    { a pointer to each successive component of the multiplier, from
      least significant to most significant }
  partial: natural;
    { a partial product: the product of the multiplicand and one digit
      of the multiplier }
  zero_tallier: integer;
    { counts off trailing zeroes as they are appended to a partial
      product }
begin
  allocate_zero (result);
  if not is_zero (multiplicand) then begin
    lower_positions := 0;
    traverser := multiplier^.lowest;
    while traverser <> NIL do begin
      if traverser^.digit <> 0 then begin
        partial :=
             multiply_natural_by_digit (multiplicand, traverser^.digit);
        for zero_tallier := 1 to lower_positions do
          append (0, partial);
        increment (result, partial);
        deallocate_natural (partial)
      end;
      lower_positions := lower_positions + 1;
      traverser := traverser^.higher
    end
  end;
  product := result
end;

This function invokes the utility function multiply_natural_by_digit, which is not exported from the module, to compute each partial product. In this function, the result of the multiplication is accumulated in a natural number that is initialized to zero and then built up a digit at a time, starting with the least significant digit. To obtain each new digit, a digit is recovered from the multiplicand (starting at the low end) and multiplied by the entire multiplier; to this product is added a carry from the next lower digit position. The resulting ``place product'' is then divided by the base of numeration; the remainder is the digit to be prepended to the overall product, and the quotient becomes the carry into the next higher digit position. If there is a non-zero carry out of the highest digit position, it becomes the most significant digit of the result.

This algorithm exactly reflects the mechanics of ``short multiplication'' by hand. Here is the Pascal code:

function multiply_natural_by_digit (multiplicand: natural;
                                    multiplier: base_digits): natural;
var
  result: natural;
    { the product, as it is built up one digit at a time }
  traverser: link;
    { a pointer to each successive component of the multiplicand, from
      least significant to most significant }
  carry: base_digits;
    { the number of multiples of BASE to be carried into the next higher
      digit position }
  place_product: integer;
    { the product of one digit of the multiplicand and the multiplier,
      plus any carry from the next lower digit position }
begin
  allocate_zero (result);
  if multiplier <> 0 then begin
    carry := 0;
    traverser := multiplicand^.lowest;
    while traverser <> NIL do begin
      place_product := traverser^.digit * multiplier + carry;
      prepend (place_product mod BASE, result);
      carry := place_product div BASE;
      traverser := traverser^.higher
    end;
    if carry <> 0 then
      prepend (carry, result)
  end;
  multiply_natural_by_digit := result
end;

Furthermore, product invokes the increment procedure and the is_zero function, which are exported from the naturals module and will be discussed in their proper places below.

Division

This algorithm finds the quotient and the remainder simultaneously. There are a fair number of cases in which the user of our natural-number library module will need both, so we might as well export the procedure that computes them together:

procedure divide (dividend, divisor: natural; var quot, rem: natural);

function quotient (dividend, divisor: natural): natural;
var
  quot: natural;
    { the quotient resulting from the division }
  wastebasket: natural;
    { the remainder -- to be discarded }
begin
  divide (dividend, divisor, quot, wastebasket);
  deallocate_natural (wastebasket);
  quotient := quot
end;

function remainder (dividend, divisor: natural): natural;
var
  rem: natural;
    { the remainder resulting from the division }
  wastebasket: natural;
    { the quotient -- to be discarded }
begin
  divide (dividend, divisor, wastebasket, rem);
  deallocate_natural (wastebasket);
  remainder := rem
end;

procedure divide_natural_by_digit (dividend: natural;
                                   divisor: base_digits;
                                   var quot: natural;
                                   var rem: base_digits);
var
  traverser: link;
    { points successively to the digits of the dividend, from most
      significant to least significant }
  short_dividend: integer;
    { a partial dividend, composed of the remainder left over from the
      division in the previous position, scaled up by a factor of BASE,
      plus the digit in the current position }
begin
  assert (divisor <> 0, DIVISION_EXCEPTION, natural_handler);
  allocate_zero (quot);
  traverser := dividend^.highest;
  if traverser = NIL then
    rem := 0
  else if traverser^.digit < divisor then begin
    rem := traverser^.digit;
    traverser := traverser^.lower
  end
  else
    rem := 0;
  while traverser <> NIL do begin
    short_dividend := rem * BASE + traverser^.digit;
    append (short_dividend div divisor, quot);
    rem := short_dividend mod divisor;
    traverser := traverser^.lower
  end
end;

On each iteration of the division loop, the procedure constructs a ``short dividend,'' performs the short division, attaches the quotient (which will always be a single digit) to the low end of quot, and changes rem to reflect the remainder of the division. When the loop terminates, all the digits of the quotient will have been obtained, and the final remainder is the remainder for the entire division.

Here, then, is the general structure for the divide procedure:

procedure divide (dividend, divisor: natural; var quot, rem: natural);
var
  irem: base_digits;
    { the remainder of a short division -- used only if divisor is a
      single digit }
begin
  assert (not is_zero (divisor), DIVISION_EXCEPTION, natural_handler);
  if less (dividend, divisor) then begin
    allocate_zero (quot);
    assign (rem, dividend)
  end
  else if divisor^.size <= 1 then begin
    divide_natural_by_digit (dividend, divisor^.highest^.digit,
                             quot, irem);
    rem := integer_to_natural (irem)
  end
  else
    { Do the really difficult case in which the divisor has at least
      two digits and the dividend is not less than the divisor. }
end;

In pencil-and-paper computation, we begin the ``really difficult case'' by guessing the first digit of the quotient, multiplying that digit by the divisor, and subtracting the result from the first few digits of the dividend. If the result of the subtraction is either negative or greater than the divisor, we know that we've guessed wrong, and we back up and make another guess.

Computers cannot guess, so we need to automate the process. The brute-force method would be to initialize a variable trial_digit to BASE - 1 and then decrement it as long as the result of multiplying it by the divisor is greater than the first few digits of the dividend. This works quite well if BASE is 2, but it's terribly inefficient when BASE is large, as it is here.

A better approach is to construct another, much simpler division problem in which the first digit of the quotient is the same as in ours. For instance, if we take the first three digits of the dividend and treat them as if they were a separate number first_three (which will be less than or equal to MAXINT -- we chose our base of numeration precisely to ensure this), and similarly split off the first two digits of the divisor to form the number first_two, then dividing first_three by first_two should yield a quotient in which the first digit is a close approximation of the first digit of the original division problem. In fact, Brinch Hansen shows by elementary algebra that the trial digit obtained in this way is either correct or, at worst, one unit too large.

The cases in which the trial digit is even one unit too large are expensive, however, since correcting the approximation involves discarding all the work that was done in constructing the product of the multiplicand and the trial digit. It turns out that the frequency of an erroneous trial digit decreases as the value of first_two increases, since the relative error involved in discarding the subsequent digits becomes less. A convenient way to make first_two large is to begin the entire division process by normalizing the dividend and the divisor -- scaling them up by the same (one-digit integer) factor. This obviously will not affect the quotient at all; and, although it will have the effect of scaling up the remainder by the same factor, we can easily scale it back down again at the end of the division process.

The scaling factor is computed as BASE div (first + 1), where first is the leading digit of the original divisor. Before scaling, a leading zero is prepended to the dividend; this ensures that the quotient obtained from the first three digits of the scaled dividend and the first two of the scaled divisor will always be a single digit. (Actually, it is still possible for the quotient to be BASE, in the rare case that the correct digit is BASE - 1 and the quotient is too high in spite of the normalization. The function that we will use to determine the trial digit will test for this case and correct it by subtracting one before returning the quotient.)

Unfortunately, it is possible for the trial digit computed from the first three digits of the scaled dividend to be 0, so we may still have to remove a leading zero from the quotient at the end of the division process.

The procedure for scaling up a natural number is similar to the algorithm used in the multiply_natural_by_digit procedure, except that it destroys or overwrites the original natural number:

procedure scale_up_by_digit (var given: natural; scale: base_digits);
var
  traverser: link;
    { a pointer to each successive component of given, from least
      significant to most significant }
  carry: base_digits;
    { the number of multiples of BASE to be carried into the next
      higher digit position }
  place_product: integer;
    { the product of one digit of given and the scale, plus any carry
      from the next lower digit position }
begin
  if scale = 0 then begin
    deallocate_natural (given);
    allocate_zero (given)
  end
  else begin
    carry := 0;
    traverser := given^.lowest;
    while traverser <> NIL do begin
      place_product := traverser^.digit * scale + carry;
      traverser^.digit := place_product mod BASE;
      carry := place_product div BASE;
      traverser := traverser^.higher
    end;
    if carry <> 0 then
      prepend (carry, given)
  end
end;

function find_trial_digit (residue: natural): base_digits;
var
  first_three: integer;
    { the value of the first three digits of residue; it will always have
      at least three, because the scaled divisor has at least two and
      residue always has one more than the scaled divisor }
  traverser: link;
    { a pointer to successive components of residue, starting with the
      most significant digit }
  position: integer;
    { counts off digits in residue }
  trial: integer;
    { an approximate value for the next quotient digit }
begin
  first_three := 0;
  traverser := residue^.highest;
  for position := 1 to 3 do begin
    first_three := first_three * BASE + traverser^.digit;
      traverser := traverser^.lower
  end;
  trial := first_three div first_two;
  if trial = BASE then
    find_trial_digit := BASE - 1
  else
    find_trial_digit := trial
end;

Here's how it looks in Pascal:

begin

  { Compute the appropriate scaling factor. }

  scaling_factor := BASE div (divisor^.highest^.digit + 1);

  { Attach a leading zero to the dividend and scale it up; scale up the
    divisor. }

  assign (scaled_dividend, dividend);
  prepend (0, scaled_dividend);
  scale_up_by_digit (scaled_dividend, scaling_factor);
  assign (scaled_divisor, divisor);
  scale_up_by_digit (scaled_divisor, scaling_factor);

  { Compute the value of the first two digits of the scaled divisor,
    to be used in computing trial digits for the quotient. }

  first_two := BASE * scaled_divisor^.highest^.digit +
               scaled_divisor^.highest^.lower^.digit;

  { Initialize residue by copying into it one more digit than the
    scaled divisor has. }

  allocate_zero (residue);
  traverser := scaled_dividend^.highest;
  for position := 1 to scaled_divisor^.size + 1 do begin
    append (traverser^.digit, residue);
    traverser := traverser^.lower
  end;

  { Construct the quotient one digit at a time, from most significant
    to least significant. }

  allocate_zero (quot);
  finished := FALSE;
  while not finished do begin

    { Call find_trial_digit to obtain a good estimate of the next
      digit of the quotient. }

    trial_digit := find_trial_digit (residue);

    { Multiply that trial digit by the scaled divisor; if the result
      exceeds the value of residue, decrease the trial digit by 1 at
      multiply again. }

    trial_product :=
        multiply_natural_by_digit (scaled_divisor, trial_digit);
    strip_leading_zeroes (residue) { so that the less function will work }
    if less (residue, trial_product) then begin
      trial_digit := trial_digit - 1;
      deallocate_natural (trial_product);
      trial_product :=
          multiply_natural_by_digit (scaled_divisor, trial_digit);
    end;

    { Attach the trial digit at the low end of the quotient and
      subtract the product of the trial digit and the scaled divisor
      from the value of residue.  Ensure that residue has the correct
      number of digits. }

    append (trial_digit, quot);
    decrement (residue, trial_product);
    deallocate_natural (trial_product);
    while residue^.size < scaled_divisor^.size do
      prepend (0, residue);

    { If there are any remaining digits in the scaled dividend,
      attach the first of them at the low end of the residue and repeat
      the loop. }

    if traverser = NIL then
      finished := TRUE
    else begin
      append (traverser^.digit, residue);
      traverser := traverser^.lower
    end

  end;

  { Strip any leading zeroes from the quotient and divide the value of
    residue by the scaling factor to obtain the correct remainder for
    the original division. } 

  strip_leading_zeroes (quot);
  deallocate_natural (scaled_dividend);
  deallocate_natural (scaled_divisor);
  strip_leading_zeroes (residue);
  divide_natural_by_digit (residue, scaling_factor, rem, wastebasket);
  deallocate_natural (residue)

end

Exponentiation

function power (base, exponent: natural): natural;
begin
  if is_zero (exponent) then
    power := integer_to_natural (1)
  else if is_even (exponent) then
    power := power (square (base), half (exponent))
  else
    power := product (base, power (square (base), half (exponent)))
end;

This reflection suggests the following alternative structure, in which the recursive function takes a third argument, accumulator, to hold the product of all of the factors of the final result that have so far been encountered:

function power (base, exponent: natural): natural;
var
  copy_of_base, copy_of_exponent: natural;
    { duplicates of the base and exponent that can be changed destructively
      during the execution of this function }

  function power_helper (base, exponent, accumulator: natural): natural;
  begin
    if is_zero (exponent) then
      power_helper := accumulator
    else begin
      if is_odd (exponent) then
        scale_up (accumulator, base);
      scale_up (base, base);
      halve (exponent);
      power_helper := power_helper (base, exponent, accumulator)
    end
  end;

begin { function power }
  assign (copy_of_base, base);
  assign (copy_of_exponent, exponent);
  power := power_helper (copy_of_base, copy_of_exponent,
                         integer_to_natural (1));
  deallocate_natural (copy_of_base);
  deallocate_natural (copy_of_exponent)
end;

function power (base, exponent: natural): natural;
var
  copy_of_base, copy_of_exponent: natural;
    { duplicates of the base and exponent that can be changed destructively
      during the execution of this function }
  accumulator: natural;
    { the product of various powers of the base }
begin
  assign (copy_of_base, base);
  assign (copy_of_exponent, exponent);
  accumulator := integer_to_natural (1);
  while not is_zero (copy_of_exponent) do begin
    if is_odd (copy_of_exponent) then
      scale_up (accumulator, copy_of_base);
    scale_up (copy_of_base, copy_of_base);
    halve (copy_of_exponent)
  end;
  deallocate_natural (copy_of_base);
  deallocate_natural (copy_of_exponent);
  power := accumulator
end;

The next few functions are special cases of these basic arithmetic operations, and could be coded as calls to the corresponding functions with one operand fixed. However, since the fixed operands are small, we can write much faster versions of the specialized functions.

Successor and predecessor

function successor (given: natural): natural;
var
  carry: bit;
    { initially, the 1 that must be added; subsequently, a 1 so long as
      a carry must be propagated, and 0 thereafter }
  result: natural;
    { the successor natural number as it is constructed, digit by digit }
  traverser: link;
    { a pointer to successive components of the given natural number,
      from least significant to most significant }
begin
  carry := 1;
  allocate_zero (result);
  traverser := given^.lowest;
  while traverser <> NIL do begin
    if carry = 0 then
      prepend (traverser^.digit, result)
    else if traverser^.digit = BASE - 1 then
      prepend (0, result)
    else begin
      prepend (traverser^.digit + 1, result);
      carry := 0
    end;
    traverser := traverser^.higher
  end;
  if carry = 1 then
    prepend (1, result);
  successor := result
end;

The twice and half functions

function twice (given: natural): natural;
var
  result: natural;
    { the double of the given natural number, as it is constructed a
      digit at a time }
  carry: bit;
    { indicates whether a carry was generated by the preceding
      digit-doubling }
  traverser: link;
    { a pointer to successive components of the given natural number,
      from least significant to most significant }
  place_product: integer;
    { the double of one digit of the given natural number, plus any carry
      from the next lower digit position } 
begin
  allocate_zero (result);
  carry := 0;
  traverser := given^.lowest;
  while traverser <> NIL do begin
    place_product := traverser^.digit * 2 + carry;
    if place_product < BASE then begin
      prepend (place_product, result);
      carry := 0
    end
    else begin
      prepend (place_product - BASE, result);
      carry := 1
    end;
    traverser := traverser^.higher
  end;
  if carry = 1 then
    prepend (1, result);
  twice := result
end;

The square and cube functions

function square (given: natural): natural;
begin
  square := product (given, given)
end;

function cube (given: natural): natural;
var
  sq: natural;
    { the square of the given number }
begin
  sq := product (given, given);
  cube := product (sq, given);
  deallocate_natural (sq)
end;

Destructive arithmetic

In coding the increment and decrement procedures, one should give careful consideration to what will happen when the procedure is invoked with the same natural number in both argument positions.

To illustrate the actual code, here are the increment and scale_up procedures:

procedure increment (var given: natural; amount: natural);
var
  given_traverser: link;
    { a pointer to successive components of the natural number to be
      incremented }
  amount_traverser: link;
    { a pointer to successive components of the increment }
  carry: bit;
    { a carry from one digit position to the next }
  given_digit, amount_digit: base_digits;
    { digits from corresponding positions in given and amount }
  place_sum: integer;
    { the sum of the digits in the same position in given and amount,
      plus any carry from the next lower digit position }
begin
  given_traverser := given^.lowest;
  amount_traverser := amount^.lowest;
  carry := 0;
  while (given_traverser <> NIL) or (amount_traverser <> NIL) or
        (carry <> 0) do begin
    if given_traverser = NIL then
      given_digit := 0
    else
      given_digit := given_traverser^.digit;
    if amount_traverser = NIL then
      amount_digit := 0
    else begin
      amount_digit := amount_traverser^.digit;
      amount_traverser := amount_traverser^.higher
    end;
    place_sum := given_digit + amount_digit + carry;
    if place_sum < BASE then
      carry := 0
    else begin
      place_sum := place_sum - BASE;
      carry := 1
    end;
    if given_traverser = NIL then
      prepend (place_sum, given)
    else begin
      given_traverser^.digit := place_sum;
      given_traverser := given_traverser^.higher
    end
  end
end;

procedure scale_up (var given: natural; scale: natural);
var
  delendum: natural;
    { a copy of the pointer to the header node of the old value of given }
begin
  delendum := given;
  given := product (given, scale);
  deallocate_natural (delendum)
end;

Comparison functions

For example, in the less function, the first operand must be less than the second if it contains fewer digits, and cannot be less than the second if it contains more. In the remaining case, the function can simply start at the high end of both numbers and compare corresponding digits until either they fail to match or both lists are exhausted.

Here's what it looks like in Pascal:

function less (first, second: natural): Boolean;
var
  first_traverser, second_traverser: link;
    { pointers to corresponding components of first and second, starting
      with the most significant digits }
  finished: Boolean;
    { indicates whether the value to be returned has been conclusively
      established yet }
begin
  if first^.size < second^.size then
    less := TRUE
  else if second^.size < first^.size then
    less := FALSE
  else begin
    first_traverser := first^.highest;
    second_traverser := second^.highest;
    finished := FALSE;
    repeat
      if first_traverser = NIL then begin
        less := FALSE;
        finished := TRUE
      end
      else if first_traverser^.digit < second_traverser^.digit then begin
        less := TRUE;
        finished := TRUE
      end
      else if second_traverser^.digit < first_traverser^.digit then begin
        less := FALSE;
        finished := TRUE
      end
      else begin
        first_traverser := first_traverser^.lower;
        second_traverser := second_traverser^.lower;
      end
    until finished
  end
end;

function equal (first, second: natural): Boolean;
var
  first_traverser, second_traverser: link;
    { pointers to corresponding components of first and second, starting
      with the most significant digits }
  finished: Boolean;
    { indicates whether the value to be returned has been conclusively
      established yet }
begin
  if first^.size <> second^.size then
    equal := FALSE
  else begin
    first_traverser := first^.highest;
    second_traverser := second^.highest;
    finished := FALSE;
    repeat
      if first_traverser = NIL then begin
        equal := TRUE;
        finished := TRUE
      end
      else if first_traverser^.digit <> second_traverser^.digit then begin
        equal := FALSE;
        finished := TRUE
      end
      else begin
        first_traverser := first_traverser^.lower;
        second_traverser := second_traverser^.lower;
      end
    until finished
  end
end;

The is_zero function

function is_zero (given: natural): Boolean;
begin
  is_zero := (given^.size = 0)
end;

The max and min functions

function max (first, second: natural): natural;
begin
  if not_less (first, second) then
    max := first
  else
    max := second
end;

The is_multiple function

function is_multiple (first, second: natural): Boolean;
var
  rem: natural;
    { the remainder when first is divided by second }
begin
  if is_zero (second) then
    is_multiple := is_zero (first)
  else begin
    rem := remainder (first, second);
    is_multiple := is_zero (rem);
    deallocate_natural (rem)
  end
end;

is_multiple := is_zero (remainder (first, second))

Parity functions

Here is the Pascal code for is_even; the is_odd function is similar. Note that, since zero has no last digit (in our implementation), it must be handled as a separate case.

function is_even (given: natural): Boolean;
begin
  if given^.size = 0 then
    is_even := TRUE
  else
    is_even := not odd (given^.lowest^.digit)
end;

Input and output procedures

A variable that will hold the natural number that is being read in is initialized to zero at the beginning of the process. Then, as each digit is read in, the value of this variable is scaled up by a factor of ten and incremented by the decimal value of that digit.

Pascal's numeric input routines tacitly skip over any whitespace characters (tabs, newlines, vertical tabs, form feeds, carriage returns, and spaces) that precede the first digit. For consistency, and as an aid to the programmer, the read_natural procedure also consumes such characters, by invoking the skip_white_space procedure:

procedure skip_white_space (var source: text);
const
  TAB = #I;
  NEWLINE = #J;
  VERTICAL_TAB = #K;
  FORMFEED = #L;
  RETURN = #M;
  SPACE = #32;
var
  finished: Boolean;
    { indicates whether it is necessary to keep looking for white space
      to skip }
begin
  finished := FALSE;
  while not finished do
    if eof (source) then
      finished := TRUE
    else if source^ in [TAB, NEWLINE, VERTICAL_TAB, FORMFEED, RETURN,
                        SPACE] then
      get (source)
    else
      finished := TRUE
end;

The read_natural procedure looks more or less like this:

procedure read_natural (var source: text; var legendum: natural);
var
  finished: Boolean;
    { indicates whether all of the digits of the numeral have been read
      in }
begin
  skip_white_space (source);
  allocate_zero (legendum);
  finished := FALSE;
  repeat
    scale_up_by_digit (legendum, 10);
    increment_by_digit (legendum, ord (source^) - ORD_ZERO);
    get (source);
    if eof (source) then
      finished := TRUE
    else
      finished := ((source^ < '0') or ('9' < source^))
  until finished
end;

The input_natural procedure is implemented by a call to read_natural, with the text file argument specialized to input.

The main difficulty with the write_natural procedure is that it is easiest to compute the decimal digits of the number to be written out in reverse order -- least significant digit first -- while the first digit that is written out must be the most significant digit. To cope with this problem, the naturals module imports another module, stacks.

A stack is a data structure that can receive and store any number of values of a single data type, in this case integer, and can recover and return those values one by one, subject to the constraint that the only value that can be recovered from the stack at any given time is the one that was most recently added to the stack. The write_natural procedure uses four routines from the stacks module: initialize_stack, which assigns an empty stack to its argument; push_to_stack, which stores a given integer value in a given stack; pop_stack, which recovers an integer value from a given stack and returns it; and is_empty_stack, which determines whether a given stack is empty.

The plan of the write_natural procedure is to create an empty stack and then to repeatedly divide the given natural number by ten, pushing each remainder on the stack as it is recovered. The successive remainders are the digits of the decimal numeral for the given natural number, from least significant to most significant. Once the natural number has been reduced to zero by these successive divisions, the output phase begins: The stack is popped repeatedly until it is empty, and each recovered remainder is written out as a one-digit number.

As usual, zero must be handled as a special case, since the method just described would print nothing at all.

Here's the Pascal code:

procedure write_natural (var target: text; scribendum: natural);
var
  digits: stack;
    { It is easiest to recover the digits of the number to be printed
      least significant first.  Consequently, we push them onto this
      stack as they are recovered, then pop them off this stack to
      print them one by one. }
  working: natural;
    { a throwaway copy of scribendum; later, the part of the number not
      yet reduced to decimal digits }
  quot: natural;
    { the quotient resulting from a division of working by 10 }
  rem: base_digits;
    { the remainder resulting from such a division }
begin
  if is_zero (scribendum) then
    write (target, '0')
  else begin
    initialize_stack (digits);
    assign (working, scribendum);
    while not is_zero (working) do begin
      divide_natural_by_digit (working, 10, quot, rem);
      deallocate_natural (working);
      working := quot;
      push_to_stack (rem, digits)
    end;
    deallocate_natural (working);
    while not is_empty_stack (digits) do
      write (target, pop_stack (digits) : 1)
  end
end;

Conversion to and from strings

function string_to_natural (var s: string): natural;
var
  result: natural;
    { the natural number expressed by s, constructed by evaluating s
      digit by digit }
  position: integer;
    { counts off the character positions in s }
begin
  allocate_zero (result);
  for position := 1 to strlen (s) do begin
    scale_up_by_digit (result, 10);
    increment_by_digit (result, ord (s[position]) - ORD_ZERO)
  end;
  string_to_natural := result
end;

procedure natural_to_string (n: natural; var s: string);
var
  digits: stack;
    { It is easiest to recover the digits of the number least significant
      first.  Consequently, we push them onto this stack as they are
      recovered, then pop them off this stack to store them from left
      to right in s. }
  working: natural;
    { a throwaway copy of n; later, the part of the number not
      yet reduced to decimal digits }
  quot: natural;
    { the quotient resulting from a division of working by 10 }
  rem: base_digits;
    { the remainder resulting from such a division }
  digit_count: integer;
    { tallies the number of digits in the decimal numeral for n }
  position: integer;
    { counts off the character positions in s }
begin
  if is_zero (n) then begin
    setstrlen (s, 1);
    s[1] := '0'
  end
  else begin
    initialize_stack (digits);
    digit_count := 0;
    assign (working, n);
    while not is_zero (working) do begin
      divide_natural_by_digit (working, 10, quot, rem);
      deallocate_natural (working);
      working := quot;
      push_to_stack (rem, digits);
      digit_count := digit_count + 1
    end;
    deallocate_natural (working);
    setstrlen (s, digit_count);
    for position := 1 to digit_count do
      s[position] := chr (ORD_ZERO + pop_stack (digits));
  end
end;

The assign procedure

procedure assign (var target: natural; source: natural);
var
  traverser: link;
    { a pointer to each successive digit of source }
begin
  allocate_zero (target);
  traverser := source^.highest;
  while traverser <> NIL do begin
    append (traverser^.digit, target);
    traverser := traverser^.lower
  end
end;

Transfer functions

function integer_to_natural (n: integer): natural;
var
  result: natural;
    { the natural number of value equal to n, constructed digit by digit
      through successive divisions of n }
begin
  allocate_zero (result);
  while n <> 0 do begin
    prepend (n mod BASE, result);
    n := n div BASE
  end;
  integer_to_natural := result
end;

function natural_to_integer (n: natural): integer;
var
  result: integer;
    { the integer of value equal to n, constructed by evaluating n
      digit by digit }
  traverser: link;
    { a pointer to each successive digit of n, from most significant to
      least significant }
begin
  result := 0;
  traverser := n^.highest;
  while traverser <> NIL do begin
    result := result * BASE + traverser^.digit;
    traverser := traverser^.lower
  end;
  natural_to_integer := result
end;

Additions

(2) Although the library procedures are themselves careful to recycle any storage they allocate for natural numbers that are used only internally, we have not yet provided any way for the calling program to declare that it has finished with a natural number and to deallocate the storage it occupies.

Ideally, the run-time system that is linked to executable programs would provide a garbage collector, that is, a routine that would keep track of allocated storage and recycle any pieces of it that the program was no longer able to access. Many programming languages (such as Common Lisp, Icon, and Java) have a garbage collector, but Pascal does not. The alternative is to add the deallocate_natural procedure to the interface, requiring (or rather inviting) programmers to invoke it every time they want to dispose of a natural number. Unfortunately, this makes many computations much clumsier to write. For example, if one wants to compute and display the value of ax + by, where a, x, b and y are natural numbers, one would like to be able to write

output_natural (sum (product (a, x), product (b, y)))

temp1 := product (a, x);
temp2 := product (b, y);
temp3 := sum (temp1, temp2);
deallocate_natural (temp1);
deallocate_natural (temp2);
output_natural (temp3);
deallocate_natural (temp3)

With this introduction, you're now ready to read the naturals module, in which all the procedures and functions are spelled out completely. The major difference in the code is that the preconditions of the various procedures and functions, which are simply taken for granted in this handout, are carefully tested, and errors are reported when those preconditions are violated.

This document is available on the World Wide Web as

http://www.math.grin.edu/~stone/courses/fundamentals/bignums.html