### Natural numbers

Let's return now to the limitations of Pascal's `Integer` data type. When we initially considered integers as an abstract data type, our implementation was hampered by the fact that integers larger than `MaxInt` were not representable. Now that we have a variety of pointer structures -- specifically, a module for bidirectional lists -- we can easily get around this limitation. Unsigned integers of unlimited size that are implemented as dynamically allocated structures are often called bignums; the main objective of this handout is to indicate how the usual arithmetic operations can be performed on bignums.

We begin by considering natural numbers, which are unsigned integers beginning with 0. 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 bidirectional list of digits is 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, 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 cube root of `MaxInt + 1`, 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.

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

```const
BaseOfNumeration = 1290;
{ = largest even integer less than the cube root of MaxInt + 1 }
type
Element = 0 .. BaseOfNumeration - 1;
Natural = BidirectionalList; { of Element }
```
Since our Naturals module will be implemented in a language that already has a data type that includes some natural numbers, we should add transfer functions to the interface so that we can convert between our data type and the one that is built in:

```function PascalIntegerToNatural (N: Integer): Natural;
function NaturalToPascalInteger (N: Natural): Integer;
```
Of course, there are preconditions for both of these functions: `PascalIntegerToNatural` requires a non-negative argument, and `NaturalToPascalInteger` requires an argument less than or equal to `MaxInt`.

Although the module 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 a `DeallocateNatural` 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

```WriteNatural (Output, AddNatural (MultiplyNatural (A, X),
MultiplyNatural (B, Y)))
```
But without automatic garbage collection, one must instead write

```Temp1 := MultiplyNatural (A, X);
Temp2 := MultiplyNatural (B, Y);
DeallocateNatural (Temp1);
DeallocateNatural (Temp2);
WriteNatural (Output, Temp3);
DeallocateNatural (Temp3)
```
This is a strong argument for garbage collection.

Finally, a somewhat subtle problem arises from the fact that this implementation uses a data structure accessed through a pointer for a type that many programmers don't think of as having an internal structure. The application programmer will have to be aware of the fact that the the operators `=`, `<>`, and `:=`, though applicable to the `Natural` type, do not have their accustomed meanings. The `=` operator, applied to two values of type `Natural`, determines whether they are equal as pointers -- that is, whether the same storage location is accessed through both -- not whether the values stored there are numerically equal. Similarly, an assignment using `:=` will copy a pointer instead of building a fresh copy of the number. The difference is usually invisible, but it appears when the number is deallocated; application programmers will be disconcerted to find that calling `DeallocateNumber` to recycle the value stored in one variable may also invalidate other variables -- those containing copies of the pointer.

To take the place of the `:=` operator for natural numbers, therefore, we supply a procedure that performs assignment by copying:

```procedure AssignNatural (var Target: Natural; Source: Natural);
```
This should be invoked whenever a freshly allocated copy of an existing natural number is needed. Similarly, the application programmer should use the `EqualNaturals` and `UnequalNaturals` functions (instead of `=` and `<>`) when comparingz values.

For the same reason, the documentation for the `Naturals` module should carefully distinguish the operations that always return a freshly allocated natural number (such as add) from those that always return one of their arguments unchanged (such as major); it would be a serious design mistake for the implementation to do sometimes the one thing and sometimes the other, even though such an approach might be tempting for reasons of efficiency. (For example, if the implementation of the add operation determines that the addend is zero, it should nevertheless build and return a fresh copy of the augend rather than just handing back a pointer to the augend. Since add must usually return a freshly allocated number in any case, it should always do so.) To keep the bidirectional 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 an empty bidirectional list.

### Implementing the functions and procedures

Some of the functions and procedures present interesting implementation problems. The algorithms for working with large natural numbers digit by digit are often similar to those taught in elementary school for hand computation using decimal numeration. (It does not follow that they are easy to program.)

The basic idea of the addition algorithm is to start at the low end of both operands and to add corresponding digits together, also adding in any carry that is generated from the next lower digit position. If the result is less than the base, put it in the corresponding position in the sum and carry 0; otherwise, subtract the base from it, put the difference in the corresponding position in the sum, and carry 1. Move to the next higher position and repeat. When you run out of digits in either operand, fill in with zeroes; when you run out of digits in both operands, stop as soon as the carry is 0.

The `MakeZero` function that is invoked to initialize the sum of the addition simply calls `EmptyBidirectionalList` to get a bidirectional list representing the number zero. It is one of several functions and procedures that are not exported from the module, but are present only to simplify the writing of the functions and procedures that will be exported.

Notice that the function checks whether `Augend` and `Addend` are identical before trying to manipulate them as bidirectional lists; if it turns out that the two are identical, a different algorithm is used for the computation. The problem is that the algorithm described above manipulates the cursors of `Augend` and `Addend` as if they were independent, which is not the case if `Augend` and `Addend` are equal as pointers.

#### Subtraction

The subtraction algorithm has the same general structure as the addition algorithm. A borrow from the next higher column is generated whenever the current digit of the minuend is less than the current digit of the subtrahend (or even if they are equal, if the previous digit required a borrow), and `BaseOfNumeration` is then added to the minuend digit to allow the subtraction to be performed. The process stops when all the digits of the minuend have been examined.

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 `StripLeadingZeroes` is invoked to remove any zero digits from the result before it is returned.

For natural numbers, subtraction is defined only when the minuend is greater than or equal to the subtrahend, and this procedure uses an assertion to check that this precondition is satisfied. An alternative design would be to return zero (in which case the operation is sometimes called monus) or the absolute value of the difference (disparity) in such a case.

#### Multiplication

The multiplication algorithm is the same one that is used in pencil-and-paper computation with decimal numerals: Starting from the low end of the multiplier, take one digit at a time and compute the partial product of the multiplicand and that digit; shift the partial product upwards by the number of places that come after the current digit of the multiplier. The sum of the shifted partial products is the complete product. The only difference is that on the computer it is simpler to add each shifted partial product to a running total than to save all the additions for the end.

A zero multiplicand is treated as a special case, so that a lot of generation and shifting of partial products that are all equal to 0 is avoided. The generation of a partial product is also skipped when a zero digit is encountered in the multiplier.

This function invokes the utility function `MultiplyNaturalByDigit`, 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.

#### Division

This is the tricky one. The implementation presented here closely follows the algorithm that Per Brinch Hansen carefully develops and proves correct in ``Multiple-length division: a tour of the minefield'' (Software -- practice and experience 24 [1994], 579--601).

There are four special cases that can be handled without any hard work: (1) Division by zero should be trapped and dismissed as an error. (2) If the dividend and the divisor are identical, the quotient is 1 and the remainder is 0; return these values and avoid doing any cursor manipulations! (3) When the dividend is less than the divisor, the quotient is 0 and the remainder is a copy of the dividend. (4) When the divisor contains only one digit, we can use the much simpler short-division algorithm provided by the (unexported) `DivideNaturalByDigit` procedure.

