/* arith.c:  Number-Theoretic arithmetic routines
** Written by Michael J. Gourlay
** Fall 1996, University of Colorado at Boulder
** Algorithms
*/

/*

TOP LEVEL DESCRIPTION:
----------------------

The code in this module implements some modular arithmetic routines
for arbitrarily large numbers.  The "large numbers" are described by
an abstract data type, which is called a NumberT, which contains an
array of "units", the length of the big number (where the length is
the number of "units", where a "unit" is analogous to a "digit"
except that it takes several characters to print out the value of a
"unit"), and the "base" of the number (where each "unit" can take a
value between 0 and the base value, including zero, but not including
the base value).

The values of "base" and "length" can take on a large range of values.
In particular, the value of "base" can be any number representable by
the "unit" primative type.  The "base" value is set at run time.
However, see the note below about the arithmetic operations and their
assumptions about the "base".

The "unit" primative type is easily modified to hold numbers of
different sizes, but this modification must be done at compile time.

The modular arithmetic routines implemented are addition, subtraction,
and multiplication.  These routines are implemented using more
primative arithmetic routines which perform non-modular arithmetic.

In the case of modular subtraction, and modular addition, the input
values are assumed to be modular already.  This speeds up those
routines, because making their intermediate results modular is easier
if the input values are modular.  In the case of subtraction, the
result of the non-modular subtraction is equal to the modular value,
so no modulus finding is done at the end.  In the case of addition,
the result of the non-modular addition can be at most larger than the
modulus by the value of the modulus, so finding the modulus of an
addition operation is at worst obtainable by a single subtraction, as
long as the input values are modular.

In the case of multiplication, no reasonable conditions on the input
values guarentee that the result will be anywhere close to the modular
value.  Therefore, a long division routine, which finds the quotient
and the remainder, was implemented.  The quotient value is never used
in modular arithmetic, so that value is discarded.  The remainder is
the modular value.

It is conceivable that different big numbers could have different
"base" values, and the data structure allows for that.  However, in
the interest of speed and simplicity of coding, the arithmetic
routines assume that the input values have the same "base".  No check
is done to verify this, however, so caveat emptor.

There are also auxiliary routines that create, destroy, read, print,
copy, and zero large numbers.

In order to pre-condition the input values (for the benefit of the
addition and subtraction routines), the values read in from the input
file are modularized by the "read" routine.

The "print" routine assumes that the "units" are most easily
displayed as decimal values (as opposed to, e.g., hexadecimal or
octal).

In addition to the above routines, this module contains an example
program that makes use of all of the features of this abstract data
type.  The program reads in a file containing three big numbers:  a
modulus, and two operands.  (The program defines the example numbers
to be no larger than 20 units, but this is not a restriction of the
data type itself -- it is merely a choice made for testing purposes.)
Then, the program performs each modular arithmetic operation on the
operands, in both directions, i.e.  n1(op)n2 and n2(op)n1.  In the
cases of addition and multiplication, the two results should be the
same, since addition and multiplication are commutative operators.
However, subtraction is not commutative, so the two subtraction
results should be different.

*/




#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>




/* VERBOSE: precompiler token for verbose reporting */
/* #define VERBOSE /**/

/* PARANOID: Check assumptions.  Slows code some. */
#define PARANOID




/* NumberUnitT is the primative type that stores a "unit" of a NumberT
*/
typedef unsigned int NumberUnitT;


/* BITS_PER_UNIT: number of bits used per "unit" of a NumberT
** Assumes bytes returned by sizeof() are 8 bits.
*/
#define BITS_PER_UNIT ((4 * sizeof(NumberUnitT)))


/* NUMBER_BASE_MAX is the maximum base value of a "unit" in a NumberT.
*/
#define NUMBER_BASE_MAX (1<<BITS_PER_UNIT)




/* NumberT stores a "big natural number":
**
** NumberT has an array of some primative integer type.  Each element in
** that array contains a single "unit" for the number.
**
** units[0] is the least significant unit.
** units[len-1] is the most significant unit.
*/
typedef struct {
  NumberUnitT* units;  /* array of units */
  int len;             /* length (number of units) */
  int base;            /* base of a unit, where  0 <= unit value < base */
} NumberT;