The `DivideNaturalByDigit` procedure initializes the quotient to zero and gets ready to recover its first digit. The quotient will be one digit shorter than the dividend if, and only if, the first digit of the dividend is less than the divisor, so the procedure tests to see whether this is the case; if so, it suppresses the leading zero of the quotient, arranges for the remainder of the first-digit division to be equal to the first digit of the dividend, and then enters the main division loop with the dividend's cursor pointing to its second digit. This ensures that the first pass through the loop will generate a non-zero leading digit for the quotient.

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 the quotient, and changes `Remainder` 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 DivideNatural (Dividend, Divisor: Natural;
var Quotient, Remainder: Natural);
begin
Assert (not ZeroNatural (Divisor), DivisionException,
NaturalExceptionHandler);
if Dividend = Divisor then begin
Quotient := PascalIntegerToNatural (1);
Remainder := MakeZero
end
else if LessNatural (Dividend, Divisor) then begin
Quotient := MakeZero;
Remainder := CopyBidirectionalList (Dividend)
end
else if LengthOfBidirectionalList (Divisor) <= 1 then begin
DivideNaturalByDigit (Dividend, FirstOfBidirectionalList (Divisor),
Quotient, Irem);
Remainder := PascalIntegerToNatural (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 `TrialDigit` to `BaseOfNumeration - 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 `BaseOfNumeration` is 2, but it's terribly inefficient when `BaseOfNumeration` 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 `FirstThree` (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 `FirstTwo`, then dividing `FirstThree` by `FirstTwo` 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 `FirstTwo` increases, since the relative error involved in discarding the subsequent digits becomes less. A convenient way to make `FirstTwo` 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 ```BaseOfNumeration 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 `BaseOfNumeration`, in the rare case that the correct digit is `BaseOfNumeration - 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.

Thus the plan for the really hard case looks like this: (1) Compute the appropriate scaling factor. (2) Attach a leading zero to the dividend and scale it up; scale up the divisor. (3) Initialize the variable (call it `Residue`) that will hold the appropriate number of digits of the scaled dividend by copying into it one more digit than the scaled divisor has. (4) Call `FindTrialDigit` to obtain a good estimate of the next digit of the quotient. (5) 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. (6) 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`. (7) If there are any remaining digits in the scaled dividend, attach the first of them at the low end of the residue and repeat steps 4-7. (8) Otherwise, 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.

#### Exponentiation

To raise a number to a power, one approach would be to use repeated multiplication by the base, but this is unpleasantly slow when the exponent is large. To reduce the number of intermediate results, we can use two laws of exponents: y^(2x) = (y^2)^x and y^(2x + 1) = y x (y^2)^x. Each of these laws shows a way to compute the result of raising a given base y to a large power by computing a larger base y^2 raised to a smaller power. (We use the first law to reduce even exponents, the second to reduce odd ones.) By repeating the process, we can always reduce the exponent to 0, at which point the result is known to be 1 regardless of the value of the base. This suggests the following recursive implementation:

```function RaiseNatural (Base, Exponent: Natural): Natural;
begin
if ZeroNatural (Exponent) then
RaiseNatural := PascalIntegerToNatural (1)
else if EvenNatural (Exponent) then
RaiseNatural := RaiseNatural (SquareNatural (Base), Half (Exponent))
else
RaiseNatural := MultiplyNatural (Base,
RaiseNatural (SquareNatural (Base),
Half (Exponent)))
end;
```
As soon as one considers how to recycle the storage used by this function, however, it becomes apparent that using recursion is going to require a lot of storage if either the base or the exponent is at all large, and the recursion aggravates the problem by retaining all of the storage for intermediate results until the last recursive call is reached. The key difficulty is that the multiplications performed in the various calls to `MultiplyNatural` cannot be performed until after the recursive calls to `RaiseNatural` have produced the multiplier. However, since multiplication is commutative, it actually would not matter whether these multiplications were done in reverse order, as the recursive calls are exiting, or in forwards order, as the exponent is being reduced; in either case, the current value of the base will be a factor of the product if, and only if, the current exponent is odd. But if the multiplications are performed on the way into each recursive call, then the multiplicand can be discarded -- or, better yet, overwritten -- as soon as its value has been used for the last time.

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 RaiseNatural (Base, Exponent: Natural): Natural;

function RaiseHelper (Base, Exponent, Accumulator: Natural): Natural;
begin
if ZeroNatural (Exponent) then
RaiseHelper := Accumulator
else begin
if OddNatural (Exponent) then
Accumulator := MultiplyNatural (Accumulator, Base);
Base := SquareNatural (Base);
Exponent := Half (Exponent);
RaiseHelper := RaiseHelper (Base, Exponent, Accumulator)
end
end;

begin { function RaiseNatural }
RaiseNatural := RaiseHelper (Base, Exponent, PascalIntegerToNatural (1))
end;
```
Note now that the inner function `RaiseHelper` is ``tail-recursive'': When it receives the value of a recursive call, it simply returns it, without performing any further computation on it. It is particularly simple to convert a tail-recursive function to an iterative form, and in Pascal the resulting function will generally run faster and consume less memory than the recursive version. The code in the completed `Naturals` module shows the iterative form of this algorithm. All of these versions of `RaiseNatural` return a result of 1 when given a base of 0 and an exponent of 0. This accords with the usual combinatorial interpretation of the expression 0^0. Some mathematicians regard this value as undefined, since the limiting value of y^0 as y approaches 0 is 1, while the limiting value of 0^x as x approaches 0 is 0; but this line of reasoning, based on the continuity of the real numbers, seems rather less forceful when we are considering only natural numbers. It would be easy enough to treat this special case as an error, but there isn't any harm in letting it pass, so long as the function's behavior is well documented.

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

In the `SuccessorOfNatural` function, the idea is to copy the operand digit by digit, starting at the low end and propagating a carry exactly as far as necessary. If the carry is propagated all the way through the operand, an extra place must be added to the result.

The algorithm for the `PredecessorOfNatural` function is similar, but it differs in three important ways: (1) It must test the precondition that the given natural number is non-zero; (2) a borrow is propagated instead of a carry; and (3) instead of testing whether a new leading digit must be added, the function must check whether the leading digit of the result has been reduced to zero and remove it if it has.

#### Twice

Finding the double of a given natural number is quite similar to multiplying a natural number by a digit (as in the procedure `MultiplyNaturalByDigit`, above), except for a few small optimizations resulting from the fact that the carry is known to be either zero or one.

#### Square and cube

The most straightforward way to do these is simply to multiply.

#### Comparison functions

Since the header node for a natural number includes a count of the number of digits in that number (that is, the number of components in the doubly-linked list), one can make a first cut at comparing two numbers by comparing their digit counts, proceeding to examine the digit values only in the special case where the two numbers are equal in length; this considerably simplifies the coding of the traversals of the doubly-linked lists.

For example, in the `LessNatural` 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.

The equality test is somewhat similar, but simpler, since when an inequality is detected it is not necessary to determine which operand is smaller.

#### Zero

In the special case of comparison with zero, there is an easy optimization: A natural number is equal to zero if, and only if, the bidirectional list is empty.

#### Major and minor

These functions use the obvious algorithm. Note, however, that each of them returns one of the operands itself, not a copy.

#### Multiple

The basic idea is to determine whether dividing the first operand by the second leaves a remainder of zero. The special case in which the second operand is zero is split off and handled separately before the division is attempted.

Once again, note that this function is careful to recycle the storage associated with the natural number that it generates. There is a temptation to write the `else`-clause simply as

```MultipleNatural := ZeroNatural (RemainderNatural (Candidate, Unit))
```
but this would create a memory leak. The storage allocated for the remainder could never be deallocated, because no pointer to it is retained.

#### Parity functions

Since the base of numeration is even, one can tell whether a given natural number is even or odd by looking at its last digit. (The other digits signify multiples of positive powers of the base and therefore contribute only even values to the number represented; only the last digit can contribute an even value or an odd one.) Note, however, that since zero has no last digit (in our implementation) it must be handled as a separate case.

#### Input and output procedures

The `ReadNatural` procedure reads in decimal digits one by one from a text file, continuing until either the end of the file or a non-digit character is reached.

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 `ReadNatural` procedure also consumes such characters, by invoking the `SkipWhiteSpace` procedure.

The main difficulty with the `WriteNatural` 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 `BaseDigits`, 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 `WriteNatural` procedure uses five routines from the `Stacks` module: `CreateStack`, which builds and returns an empty stack; `PushToStack`, which stores a given integer value in a given stack; `PopFromStack`, which recovers an integer value from a given stack and returns it; and `EmptyStack`, which determines whether a given stack is empty; and `DeallocateStack`, which retires the stack when the procedure is finished with it.

The plan of the `WriteNatural` 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.

#### Assign

The procedure for assigning a natural number by copying it into new storage simply invokes the `CopyBidirectionalList` procedure.

#### Transfer functions

To convert an integer to a natural number, repeatedly divide it by `BaseOfNumeration` and add each successive remainder at the front of the natural number.

To convert a natural number to an integer, initialize a result variable to zero, then take up its digits one at a time, from most significant to least significant; at each digit, multiply the current result by `BaseOfNumeration` and add the value of the current digit.

While testing and debugging the `Naturals` module, I found it helpful to be able to invoke a procedure (`DebugOutputNatural`) that displayed the internal structure of natural numbers. There is an argument for including this procedure in the public interface, even though in a way it breaks the abstraction boundary, because it might allow the user of the library to supply more informative bug reports.

### The `Naturals` module

```{ This module defines an interface for a natural-number data type with no
upper bound and implements it for HP 9000 Series 700 workstations under
HP-UX 9.x, using HP Pascal.

Programmer: John Stone, Grinnell College.
Original version: January 10-24, 1996.
Last revised: November 5-6, 1996.
}

\$heap_dispose on\$

module Naturals;

\$search 'natural-elements.o, bidirectional-lists.o'\$

import
Elements, BidirectionalLists;

export

type
Natural = BidirectionalList;

{ The AddNatural function adds any two natural numbers and returns their
sum. }

{ The SubtractNatural function subtracts one natural number (the
subtrahend) from another (the minuend) and returns their difference.
It presupposes that the subtrahend is less than or equal to the
minuend. }

function SubtractNatural (Minuend, Subtrahend: Natural): Natural;

{ The MultiplyNatural function multiplies one natural number (the
multiplicand) by another (the multiplier) and returns their product. }

function MultiplyNatural (Multiplicand, Multiplier: Natural): Natural;

{ The DivideNatural procedure divides one natural number (the dividend)
by another (the divisor), returning both the quotient and the remainder
as natural numbers.  It presupposes that the divisor is not zero. }

procedure DivideNatural (Dividend, Divisor: Natural;
var Quotient, Remainder: Natural);

{ The QuotientNatural function divides one natural number (the dividend)
by another (the divisor) and returns the quotient as a natural number,
discarding the remainder.  It presupposes that the divisor is not
zero. }

function QuotientNatural (Dividend, Divisor: Natural): Natural;

{ The RemainderNatural function divides one natural number (the dividend)
by another (the divisor) and returns the remainder as a natural number,
discarding the quotient.  It presupposes that the divisor is not
zero. }

function RemainderNatural (Dividend, Divisor: Natural): Natural;

{ The RaiseNatural function raises one natural number (the base) to the
power specified by another natural number (the exponent) and returns
the result as a natural number.  It returns 1 whenever the exponent is
zero, even if the base is also zero. }

function RaiseNatural (Base, Exponent: Natural): Natural;

{ The SuccessorOfNatural function returns the natural number that
immediately follows a given natural number. }

function SuccessorOfNatural (Operand: Natural): Natural;

{ The PredecessorOfNatural function returns the natural number that
immediately precedes a given natural number.  It presupposes that its
operand is not zero. }

function PredecessorOfNatural (Operand: Natural): Natural;

{ Given a natural number, the TwiceNatural function computes and returns
its double -- that is, the natural number that is twice as large. }

function TwiceNatural (Operand: Natural): Natural;

{ Given a natural number, the SquareNatural function computes and returns
its square. }

function SquareNatural (Operand: Natural): Natural;

{ Given a natural number, the CubeNatural function computes and returns
its cube. }

function CubeNatural (Operand: Natural): Natural;

{ The EqualNaturals function determines whether its arguments are
numerically equal -- not necessarily identical as storage structures,
but equal in value. }

function EqualNaturals (LeftOperand, RightOperand: Natural): Boolean;

{ The UnequalNaturals function determines whether the first of its
arguments differs in value from its second -- not whether they differ
as storage structures, but whether their numerical values differ. }

function UnequalNaturals (LeftOperand, RightOperand: Natural): Boolean;

{ The LessNatural function determines whether the first of its arguments
is numerically less than the second. }

function LessNatural (LeftOperand, RightOperand: Natural): Boolean;

{ The LessOrEqualNatural function returns True if its first argument is
numerically less than or equal to its second, False otherwise. }

function LessOrEqualNatural (LeftOperand, RightOperand: Natural):
Boolean;

{ The GreaterNatural function determines whether the first of its
arguments is numerically greater than the second. }

function GreaterNatural (LeftOperand, RightOperand: Natural): Boolean;

{ The GreaterOrEqualNatural returns True if its first argument is
numerically greater than or equal to its second, False otherwise. }

function GreaterOrEqualNatural (LeftOperand, RightOperand: Natural):
Boolean;

{ The MajorNatural function returns the greater of its two arguments; if
they are equal, it returns the first of the two. }

function MajorNatural (LeftOperand, RightOperand: Natural): Natural;

{ The MinorNatural function returns the lesser of its two arguments; if
they are equal, it returns the first of the two. }

function MinorNatural (LeftOperand, RightOperand: Natural): Natural;

{ The ZeroNatural function determines whether a given natural number is
0, returning True if it is and False if it is not. }

function ZeroNatural (Operand: Natural): Boolean;

{ The MultipleNatural function determines whether its first argument is
an exact multiple of its second.  (The only exact multiple of 0 is
0.) }

function MultipleNatural (Candidate, Unit: Natural): Boolean;

{ The EvenNatural function determines whether its argument is even --
that is, whether it is an integer multiple of 2. }

function EvenNatural (Operand: Natural): Boolean;

{ The OddNatural function determines whether its argument is odd -- that
is, whether dividing it by 2 leaves a remainder of 1. }

function OddNatural (Operand: Natural): Boolean;

{ The ReadNatural procedure collects a sequence of decimal digits,
optionally preceded by any number of whitespace characters, from a
specified text file and stores the natural number that they represent
into a given variable.  It presupposes that the text file has already
been opened for input. }

procedure ReadNatural (var Source: Text; var Legend: Natural;
var Success: Boolean);

{ The WriteNatural procedure writes the decimal numeral for a given
natural number to a specified text file, with no leading zeroes or
spaces. (But the single digit '0' is written if the given natural
number is 0.)  It presupposes that the text file has already been
opened for output. }

procedure WriteNatural (var Target: Text; Scribend: Natural);

{ The AssignNatural procedure creates a natural number equal in value to
its second argument, but in separately allocated storage, and returns
it through its first argument. }

procedure AssignNatural (var Target: Natural; Source: Natural);

{ Given any non-negative integer argument, the PascalIntegerToNatural
function constructs and returns a natural number of equal value. }

function PascalIntegerToNatural (N: Integer): Natural;

{ Given any natural number not exceeding MaxInt, the
NaturalToPascalInteger function returns the integer of equal value. }

function NaturalToPascalInteger (N: Natural): Integer;

{ The DebugOutputNatural procedure writes out, to standard error, a
representation of a natural number that indicates its internal
structure.  This procedure is provided for debugging purposes only. }

procedure DebugOutputNatural (var Target: Text; Scribend: Natural);

{ The DeallocateNatural procedure frees all of the storage allocated for
a given natural number and changes the value of its argument to nil.
It presupposes that a natural number has been stored in N and not
previously deallocated. }

procedure DeallocateNatural (var N: Natural);

implement

\$search 'natural-elements.o, stacks.o'\$
import
StdErr, Stacks;

const
{ The following constants are more or less arbitrary integers
signifying various kinds of exceptions that can occur within this
module. }

FirstExceptionCode = 1;
NegativeException = 1;
IntegerOverflowException = 2;
SubtractionException = 3;
DivisionException = 4;
PredecessorException = 5;
ExceptionException = 6;
LastExceptionCode = ExceptionException;

type
Bit = 0 .. 1;
BaseDigits = 0 .. BaseOfNumeration - 1;
{ possible values of one digit in the chosen base of numeration }

{ The NaturalExceptionHandler procedure is invoked whenever one of the
preconditions for the successful execution of a procedure is found to
be false.  It prints out an appropriate explanation of the exception
just before the program is halted. }

procedure NaturalExceptionHandler (ExceptionCode: Integer);
begin
if (ExceptionCode < FirstExceptionCode) or
(LastExceptionCode < ExceptionCode) then
ExceptionCode := ExceptionException;
Write (StdErr, 'Exception #', ExceptionCode : 1,
' in module Naturals: ');
case ExceptionCode of
NegativeException:
WriteLn (StdErr, 'negative number as argument to function ',
'PascalIntegerToNatural');
IntegerOverflowException:
WriteLn (StdErr, 'value of numeral exceeds MaxInt in function ',
'NaturalToPascalInteger');
SubtractionException:
WriteLn (StdErr, 'minuend less than subtrahend in function ',
'SubtractNatural');
DivisionException:
WriteLn (StdErr, 'division by zero attempted in procedure ',
'DivideNatural');
PredecessorException:
WriteLn (StdErr, 'zero operand in function PredecessorOfNatural');
ExceptionException:
WriteLn (StdErr, 'argument out of range in procedure ',
'NaturalExceptionHandler.');
end
end;

{ The MakeZero function creates and returns a newly allocated instance
of the natural number 0. }

function MakeZero: Natural;
begin
MakeZero := MakeEmptyBidirectionalList
end;

{ The StripLeadingZeroes procedure removes any zero digits from the
beginning of the digit list representing a natural number. }

var
Continue: Boolean;
{ indicates whether the search for a non-zero digit can and should
continue }
begin
Continue := True;
while Continue do
if EmptyBidirectionalList (N) then
Continue := False
else if FirstOfBidirectionalList (N) = 0 then
DeleteFirstOfBidirectionalList (N)
else
Continue := False
end;

{ The MultiplyNaturalByDigit function finds the product of a given
natural number (the multiplicand) and a one-digit number in the bignum
base of numeration and returns that product as a natural number. }

function MultiplyNaturalByDigit (Multiplicand: Natural;
Multiplier: BaseDigits): Natural;
var
Result: Natural;
{ the product, as it is built up one digit at a time }
Carry: BaseDigits;
{ the number of multiples of BaseOfNumeration to be carried into the
next higher digit position }
PlaceProduct: integer;
{ the product of one digit of the multiplicand and the multiplier,
plus any Carry from the next lower digit position }
begin
Result := MakeZero;
if Multiplier <> 0 then begin
Carry := 0;
CursorToEndOfBidirectionalList (Multiplicand);
while not NullCursorInBidirectionalList (Multiplicand) do begin
PlaceProduct := Carry +
ElementAtCursorInBidirectionalList (Multiplicand) * Multiplier;
PrependToBidirectionalList (PlaceProduct mod BaseOfNumeration, Result);
Carry := PlaceProduct div BaseOfNumeration;
RetreatCursorAlongBidirectionalList (Multiplicand)
end;
if Carry <> 0 then
PrependToBidirectionalList (Carry, Result)
end;
MultiplyNaturalByDigit := Result
end;

{ The ScaleUpByDigit procedure scales up the value of a given
natural-number variable by a given factor. }

procedure ScaleUpByDigit (var Given: Natural; Scale: BaseDigits);
var
Carry: BaseDigits;
{ the number of multiples of BaseOfNumeration to be carried into the
next higher digit position }
PlaceProduct: 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
DeallocateNatural (Given);
Given := MakeZero
end
else begin
Carry := 0;
CursorToEndOfBidirectionalList (Given);
while not NullCursorInBidirectionalList (Given) do begin
PlaceProduct := Carry +
ElementAtCursorInBidirectionalList (Given) * Scale;
AssignAtCursorInBidirectionalList (Given,
PlaceProduct mod BaseOfNumeration);
Carry := PlaceProduct div BaseOfNumeration;
RetreatCursorAlongBidirectionalList (Given)
end;
if Carry <> 0 then
PrependToBidirectionalList (Carry, Given)
end
end;

{ The DivideNaturalByDigit procedure determines the quotient (as a
natural number) and the remainder (as an integer) when a given natural
number is divided by a divisor less than BaseOfNumeration. }

procedure DivideNaturalByDigit (Dividend: Natural;
Divisor: BaseDigits; var Quotient: Natural; var Remainder: BaseDigits);
var
ShortDividend: Integer;
{ a partial dividend, composed of the remainder left over from the
division in the previous position, scaled up by a factor of
BaseOfNumeration, plus the digit in the current position }
NextDigit: BaseDigits;
{ one digit at a time from the dividend }
begin

{ Prevent division by zero. }

Assert (Divisor <> 0, DivisionException, NaturalExceptionHandler);

{ Set up to recover the first digit of the quotient.  The quotient will
be one digit shorter than the dividend if, and only if, the first
digit of the dividend is less than the divisor, so test for this
case; if it is found, suppress the leading zero of the quotient,
arrange for the remainder of the first-digit division to be equal
to the first digit of the divisor, and then enter the main division
loop with traverser pointing to the second digit of the divisor.
This will ensure that the first pass through the loop will generate
a non-zero leading digit for the quotient. }

Quotient := MakeZero;
CursorToStartOfBidirectionalList (Dividend);
if NullCursorInBidirectionalList (Dividend) then
Remainder := 0
else begin
NextDigit := ElementAtCursorInBidirectionalList (Dividend);
if NextDigit < Divisor then begin
Remainder := NextDigit;
end
else
Remainder := 0
end;

{ On each iteration of the loop, construct the "short dividend",
perform short division, attach the quotient (which will always be
a single BaseDigit) to the right end of the quotient, and change
Remainder to reflect the remainder of the division. }

while not NullCursorInBidirectionalList (Dividend) do begin
ShortDividend := Remainder * BaseOfNumeration +
ElementAtCursorInBidirectionalList (Dividend);
AppendToBidirectionalList (Quotient, ShortDividend div Divisor);
Remainder := ShortDividend mod Divisor;
end

end;

{*** Exported procedures and functions ***}

var
Result: Natural;
{ the sum of the arguments, as it is constructed }
Carry: Bit;
{ a carry from one digit position to the next }
{ digits from corresponding positions in the augend and addend }
PlaceSum: 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
else begin
Result := MakeZero;
CursorToEndOfBidirectionalList (Augend);
Carry := 0;
while not NullCursorInBidirectionalList (Augend) or
not NullCursorInBidirectionalList (Addend) or (Carry <> 0) do begin
if NullCursorInBidirectionalList (Augend) then
AugendDigit := 0
else begin
AugendDigit := ElementAtCursorInBidirectionalList (Augend);
RetreatCursorAlongBidirectionalList (Augend)
end;
else begin
end;
PlaceSum := AugendDigit + AddendDigit + Carry;
if PlaceSum < BaseOfNumeration then begin
PrependToBidirectionalList (PlaceSum, Result);
Carry := 0
end
else begin
PrependToBidirectionalList (PlaceSum - BaseOfNumeration, Result);
Carry := 1
end
end;
end
end;

function SubtractNatural (Minuend, Subtrahend: Natural): Natural;
var
Result: Natural;
{ the difference of the arguments, as it is constructed }
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 }
MinuendDigit, SubtrahendDigit: BaseDigits;
{ digits from corresponding positions in the minuend and subtrahend }
begin
Result := MakeZero;
if Minuend <> Subtrahend then begin
CursorToEndOfBidirectionalList (Minuend);
CursorToEndOfBidirectionalList (Subtrahend);
Borrow := 0;
while not NullCursorInBidirectionalList (Minuend) do begin
MinuendDigit := ElementAtCursorInBidirectionalList (Minuend);
RetreatCursorAlongBidirectionalList (Minuend);
if NullCursorInBidirectionalList (Subtrahend) then
SubtrahendDigit := 0
else begin
SubtrahendDigit := ElementAtCursorInBidirectionalList (Subtrahend);
RetreatCursorAlongBidirectionalList (Subtrahend);
end;
if MinuendDigit - Borrow < SubtrahendDigit then begin
PrependToBidirectionalList (BaseOfNumeration + MinuendDigit -
Borrow - SubtrahendDigit, Result);
Borrow := 1
end
else begin
PrependToBidirectionalList (MinuendDigit - Borrow -
SubtrahendDigit, Result);
Borrow := 0
end
end;
Assert (NullCursorInBidirectionalList (Subtrahend) and (Borrow = 0),
SubtractionException, NaturalExceptionHandler);
end;
SubtractNatural := Result
end;