NumberT *debug = NULL;




/* ----------------------- */
/*    Utility functions    */
/* ----------------------- */


#define MIN(x,y) ((x)<(y)?(x):(y))


/* my_calloc: memory allocation with error checking
** Arguments are the same as for standard library function "calloc"
** Exit program if memory allocation fails.
*/
void *
my_calloc(size_t nelem, size_t elsize)
{
  char *mem;
  if((mem=(char *)calloc(nelem, elsize))==NULL) {
    fprintf(stderr, "my_calloc: out of memory.  Tried %i x %i\n",
            nelem, elsize);
    abort();
    exit(1);
  }
  return mem;
}


/* MY_CALLOC: macro to do memory allocation easily
** ne (in): number of elements.  Must be integer.
** type (in): type token (not a value).  Must be a valid variable type.
*/
#define MY_CALLOC(ne,type) ((type *)my_calloc(ne,sizeof(type)))




/* -------------------- */
/*    Number methods    */
/* -------------------- */

/* numAlloc: allocate and initialize to zero a NumberT
** num_len (in): length (number of units) to make the NumberT
** num_base (in): base of the Number's Units
*/
void
numAlloc(NumberT *num, int num_len, int num_base)
{
  num->len = num_len;
  num->base = num_base;
  num->units = MY_CALLOC(num_len, NumberUnitT);
}




/* numFree: Free memory for NumberT */
void
numFree(NumberT *num)
{
  free(num->units);
  num->units = NULL;
  num->len = 0;
  num->base = 0;
}




/* numZero: set a NumberT to zero
*/
void
numZero(NumberT *num)
{
  memset(num->units, 0, sizeof(NumberUnitT) * num->len);
}




/* numCopy: Copy units of in into out
** Truncation is possible if in->len != out->len
*/
void
numCopy(NumberT *out, const NumberT *in)
{
  memcpy(out->units, in->units, sizeof(NumberUnitT) * MIN(out->len,in->len));
}




/* numPrintContiguous: print a NumberT withOUT spaces between units
*/
void
numPrintContiguous(const NumberT *num)
{
  int ni;

  if(num == NULL) return;

  for(ni = num->len-1; ni >= 0; ni--)
  {
    printf("%04i", num->units[ni]);
  }
  printf("\n");
}




/* numPrint: print a NumberT with spaces between units
*/
void
numPrint(const NumberT *num)
{
  int ni;

  if(num == NULL) return;

  for(ni = num->len-1; ni >= 0; ni--)
  {
    printf("%04i ", num->units[ni]);
  }
  printf("\n");
}




/* forward declaration */
void numLongDiv(const NumberT *num,const NumberT *den,NumberT *quo,NumberT *rem);




/* numRead: read a NumberT from a file and take its modulus
**
** If mod == NULL then do not find modulus.
** (This allows for reading in the modulus)
**
** File must already be open for reading.
** File position is advanced by one line during this routine.
** Therefore, there can be only one NumberT per line.
** Also, a NumberT can occupy no more than one line.
**
** number in file may have up to twice as many units as num->len.
**
** The format of the input NumberT must have units separated by white space.
** Also, units must not begin with "0".
**
** Assumes that the input num is large enough to contain the modulus.
*/
void
numRead(NumberT *num, FILE *fileP, const NumberT *mod)
{
  NumberT num_in;   /* temp space for input NumberT */
  char line[4096];  /* input line buffer */
  int ii;           /* index of unit being read */
  int si;           /* pointer into buffer line being scanned */
  int nc;           /* number of chars read in during a scanf */

  /* Set input number to zero */
  numZero(num);

  /* Allocate a larger-than-needed NumberT */
  numAlloc(&num_in, num->len * 2, num->base);

  /* Read in a single line from the file */
  fgets(line, 4095, fileP);

  /* Read the units of that line */
  for(ii=0, si=0;
      (ii < (num_in.len))
        && (sscanf(&line[si], "%i%n", &num_in.units[ii], &nc) >0) ;
      ii++, si += nc)
  {
    if(feof(fileP)) {
      printf("numRead: end of file\n");
      exit(2);
    }
    if(num_in.units[ii] >= num->base) {
      printf("numRead: unit %i bigger than base %i\n",
        num_in.units[ii], num->base);
    }
#ifdef DEBUG
    printf("at %i val[%i]=%i nc=%i\n", si, ii, num_in.units[ii], nc);
#endif
  }

  /* Make ii be the value of the index of the input number's last unit */
  ii --;

  /* Now units in num_in are out of order to re-order them in place */
  for(si=0; si <= ii/2; si ++) {
    /* nc is used here as temp space */
    nc = num_in.units[si];
    num_in.units[si] = num_in.units[ii-si];
    num_in.units[ii-si] = nc;
  }

  /* Now (maybe) find the modulus of the input NumberT */
  if(mod != NULL) {
    /* Find modulus */
    numLongDiv(&num_in, mod, NULL, num);
  } else {
    /* Don't find modulus  -- just use raw input NumberT */
    numCopy(num, &num_in);
  }

  numFree(&num_in);
}




/* numLeadNonZero: Find most significant non-zero unit
** Return index of that unit.
** Returns -1 if NumberT is zero.
*/
int
numLeadNonZero(const NumberT *num)
{
  int ni;

  for(ni = num->len - 1; (ni >= 0) && (!num->units[ni]) ; ni --)
    ;

  return ni;
}




/* numCompare: compare two NumberT's
**
** if n1 > n2, return positive
** if n1 == n2, return 0
** if n1 < n2, return negative
**
** Assumes n1->base == n2->base
*/
int
numCompare(const NumberT *n1, const NumberT *n2)
{
  int ni;
  int lnz1, lnz2;

  /* If lengths are different, then values are different */
  if((lnz1=numLeadNonZero(n1)) > (lnz2=numLeadNonZero(n2))) {
    /* n1 > n2 */
    return 1;
  } else if(lnz1 < lnz2) {
    /* n1 < n2 */
    return -1;
  } else {
    /* Lengths are the same -- do more detailed comparison */
    for(ni = lnz1; ni >=0 ; ni--) {
      if(n1->units[ni] != n2->units[ni])
        return n1->units[ni] - n2->units[ni];
    }
  }
  return 0;
}




/* numAddPrim: Add (primative) two NumberT's , no MODULUS
**
** Return sum = n1 + n2
** where most significant unit returned in sum can possibly be larger
** than strictly allowed.
**
** Assumes n1->len == n2->len
** Assumes sum->len >= n1->len
** Assumes n1->base == n2->base
**
** This allows the result of an addition to overflow slightly but
** still be handled as modulus arithmetic even in the case where the
** modulus is the largest representable integer.
**
** Addition is performed from least significant unit to most
** significant unit, with carry.
*/
void
numAddPrim(const NumberT *n1, const NumberT *n2, NumberT *sum)
{
  int ni;
  unsigned long carry;
  unsigned long tmp;

#ifdef VERBOSE
  printf("addend 1 = "); numPrintContiguous(n1);
  printf("addend 2 = "); numPrintContiguous(n2);
#endif /* VERBOSE */

#ifdef PARANOID
  if(n1->len != n2->len) {
    printf("numAddPrim: n1->len=%i != n2->len=%i\n", n1->len, n2->len);
    abort();
  }
  if(sum->len < n1->len) {
    printf("numAddPrim: sum->len=%i < n1->len=%i\n", sum->len, n1->len);
    abort();
  }
#endif /* PARANOID */

  sum->base = n1->base;

  for(ni=0, carry=0; ni < n1->len; ni++) {
    tmp = n1->units[ni] + n2->units[ni] + carry;
    carry = 0;

    /* Unless we're at the most significant unit, find the "carry" amount */
    if((ni<sum->len-1) && (tmp >= n1->base)) {
      carry = tmp / n1->base;
      tmp %= n1->base;
    }
    sum->units[ni] = tmp;
  }
}