function MultiplyNatural (Multiplicand, Multiplier: Natural): Natural;
var
DuplicateFlag: Boolean;
{ indicates whether the arguments are identical; if so, one must
be duplicated to avoid conflicting cursor movements }
Result: Natural;
{ the product, as it is built up from partial products }
LowerPositions: Integer;
{ the number of digit positions below the current one }
MultiplierDigit: BaseDigits;
{ one digit extracted from the multiplier }
Partial: Natural;
{ a partial product: the product of the multiplicand and one digit
of the multiplier }
ZeroTallier: Integer;
{ counts off trailing zeroes as they are appended to a partial
product }

procedure IncrementNatural (var Given: Natural; Amount: Natural);
var
Carry: Bit;
{ a carry from one digit position to the next }
GivenDigit, AmountDigit: BaseDigits;
{ digits from corresponding positions in Given and Amount }
PlaceSum: Integer;
{ the sum of the digits in the same position in Given and Amount,
plus any Carry from the next lower digit position }
begin
CursorToEndOfBidirectionalList (Given);
CursorToEndOfBidirectionalList (Amount);
Carry := 0;
while not NullCursorInBidirectionalList (Given) or
not NullCursorInBidirectionalList (Amount) or
(Carry <> 0) do begin
if NullCursorInBidirectionalList (Given) then
GivenDigit := 0
else
GivenDigit := ElementAtCursorInBidirectionalList (Given);
if NullCursorInBidirectionalList (Amount) then
AmountDigit := 0
else begin
AmountDigit := ElementAtCursorInBidirectionalList (Amount);
RetreatCursorAlongBidirectionalList (Amount)
end;
PlaceSum := GivenDigit + AmountDigit + Carry;
if PlaceSum < BaseOfNumeration then
Carry := 0
else begin
PlaceSum := PlaceSum - BaseOfNumeration;
Carry := 1
end;
if NullCursorInBidirectionalList (Given) then
PrependToBidirectionalList (PlaceSum, Given)
else begin
AssignAtCursorInBidirectionalList (Given, PlaceSum);
RetreatCursorAlongBidirectionalList (Given)
end
end
end;