/* numSubPrim: Subtract (primative) two NumberT's
**
** return dif = n1 - n2  for n1 > n2
**
** i.e. n1 is the minuend, n2 is the subtrahend
**
** Assumes n1->len == n2->len <= dif->len
** Assumes n1->base == n2->base
**
** minuend must be larger than the subtrahend,
** otherwise result is undefined.
**
** Subtraction is performed from least significant unit to most
** significant unit, using complementary subtraction with borrow.
*/
void
numSubPrim(const NumberT *n1, const NumberT *n2, NumberT *dif)
{
  int ni;
  int borrow;

#ifdef VERBOSE
  printf("minuend    = "); numPrintContiguous(n1);
  printf("subtrahend = "); numPrintContiguous(n2);
#endif /* VERBOSE */

#ifdef PARANOID
  if(n1->len > n2->len) {
    printf("numSubPrim: n1->len=%i > n2->len=%i\n", n1->len, n2->len);
    abort();
  }
  if(numCompare(n1,n2) < 0) {
    printf("numSubPrim: n1 < n2\n");
    printf("n1 = "); numPrint(n1);
    printf("n2 = "); numPrint(n2);
    abort();
  }
#endif

  dif->base = n1->base;

  for(ni=0, borrow=0;  ni < n1->len;  ni++) {
    if(((signed)n1->units[ni] - borrow) < (signed)n2->units[ni]) {
      /* Minuend unit is smaller than subtrahend unit, so
      ** find difference using complementary subtraction.
      */
      dif->units[ni] = n1->units[ni] - borrow + (n1->base - n2->units[ni]);

      /* Remember the borrow from the next minuend unit */
      borrow = 1;
    } else {
      /* Minuend unit is larger than subtrahend unit so
      ** find difference using ordinary subtraction.
      */
      dif->units[ni] = n1->units[ni] - borrow - n2->units[ni];
      borrow = 0;
    }
  }
}




/* numSubtract: Subtract two NumberT's , MODULUS mod
**
** return dif = n1 - n2 MODULUS mod
**
** i.e. n1 is the minuend, n2 is the subtrahend
**
** Input numbers are assumed modular.  No modularizing is done here.
** It's faster this way.
**
** If the minuend is smaller than the subtrahend, then the result would
** be negative in non-modular arithmetic.  However, since this routine
** performs modulus subtraction, we can avoid needing to represent a
** negative result by adding the modulus to the minuend before
** subtraction.  This requires an addition without modulus.
*/
void
numSubtract(const NumberT *n1,const NumberT *n2,NumberT *dif,const NumberT *mod)
{
  NumberT minuend;
  int allocd = 0;

  if(numCompare(n1,n2) < 0) {
    /* Difference would be negative so add modulus to minuend */
    numAlloc(&minuend, n1->len, n1->base);
    numAddPrim(n1, mod, &minuend);
    allocd = 1;
  } else {
    minuend = *n1;
  }

  /* Perform the non-modulus subtraction */
  numSubPrim(&minuend, n2, dif);

  /* dif should now be the modulus answer.  If it isn't, then the
  ** input numbers weren't already modular, in which case the result
  ** is undefined.
  */

  if(allocd) numFree(&minuend);
}




/* numAdd: Add two NumberT's , MODULUS mod
**
** return sum = n1 + n2 MODULUS mod
**
** Input numbers are assumed modular.  Not more than one subtraction
** should be done at the end to modularize the result, and that will
** only work if the input is modular.  I could just find the remainder
** of a division, but that takes longer, and it's no more imposing to
** just have the inputs be modular.
*/
void
numAdd(const NumberT *n1, const NumberT *n2, NumberT *sum, const NumberT *mod)
{
  /* Add, but allow for sum to be "too big" */
  numAddPrim(n1, n2, sum);

  /* Now find modulus of sum. */
  if(numCompare(sum, mod) >= 0) {
    numSubPrim(sum, mod, sum);
  }
}




/* numMultPrimUnit: Multiply (primative) a NumberT by a single "unit"
**
** unit (in): single "unit" multiplicand
** big (in): big number multiplicand
** prod (out): product = unit * big.  Assumes prod->len >= big->len + shift + 1
** shift (in): amount to shift left the product
*/
void
numMultPrimUnit(NumberUnitT unit, const NumberT *big,
                NumberT *prod, int shift)
{
  unsigned long tmp;
  unsigned long carry;
  int ni;

#ifdef PARANOID
  if((big->len + shift + 1)  >  prod->len) {
    printf("numMultPrimUnit: big->len=%i,+stuff > prod->len=%i\n",
           big->len, prod->len);
    abort();
  }
#endif /* PARANOID */

  prod->base = big->base;

  /* Compute product of a Number with a unit */
  for(ni=0, carry=0; ni < big->len; ni++) {
    tmp = (unit * big->units[ni]) + carry;
    carry = 0;
    /* Find the "carry" amount */
    if(tmp >= big->base) {
      carry = tmp / big->base;
      tmp %= big->base;
    }
    prod->units[ni + shift] = tmp;
  }
  if(carry) {
    prod->units[ni + shift] = carry;
  }
}




/* numMultPrim: Multiply (primative) two NumberT's
**
** return prod = n1 * n2
**
** prod(out): product.  Assumes prod->len >= n1->len + n2->len
**
** Assumes n1->len == n2->len.
** Assumes n1->base == n2->base.
** Assumes prod->len >= 2 * n1->len
**
** Product can have roughly as many digits as the sum of the number
** of digits of the multiplicands.
*/
void
numMultPrim(const NumberT *n1, const NumberT *n2, NumberT *prod)
{
  NumberT prod_tmp;
  int pi;

#ifdef VERBOSE
  printf("multiplicand 1 = "); numPrintContiguous(n1);
  printf("multiplicand 2 = "); numPrintContiguous(n2);
#endif /* VERBOSE */
#ifdef PARANOID
  if(n1->len != n2->len) {
    printf("numMultPrim: n1->len=%i != n2->len=%i\n", n1->len, n2->len);
  }
  if(prod->len < n1->len + n2->len) {
    printf("numMultPrim:prod->len=%i too short\n", prod->len);
  }
#endif /* PARANOID */

  /* Create temporary calculation workspace */
  numAlloc(&prod_tmp, prod->len, n1->base);

  /* Set product holder to zero */
  prod->base = n1->base;

  /* Loop through each unit in the bottom multiplicand */
  for(pi=0; pi < n1->len; pi++) {
    /* Only multiply if multiplicand is non-zero */
    if(n2->units[pi]) {
      /* Set prod_tmp to zero */
      numZero(&prod_tmp);

      /* Compute product of top Number with a bottom unit */
      numMultPrimUnit(n2->units[pi], n1, &prod_tmp, pi);

      /* Add this new product to the cumulative sum */
      numAddPrim(prod, &prod_tmp, prod);
    }
  }
  numFree(&prod_tmp);
}




#include "arith-update.c"