begin { procedure MultiplyNatural }
Result := MakeZero;
if not ZeroNatural (Multiplicand) then begin
DuplicateFlag := (Multiplicand = Multiplier);
if DuplicateFlag then
Multiplier := CopyBidirectionalList (Multiplicand);
LowerPositions := 0;
CursorToEndOfBidirectionalList (Multiplier);
while not NullCursorInBidirectionalList (Multiplier) do begin
MultiplierDigit := ElementAtCursorInBidirectionalList (Multiplier);
if MultiplierDigit <> 0 then begin
Partial :=
MultiplyNaturalByDigit (Multiplicand, MultiplierDigit);
for ZeroTallier := 1 to LowerPositions do
AppendToBidirectionalList (Partial, 0);
IncrementNatural (Result, Partial);
DeallocateNatural (Partial)
end;
LowerPositions := LowerPositions + 1;
RetreatCursorAlongBidirectionalList (Multiplier)
end;
if DuplicateFlag then
DeallocateBidirectionalList (Multiplier)
end;
MultiplyNatural := Result
end;

procedure DivideNatural (Dividend, Divisor: Natural;
var Quotient, Remainder: Natural);
var
Irem: BaseDigits;
{ the remainder of a short division -- used only if divisor is a
single digit }
ScalingFactor: BaseDigits;
{ a multiplier applied both to the dividend and to the divisor, to
ensure that the FindTrialDigit function will accurately guess
the correct quotient digit as often as possible }
ScaledDividend, ScaledDivisor: Natural;
{ the values of dividend and divisor, respectively, each multiplied
by the scaling factor }
FirstTwo: Integer;
{ the value of the first two digits of the scaled divisor, considered
as a natural number in its own right; this is used in the process
of estimating the next digit of the quotient }
Position: Integer;
{ counts off digits as they are transferred from the high end of the
scaled dividend to the residue }
Residue: Natural;
{ a ``partial dividend,'' formed by successively appending digits of
the scaled dividend and subtracting away partial products formed
by multiplying digits of the quotient by the scaled divisor; zeroes
are prepended as needed to ensure that the number of digits in the
residue is always one more than in the divisor }
Finished: Boolean;
{ indicates whether the division of the scaled dividend by the scaled
divisor is complete, that is, whether we have reached the low end
of the scaled divisor }
TrialDigit: BaseDigits;
{ an estimate, possibly one too high, of the next digit of the
quotient }
TrialProduct: Natural;
{ the product of the scaled divisor and the current trial digit }
{ the remainder from the short division in which the scaling factor
is divided out of the residue; this remainder is always zero and
can therefore be ignored }