/* numLongDiv: Divide two NumberT's using long division
**
** return quo = num / den  and  rem = num MOD den
**
** num (in): numerator or dividend, must be twice the size of den
** den (in): denominator or divisor
** quo (out): quotient. not stored if quo == NULL
** rem (out): remainder. not stored if rem == NULL
**
** Assumed num->len >= 2 * den->len  (is this realy true?)
** den->len == quo->len == rem->len
**
** Can NOT do in-place divisions; i.e. quo must not be num or dem.
*/
void
numLongDiv(const NumberT *num, const NumberT *den, NumberT *quo, NumberT *rem)
{
  int lnz_num;       /* index of leading non-zero unit of numerator */
  int lnz_den;       /* index of leading non-zero unit of denominator */
  NumberT un;        /* temp space numerator */
  NumberT uns;       /* subset of temp space numerator */
  NumberT tmp;       /* temp space */
  int ji;            /* loop variable */
  unsigned long qh;  /* trial quotient */

#ifdef VERBOSE
  printf("numerator   = "); numPrintContiguous(num);
  printf("denominator = "); numPrintContiguous(den);
#endif /* VERBOSE */

#ifdef DEBUG_LD
  printf("numerator   = "); numPrint(num);
  printf("denominator = "); numPrint(den);
#endif

  /* Empty the quotient and remainder holders */
  if(quo) quo->base = num->base;
  if(quo) numZero(quo);
  if(rem) rem->base = num->base;
  if(rem) numZero(rem);

  /* Check for divide by zero  */
  if(!numCompareZero(den)) {
    printf("numLongDiv: divide by zero attempted\n");
    return;
  }

  /* Find leading non-zero units of numerator and denominator */
  lnz_num = numLeadNonZero(num);
  lnz_den = numLeadNonZero(den);

  if(lnz_num < lnz_den) {
    /* Numerator smaller than denominator:
    ** quotient is zero, remainder is numerator
    */
    if(quo) numZero(quo);
    if(rem) numCopy(rem, num);
    return;
  }

  /* Allocate space for temp space variables */
  numAlloc(&un, num->len + 1, num->base);
  numAlloc(&tmp, den->len + 1, num->base);

  /* Copy numerator into temp space */
  numCopy(&un, num);

  /* uns is used to work with a few units of the un variable */
  uns.len = tmp.len;
  uns.base = tmp.base;

  /* ji is the index of the left (more significant) unit of the numerator */
  for(ji = lnz_num + 1; ji >= lnz_den+1; ji --) {

    /* -- Calculate qh -- */
    qh = (un.units[ji] * num->base + un.units[ji-1]) / den->units[lnz_den];
    if(qh >= num->base)
      qh = num->base - 1;

    /* -- Multiply and subtract -- */
multiply:
    /* Multiply to form subtrahend */
    numZero(&tmp);
    numMultPrimUnit(qh, den, &tmp, 0);

    /* Check to see if trial quotient was too big */
    uns.units = &un.units[ji-lnz_den-1];

#ifdef DEBUG_LD
    printf("ji=%i  lnz_num=%i  lnz_den=%i\n", ji, lnz_num, lnz_den);
    printf("qh=%i  un = %i %i  den = %i\n", qh, un.units[ji], un.units[ji-1],
           den->units[lnz_den]);
    printf("uns      = "); numPrint(&uns);
    printf("qh * den = "); numPrint(&tmp);
#endif

    if(numCompare(&uns, &tmp) < 0) {
      qh--;
#ifdef DEBUG_LD
      printf("loop\n");
#endif
      goto multiply;  /* See note below about the use of the dreaded goto */
    }

    /* Subtract */
    numSubPrim(&uns, &tmp, &uns);

    /* -- Store quotient -- */
    if(quo) quo->units[ji-lnz_den-1] = qh;

#ifdef DEBUG_LD
    printf("new uns  = "); numPrint(&uns);
    if(quo) { printf("quo = "); numPrint(quo); }
#endif

  }

  /* Remainder is least significant lnz_den+1 units */
  if(rem)
    memcpy(rem->units, un.units, sizeof(NumberUnitT) * (lnz_den+1));

#ifdef DEBUG_LD
    printf("final rem  = "); numPrint(rem);
    printf("final uns  = "); numPrint(&uns);
    printf("final un  = "); numPrint(&un);
    if(quo) { printf("final quo = "); numPrint(quo); }
#endif

  numFree(&un);
  numFree(&tmp);
}
/* Note to the grader about the use of the goto above:
**
** I've thought of 6 ways to rewrite this routine in order to avoid
** the use of the "goto" because it's generally considered to be bad form
** to use a goto.  I usually agree, and this is probably the third time
** I've ever used a goto in a C program.  However, all of the other ways
** of rewriting this loop, that do not use the goto, are less elegant.
**
** I've analysed the problem, and I think that the most concise
** explanation for why /not/ using a goto in this case makes the
** routine less elegant is this:  The situation handled by the goto
** statement is a sort of "exception".  The "goto" is executed in a case
** where the algorithm decided that it made a mistake at an earlier
** step.  Algorithms that don't make guesses generally do not have to
** intercept themselves and restart their process, which is why "goto"
** statements are rare.  But I think that in this case, its use is
** justified.
*/