function FindTrialDigit (Residue: Natural): BaseDigits;
var
FirstThree: 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 }
Position: Integer;
{ counts off digits in Residue }
Trial: Integer;
{ an approximate value for the next quotient digit }
begin
FirstThree := 0;
CursorToStartOfBidirectionalList (Residue);
for Position := 1 to 3 do begin
FirstThree := FirstThree * BaseOfNumeration +
ElementAtCursorInBidirectionalList (Residue);
end;
Trial := FirstThree div FirstTwo;
if Trial = BaseOfNumeration then
FindTrialDigit := BaseOfNumeration - 1
else
FindTrialDigit := Trial
end;

procedure DecrementNatural (var Given: Natural; Amount: Natural);
var
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 }
GivenDigit, AmountDigit: BaseDigits;
{ digits from corresponding positions in given and amount }
begin
if Given = Amount then
Given := MakeZero
else begin
CursorToEndOfBidirectionalList (Given);
CursorToEndOfBidirectionalList (Amount);
Borrow := 0;
while not NullCursorInBidirectionalList (Given) do begin
GivenDigit := ElementAtCursorInBidirectionalList (Given);
if NullCursorInBidirectionalList (Amount) then
AmountDigit := 0
else begin
AmountDigit := ElementAtCursorInBidirectionalList (Amount);
RetreatCursorAlongBidirectionalList (Amount)
end;
if GivenDigit - Borrow < AmountDigit then begin
AssignAtCursorInBidirectionalList (Given,
BaseOfNumeration + GivenDigit - Borrow - AmountDigit);
Borrow := 1
end
else begin
AssignAtCursorInBidirectionalList (Given,
GivenDigit - Borrow - AmountDigit);
Borrow := 0
end;
RetreatCursorAlongBidirectionalList (Given)
end;
Assert (NullCursorInBidirectionalList (Amount) and (Borrow = 0),
SubtractionException, NaturalExceptionHandler);
end
end;

begin { procedure DivideNatural }

{ Trap division by zero. }

Assert (not ZeroNatural (Divisor), DivisionException,
NaturalExceptionHandler);

{ Handle the dangerous case: Dividend is identical to Divisor. }

if Dividend = Divisor then begin
Quotient := PascalIntegerToNatural (1);
Remainder := MakeZero
end

{ Handle one easy case: Dividend < Divisor. }

else if LessNatural (Dividend, Divisor) then begin
Quotient := MakeZero;
Remainder := CopyBidirectionalList (Dividend)
end

{ Handle another easy case: one-digit divisor. }

else if LengthOfBidirectionalList (Divisor) <= 1 then begin
DivideNaturalByDigit (Dividend, FirstOfBidirectionalList (Divisor),
Quotient, Irem);
Remainder := PascalIntegerToNatural (Irem)
end

{ Now for the hard case: a divisor of two or more digits and a
dividend greater than or equal to the divisor. }

else begin

{ Compute the appropriate scaling factor. }

ScalingFactor :=
BaseOfNumeration div (FirstOfBidirectionalList (Divisor) + 1);

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

ScaledDividend := CopyBidirectionalList (Dividend);
PrependToBidirectionalList (0, ScaledDividend);
ScaleUpByDigit (ScaledDividend, ScalingFactor);
ScaledDivisor := CopyBidirectionalList (Divisor);
ScaleUpByDigit (ScaledDivisor, ScalingFactor);

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

FirstTwo := BaseOfNumeration *
FirstOfBidirectionalList (ScaledDivisor) +
RecoverByPositionFromBidirectionalList (2, ScaledDivisor);

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

Residue := MakeZero;
CursorToStartOfBidirectionalList (ScaledDividend);
for Position := 1 to LengthOfBidirectionalList (ScaledDivisor) + 1
do begin
AppendToBidirectionalList (Residue,
ElementAtCursorInBidirectionalList (ScaledDividend));
end;

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

Quotient := MakeZero;
Finished := False;
while not Finished do begin

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

TrialDigit := FindTrialDigit (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. }

TrialProduct := MultiplyNaturalByDigit (ScaledDivisor, TrialDigit);
{ so that the LessNatural function will work }
if LessNatural (Residue, TrialProduct) then begin
TrialDigit := TrialDigit - 1;
DeallocateNatural (TrialProduct);
TrialProduct := MultiplyNaturalByDigit (ScaledDivisor, TrialDigit);
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. }

AppendToBidirectionalList (Quotient, TrialDigit);
DecrementNatural (Residue, TrialProduct);
DeallocateNatural (TrialProduct);
while LengthOfBidirectionalList (Residue) <
LengthOfBidirectionalList (ScaledDivisor) do
PrependToBidirectionalList (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 NullCursorInBidirectionalList (ScaledDividend) then
Finished := True
else begin
AppendToBidirectionalList (Residue,
ElementAtCursorInBidirectionalList (ScaledDividend));
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. }

DeallocateNatural (ScaledDividend);
DeallocateNatural (ScaledDivisor);
DeallocateNatural (Residue)

end

end;

function QuotientNatural (Dividend, Divisor: Natural): Natural;
var
Quotient: Natural;
{ the quotient resulting from the division }
{ the remainder -- to be discarded }
begin
QuotientNatural := Quotient
end;

function RemainderNatural (Dividend, Divisor: Natural): Natural;
var
Remainder: Natural;
{ the remainder resulting from the division }
{ the quotient -- to be discarded }
begin
RemainderNatural := Remainder
end;

function RaiseNatural (Base, Exponent: Natural): Natural;
var
CopyOfBase, CopyOfExponent: 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 }
Temporary: Natural;
{ temporary storage for the result of a computation that will
ultimately be transferred to one of the principal variables }
{ the last bit of the current exponent, discarded once it has been
examined }
begin
CopyOfBase := CopyBidirectionalList (Base);
CopyOfExponent := CopyBidirectionalList (Exponent);
Accumulator := PascalIntegerToNatural (1);
while not ZeroNatural (CopyOfExponent) do begin
if OddNatural (CopyOfExponent) then begin
Temporary := MultiplyNatural (Accumulator, CopyOfBase);
DeallocateNatural (Accumulator);
Accumulator := Temporary
end;
Temporary := SquareNatural (CopyOfBase);
DeallocateNatural (CopyOfBase);
CopyOfBase := Temporary;
DeallocateNatural (CopyOfExponent);
CopyOfExponent := Temporary
end;
DeallocateNatural (CopyOfBase);
DeallocateNatural (CopyOfExponent);
RaiseNatural := Accumulator
end;

function SuccessorOfNatural (Operand: 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 }
Digit: BaseDigits;
{ one digit at a time from the operand }
begin
Carry := 1;
Result := MakeZero;
CursorToEndOfBidirectionalList (Operand);
while not NullCursorInBidirectionalList (Operand) do begin
Digit := ElementAtCursorInBidirectionalList (Operand);
if Carry = 0 then
PrependToBidirectionalList (Digit, Result)
else if Digit = BaseOfNumeration - 1 then
PrependToBidirectionalList (0, Result)
else begin
PrependToBidirectionalList (Digit + 1, Result);
Carry := 0
end;
RetreatCursorAlongBidirectionalList (Operand)
end;
if Carry = 1 then
PrependToBidirectionalList (1, Result);
SuccessorOfNatural := Result
end;

function PredecessorOfNatural (Operand: Natural): Natural;
var
Borrow: Bit;
{ initially, the 1 that must be subtracted; subsequently, a 1 so long
as a borrow must be propagated, and 0 thereafter }
Result: Natural;
{ the predecessor natural number as it is constructed, digit by
digit }
Digit: BaseDigits;
{ one digit at a time from the operand }
begin
Assert (not ZeroNatural (Operand), PredecessorException,
NaturalExceptionHandler);
Borrow := 1;
Result := MakeZero;
CursorToEndOfBidirectionalList (Operand);
while not NullCursorInBidirectionalList (Operand) do begin
Digit := ElementAtCursorInBidirectionalList (Operand);
if Borrow = 0 then
PrependToBidirectionalList (Digit, Result)
else if Digit = 0 then
PrependToBidirectionalList (BaseOfNumeration - 1, Result)
else begin
PrependToBidirectionalList (Digit - 1, Result);
Borrow := 0
end;
RetreatCursorAlongBidirectionalList (Operand)
end;
PredecessorOfNatural := Result
end;