/* numMultiply: Multiply two NumberT's, MODULUS mod
**
** return prod = n1 * n2 MODULUS mod
**
** Assumes n1->len == n2->len == prod->len.
** Assumes n1->base == n2->base.
**
** Input values are NOT assumed to be modular.
**
** The large intermediate result modulus is found by using a
** long division routine to get the remainder.
*/
void
numMultiply(const NumberT *n1,const NumberT *n2,NumberT *prod,const NumberT *mod)
{
  NumberT prod_big;  /* temp space -- holds large intermediate result */

  /* Allocate space for big product, which will in general have twice
  ** as many digits as the multiplicands
  */
  numAlloc(&prod_big, n1->len * 2 + 2, n1->base);

  /* Perform the primative (non-modular) multiply */
  numMultPrim(n1, n2, &prod_big);

  /* Find the modulus of the intermediate product */
  numLongDiv(&prod_big, mod, NULL, prod);

  /* Free the temp space */
  numFree(&prod_big);
}




/* NUMBER_LENGTH is the number of primatives per NumberT */
#define NUMBER_LENGTH 20


/* Read in some numbers from a file whose name is given on the command
** line.  Then perform all of the modular arithmetic routines that have
** been implemented, on those numbers, in each direction, to demonstrate
** the code.
*/
int
main(int argc, char **argv)
{
  FILE *fileP;
  char *filename;
  NumberT n1;
  NumberT n2;
  NumberT result;
  NumberT xn;
  NumberT yn;
  NumberT mod;
  NumberT max;
  int ni;
  int number_base;


  /* Find a number base that is a power of ten that fits into the
  ** primative Number Unit type.
  */
  number_base = pow(10,(int)log10(NUMBER_BASE_MAX));
  printf("number_base is %i\n", number_base);


  /* Allocate a bunch of numbers */
  numAlloc(&n1, NUMBER_LENGTH, number_base);
  numAlloc(&n2, NUMBER_LENGTH, number_base);
  numAlloc(&result, NUMBER_LENGTH, number_base);
  numAlloc(&xn, NUMBER_LENGTH, number_base);
  numAlloc(&yn, NUMBER_LENGTH, number_base);
  numAlloc(&mod, NUMBER_LENGTH, number_base);
  numAlloc(&max, NUMBER_LENGTH, number_base);


  /* Find the maximum representable number, just for reference */
  for(ni=0; ni < max.len; ni++) {
    mod.units[ni] = max.units[ni] = max.base - 1;
  }
  printf("max       = ");  numPrintContiguous(&max);


  /* Parse command line arguments */
  if(argc != 2) {
    fprintf(stderr, "usage: %s <filename>\n", argv[0]);
    exit(1);
  }

  filename = argv[1];

  if((fileP=fopen(filename, "r"))==NULL) {
    fprintf(stderr, "%s:err: file '%s' could not be opened\n",
    argv[0], filename);
    exit(1);
  }


  /* Read the modulus from the file */
  numRead(&mod, fileP, NULL);
  printf("mod       = ");
  numPrintContiguous(&mod);

  /* Read the first number from the file */
  numRead(&n1, fileP, &mod);
  printf("n1        = ");
  numPrintContiguous(&n1);

  /* Read the second number from the file */
  numRead(&n2, fileP, &mod);
  printf("n2        = ");
  numPrintContiguous(&n2);

  fclose(fileP);


  /* Perform every possible operator in both directions */

#ifdef ORIGINAL_TESTS
  numAdd(&n1, &n2, &result, &mod);
  printf("\nn1 (+) n2 = ");
  numPrintContiguous(&result);

  numAdd(&n2, &n1, &result, &mod);
  printf("n2 (+) n1 = ");
  numPrintContiguous(&result);

  numSubtract(&n1, &n2, &result, &mod);
  printf("\nn1 (-) n2 = ");
  numPrintContiguous(&result);

  numSubtract(&n2, &n1, &result, &mod);
  printf("n2 (-) n1 = ");
  numPrintContiguous(&result);

  numMultiply(&n1, &n2, &result, &mod);
  printf("\nn1 (*) n2 = ");
  numPrintContiguous(&result);

  numMultiply(&n2, &n1, &result, &mod);
  printf("n2 (*) n1 = ");
  numPrintContiguous(&result);
#endif /* ORIGINAL_TEST */

#include "main-update.c"

  /* Delete the numbers */
  numFree(&n1);
  numFree(&n2);
  numFree(&result);
  numFree(&xn);
  numFree(&yn);
  numFree(&max);
  numFree(&mod);

  /* This is to make Unix happy */
  return 0;
}