function TwiceNatural (Operand: 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 }
PlaceProduct: Integer;
{ the double of one digit of the given natural number, plus any Carry
from the next lower digit position }
begin
Result := MakeZero;
Carry := 0;
CursorToEndOfBidirectionalList (Operand);
while not NullCursorInBidirectionalList (Operand) do begin
PlaceProduct := ElementAtCursorInBidirectionalList (Operand) * 2 + Carry;
if PlaceProduct < BaseOfNumeration then begin
PrependToBidirectionalList (PlaceProduct, Result);
Carry := 0
end
else begin
PrependToBidirectionalList (PlaceProduct - BaseOfNumeration, Result);
Carry := 1
end;
RetreatCursorAlongBidirectionalList (Operand)
end;
if Carry = 1 then
PrependToBidirectionalList (1, Result);
TwiceNatural := Result
end;

function SquareNatural (Operand: Natural): Natural;
begin
SquareNatural := MultiplyNatural (Operand, Operand)
end;

function CubeNatural (Operand: Natural): Natural;
var
Sq: Natural;
{ the square of the given number }
begin
Sq := MultiplyNatural (Operand, Operand);
CubeNatural := MultiplyNatural (Sq, Operand);
DeallocateNatural (Sq)
end;

function EqualNaturals (LeftOperand, RightOperand: Natural): Boolean;
var
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
begin
if LeftOperand = RightOperand then
EqualNaturals := True
else if LengthOfBidirectionalList (LeftOperand) <>
LengthOfBidirectionalList (RightOperand) then
EqualNaturals := False
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := False;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
EqualNaturals := True;
Finished := True
end
else if ElementAtCursorInBidirectionalList (LeftOperand) =
ElementAtCursorInBidirectionalList (RightOperand) then begin
end
else begin
EqualNaturals := False;
Finished := True
end
until Finished
end
end;

function UnequalNaturals (LeftOperand, RightOperand: Natural): Boolean;
var
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
begin
if LeftOperand = RightOperand then
UnequalNaturals := False
else if LengthOfBidirectionalList (LeftOperand) <>
LengthOfBidirectionalList (RightOperand) then
UnequalNaturals := True
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := False;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
UnequalNaturals := False;
Finished := True
end
else if ElementAtCursorInBidirectionalList (LeftOperand) =
ElementAtCursorInBidirectionalList (RightOperand) then begin
end
else begin
UnequalNaturals := True;
Finished := True
end
until Finished
end
end;

function LessNatural (LeftOperand, RightOperand: Natural): Boolean;
var
LeftSize, RightSize: Integer;
{ the number of digits in the left and right operands, respectively }
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
LeftDigit, RightDigit: BaseDigits;
{ one digit at a time from the left and right operands, respectively }
begin
if LeftOperand = RightOperand then
LessNatural := False
else begin
LeftSize := LengthOfBidirectionalList (LeftOperand);
RightSize := LengthOfBidirectionalList (RightOperand);
if LeftSize < RightSize then
LessNatural := True
else if RightSize < LeftSize then
LessNatural := False
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := FALSE;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
LessNatural := False;
Finished := True
end
else begin
LeftDigit := ElementAtCursorInBidirectionalList (LeftOperand);
RightDigit := ElementAtCursorInBidirectionalList (RightOperand);
if LeftDigit < RightDigit then begin
LessNatural := True;
Finished := True
end
else if RightDigit < LeftDigit then begin
LessNatural := False;
Finished := True
end
else begin
end
end
until Finished
end
end
end;

function LessOrEqualNatural (LeftOperand, RightOperand: Natural): Boolean;
var
LeftSize, RightSize: Integer;
{ the number of digits in the left and right operands, respectively }
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
LeftDigit, RightDigit: BaseDigits;
{ one digit at a time from the left and right operands, respectively }
begin
if LeftOperand = RightOperand then
LessOrEqualNatural := True
else begin
LeftSize := LengthOfBidirectionalList (LeftOperand);
RightSize := LengthOfBidirectionalList (RightOperand);
if LeftSize < RightSize then
LessOrEqualNatural := True
else if RightSize < LeftSize then
LessOrEqualNatural := False
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := FALSE;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
LessOrEqualNatural := True;
Finished := True
end
else begin
LeftDigit := ElementAtCursorInBidirectionalList (LeftOperand);
RightDigit := ElementAtCursorInBidirectionalList (RightOperand);
if LeftDigit < RightDigit then begin
LessOrEqualNatural := True;
Finished := True
end
else if RightDigit < LeftDigit then begin
LessOrEqualNatural := False;
Finished := True
end
else begin
end
end
until Finished
end
end
end;

function GreaterNatural (LeftOperand, RightOperand: Natural): Boolean;
var
LeftSize, RightSize: Integer;
{ the number of digits in the left and right operands, respectively }
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
LeftDigit, RightDigit: BaseDigits;
{ one digit at a time from the left and right operands, respectively }
begin
if LeftOperand = RightOperand then
GreaterNatural := False
else begin
LeftSize := LengthOfBidirectionalList (LeftOperand);
RightSize := LengthOfBidirectionalList (RightOperand);
if LeftSize < RightSize then
GreaterNatural := False
else if RightSize < LeftSize then
GreaterNatural := True
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := FALSE;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
GreaterNatural := False;
Finished := True
end
else begin
LeftDigit := ElementAtCursorInBidirectionalList (LeftOperand);
RightDigit := ElementAtCursorInBidirectionalList (RightOperand);
if LeftDigit < RightDigit then begin
GreaterNatural := False;
Finished := True
end
else if RightDigit < LeftDigit then begin
GreaterNatural := True;
Finished := True
end
else begin
end
end
until Finished
end
end
end;

function GreaterOrEqualNatural (LeftOperand, RightOperand: Natural): Boolean;
var
LeftSize, RightSize: Integer;
{ the number of digits in the left and right operands, respectively }
Finished: Boolean;
{ indicates whether the value to be returned has been conclusively
established yet }
LeftDigit, RightDigit: BaseDigits;
{ one digit at a time from the left and right operands, respectively }
begin
if LeftOperand = RightOperand then
GreaterOrEqualNatural := True
else begin
LeftSize := LengthOfBidirectionalList (LeftOperand);
RightSize := LengthOfBidirectionalList (RightOperand);
if LeftSize < RightSize then
GreaterOrEqualNatural := False
else if RightSize < LeftSize then
GreaterOrEqualNatural := True
else begin
CursorToStartOfBidirectionalList (LeftOperand);
CursorToStartOfBidirectionalList (RightOperand);
Finished := FALSE;
repeat
if NullCursorInBidirectionalList (LeftOperand) then begin
GreaterOrEqualNatural := True;
Finished := True
end
else begin
LeftDigit := ElementAtCursorInBidirectionalList (LeftOperand);
RightDigit := ElementAtCursorInBidirectionalList (RightOperand);
if LeftDigit < RightDigit then begin
GreaterOrEqualNatural := False;
Finished := True
end
else if RightDigit < LeftDigit then begin
GreaterOrEqualNatural := True;
Finished := True
end
else begin
end
end
until Finished
end
end
end;

function MajorNatural (LeftOperand, RightOperand: Natural): Natural;
var
Result: Natural;
begin
if GreaterOrEqualNatural (LeftOperand, RightOperand) then
MajorNatural := LeftOperand
else
MajorNatural := RightOperand
end;

function MinorNatural (LeftOperand, RightOperand: Natural): Natural;
var
Result: Natural;
begin
if LessOrEqualNatural (LeftOperand, RightOperand) then
MinorNatural := LeftOperand
else
MinorNatural := RightOperand
end;

function ZeroNatural (Operand: Natural): Boolean;
begin
ZeroNatural := EmptyBidirectionalList (Operand)
end;

function MultipleNatural (Candidate, Unit: Natural): Boolean;
var
Remainder: Natural;
{ the remainder when Candidate is divided by Unit }
begin
if ZeroNatural (Unit) then
MultipleNatural := ZeroNatural (Candidate)
else begin
Remainder := RemainderNatural (Candidate, Unit);
MultipleNatural := ZeroNatural (Remainder);
DeallocateNatural (Remainder)
end
end;

function EvenNatural (Operand: Natural): Boolean;
begin
if EmptyBidirectionalList (Operand) then
EvenNatural := True
else
EvenNatural := not Odd (LastOfBidirectionalList (Operand))
end;

function OddNatural (Operand: Natural): Boolean;
begin
if EmptyBidirectionalList (Operand) then
OddNatural := False
else
OddNatural := Odd (LastOfBidirectionalList (Operand))
end;

procedure ReadNatural (var Source: Text; var Legend: Natural;
var Success: Boolean);
const
DigitZero = 48 { = ord ('0') };
var
Finished: Boolean;
{ indicates whether all of the digits of the numeral have been read
in }

{ The SkipWhiteSpace procedure advances through a text file until
a non-whitespace character (or the end of the file) is encountered. }

procedure SkipWhiteSpace (var Source: Text);
const
Space = ' ';
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^ <= Space) or (Source^ = Chr (127)) then
Get (Source)
else
Finished := True
end;

procedure IncrementByDigit (var Given: Natural; Amount: BaseDigits);
var
PlaceSum: Integer;
{ the sum of the last digit of the given natural number and the
amount by which it is to be incremented }
Carry: Bit;
{ a 1 so long as a carry must be propagated, and 0 thereafter }
Digit: BaseDigits;
{ one digit at a time from Given }
begin
if Amount <> 0 then begin
CursorToEndOfBidirectionalList (Given);
if NullCursorInBidirectionalList (Given) then
PrependToBidirectionalList (Amount, Given)
else begin
PlaceSum := ElementAtCursorInBidirectionalList (Given) + Amount;
if PlaceSum < BaseOfNumeration then
AssignAtCursorInBidirectionalList (Given, PlaceSum)
else begin
AssignAtCursorInBidirectionalList (Given,
PlaceSum - BaseOfNumeration);
Carry := 1;
RetreatCursorAlongBidirectionalList (Given);
while (Carry = 1) and not NullCursorInBidirectionalList (Given)
do begin
Digit := ElementAtCursorInBidirectionalList (Given);
if Digit = BaseOfNumeration - 1 then begin
AssignAtCursorInBidirectionalList (Given, 0);
RetreatCursorAlongBidirectionalList (Given)
end
else begin
AssignAtCursorInBidirectionalList (Given, Digit + 1);
Carry := 0
end
end;
if Carry = 1 then
PrependToBidirectionalList (1, Given)
end
end
end
end;

function Numeric (Ch: Char): Boolean;
begin
Numeric := ('0' <= Ch) and (Ch <= '9')
end;

SkipWhiteSpace (Source);
if EOF (Source) then
Success := False
else if not Numeric (Source^) then
Success := False
else begin
Legend := MakeZero;
Finished := False;
repeat
ScaleUpByDigit (Legend, 10);
IncrementByDigit (Legend, Ord (Source^) - DigitZero);
Get (Source);
if EOF (Source) then
Finished := True
else
Finished := not Numeric (Source^)
until Finished;
Success := True
end
end;

procedure WriteNatural (var Target: Text; Scribend: 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 Scribend; later, the part of the number not
yet reduced to decimal digits }
Quotient: Natural;
{ the quotient resulting from a division of Working by 10 }
Remainder: BaseDigits;
{ the remainder Resulting from such a division }
begin
if ZeroNatural (Scribend) then
Write (Target, '0')
else begin
Digits := CreateStack;
Working := CopyBidirectionalList (Scribend);
while not ZeroNatural (Working) do begin
DivideNaturalByDigit (Working, 10, Quotient, Remainder);
DeallocateNatural (Working);
Working := Quotient;
PushToStack (Remainder, Digits)
end;
DeallocateNatural (Working);
while not EmptyStack (Digits) do
Write (Target, PopFromStack (Digits) : 1);
DeallocateStack (Digits)
end
end;

procedure AssignNatural (var Target: Natural; Source: Natural);
begin
Target := CopyBidirectionalList (Source)
end;

function PascalIntegerToNatural (N: Integer): Natural;
var
Result: Natural;
{ the natural number of value equal to n, constructed digit by digit
through successive divisions of n }
begin
Assert (0 <= N, NegativeException, NaturalExceptionHandler);
Result := MakeZero;
while N <> 0 do begin
PrependToBidirectionalList (N mod BaseOfNumeration, Result);
N := N div BaseOfNumeration
end;
PascalIntegerToNatural := Result
end;

function NaturalToPascalInteger (N: Natural): Integer;
var
Result: Integer;
{ the integer of value equal to N, constructed by evaluating n
digit by digit }
Digit: BaseDigits;
{ one digit at a time from N }

{ The InBounds function determines, cautiously, whether the integer
that would be obtained by multiplying value and base and adding
NewDigitValue the Result is within HP Pascal's integer data type
-- that is, whether it would be less than or equal to MaxInt --
returning True if the proposed operation is safe and False if it
would result in an integer overflow. }

function InBounds (Value: Integer; NewDigitValue: Integer): Boolean;
const
AllButLast = MaxInt div BaseOfNumeration;
Last = MaxInt mod BaseOfNumeration;
begin
if Value < AllButLast then
InBounds := True
else if Value = AllButLast then
InBounds := (NewDigitValue <= Last)
else
InBounds := False
end;

begin { function NaturalToPascalInteger }
Result := 0;
CursorToStartOfBidirectionalList (N);
while not NullCursorInBidirectionalList (N) do begin
Digit := ElementAtCursorInBidirectionalList (N);
Assert (InBounds (Result, Digit), IntegerOverflowException,
NaturalExceptionHandler);
Result := Result * BaseOfNumeration + Digit;
end;
NaturalToPascalInteger := Result
end;

procedure DebugOutputNatural (var Target: Text; Scribend: Natural);
begin
CursorToStartOfBidirectionalList (Scribend);
while not NullCursorInBidirectionalList (Scribend) do begin
Write (Target, '| ',
ElementAtCursorInBidirectionalList (Scribend): DigitWidth,
' ');
end;
Write (Target, '|')
end;

procedure DeallocateNatural (var N: Natural);
begin
DeallocateBidirectionalList (N)
end;

end.
```

### Integers

A superior version of the `Integer` data type can now be constructed by representing an integer as a natural number (the magnitude) and a sign:

```record
Sign: (Positive, Negative);
Magnitude: Natural
end
```
It is straightforward to implement all of the operations described in the handout on integers for this data type in terms of the operations on natural numbers; writing up the module is left as an exercise for the reader.

This document is available on the World Wide Web as

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

created January 24, 1996
last revised November 8, 1996