Listing 1 (facbase.c) Arithmetic in Factorial-Base

#include <stdio.h>

#ifndef TRUE
#define TRUE 1
#endif
#ifndef FALSE
#define FALSE 0
#endif

#define DEBUG TRUE  /* include RPN calculator */
#define PLUS FALSE  /* sign values */
#define MINUS TRUE
#define MAXFBINTEGER 100    /* <= 180 for arrays of ints */
#define MAXFBFRACTION 100   /* <= 180 for arrays of ints */
#define HIGHGUARD MAXFBINTEGER + 1
#define LOWGUARD MAXFBFRACTION + 1
#define ARRAYSIZE 1 + HIGHGUARD + LOWGUARD

/* two constants in the proper format */
int zero[ARRAYSIZE] = { PLUS, 0 }; /* signed 0 */
int one[ARRAYSIZE] = { PLUS, 1, 0 };

/* a flag set by divide() while "estimating" */
int nowarning = FALSE;
/* assign the value of one array to a second array */
assignto(dest, source)
int *dest;
int *source;
{
   int i;

   for(i = 0; i < ARRAYSIZE; i++)
       dest[i] = source[i];
{
/* z = abs(x) */
absolute(x, z)
int *x;
int *z;
{
   if(z != x)       /* if not changing in place */
       assignto(z, x); /* copy */
   z[0] = PLUS;     /* force positive */
}

/* Assign unary negative of x to z */
negative(x, z)
int *x;
int *z;
{
   int sign;

   sign = (x[0] == PLUS) ? MINUS : PLUS;
   if(z != x)       /*if not changing in place */
       assignto(z, x); /* copy */
   z[0] = sign;     /* change sign */
}

/* compare 2 factorial-base numbers for |x| > |y| */
int greaterthan(x, y)
int *x;
int *y;
{
   int i, j;

   /*
    * check the integer part first,
    * largest subscript to smallest
    */
   for(i = MAXFBINTEGER; (x[i] == y[i]) && (i > 1); -i)
       ;
   if(i != 0)
   {
       if(x[i] > y[i])
           return(TRUE);
       if(x[i] < y[i])
           return(FALSE);
   }
   /*
    * The integer portions are equal,
    * continue with the fractional part,
    * smallest subscript to largest.
    */
   for(i = HIGHGUARD + 1, j = 1;
       (x[i] == y[i]) && (j < MAXFBFRACTION - 1);
       i++, j++)
           ;
   return((x[i] > y[i]) ? TRUE : FALSE);
}

/* Compare 2 factorial-base numbers for |x| < |y| */
int lessthan(x, y)
int *x;
int *y;
{
   int i, j;

   /*
    * check the integer part first,
    * largest subscript to smallest
    */
   for(i = MAXFBINTEGER; (x[i] == y[i]) && (i > 1); -i)
       ;
   if(i != 0)
   {
       if(x[i] < y[i])
           return (TRUE);
       if(x[i] > y[i])
           return(FALSE);
   }
   /*
    * The integer portions are equal,
    * continue with the fractional part,
    * smallest subscript to largest.
    */
   for(i = HIGHGUARD + 1, j = 1;
       (x[i] == y[i]) && (j < MAXFBFRACTION - 1);
       i++, j++)
           ;
   return((x[i] < y[i]) ? TRUE : FALSE);
}

/* Compare 2 factorial-base numbers for |x| == |y| */
int equalto(x, y)
int *x;
int *y;
{
   int i;

   /*
    * Check the whole array except the sign
    * - order doesn't matter.
    */
   for(i = ARRAYSIZE-1; (x[i] == y[i]) && (i > 1); -i)
       ;
   return((x[i] == y[i]) ? TRUE : FALSE);
}

/*
 * "Normalize" a factorial-base number. All of the
 * arithmetic functions call this routine to handle
 * carrys and borrows. A factorial-base number has a
 * proper form where every factorial position in its
 * integer part has a value between 0 and the the
 * magnitude of its position and every factorial position
 * in its fractional part has a value between 0 and one
 * less than the magnitude of its position.
 */
normalize(n)
int *n;
{
   int i, j;
   int *x;

   x = &n[HIGHGUARD];
   /*
    * First, check for loss of precision
    * during multiplication or division.
    */
   if(x[LOWGUARD])
   {
       if(nowarning != TRUE)
           fputs("UNDERFLOW\n", stderr);
       x[LOWGUARD] = 0;
   }
   /*
    * Now, work the fractional part first,
    * largest subscript to smallest.
    * The subscript j is the factorial position
    * being put into proper form.
    */
   for(j = MAXFBFRACTION, i = j - 1; i >= 1; -i, -j)
   {
       /* if(x[j] >= j) carry */
       x[i] += (x[j] / j);
       x[j] %=j;
       if(x[j] < 0)    /* borrow */
       {
           x[i] -= 1;
           /* modulo= j */
           x[j] += j;  /* make positive */
       }
   }
   /* shift any carry to integer part & clear carry */
   n[1] += x[1];
   x[1] = 0;
   /*
    * Now, normalize the integer part,
    * working from smallest subscript to largest.
    * The subscript i is the factorial position
    * being put into proper form.
    */
   x = n;
   for(i = 1, j = 2; i <= MAXFBINTEGER; i++, j++)
   {
       /* if(x[i] >= j) carry */
       x[j] += (x[i] / j);
       x[i] %= j;
       if(x[i] < 0)    /* borrow */
       {
           x[j] -= 1;
           x[i] += j;
       }
   }
   if(x[i])     /* if an entry in x[HIGHGUARD] */
   {
       fputs("OVERFLOW\n", stderr);
       x[i] = 0;
   }
}

/* Add y to x, put result in z */
add(x, y, z)
int *x;
int *y;
int *z;
{
   int sign, i;
   int copy[ARRAYSIZE];

   if(x[0] != y[0])   /* if different signs */
   {
       if(y[0] == MINUS)
       {
           /*
            * Change the sign of y
            * and subtract y from x.
            */
           negative(y, copy);
           subtract(x, copy, z);
       }
       else
       {
           /*
            * Change the sign of x
            * and subtract x from y.
            */
           negative(x, copy);
           subtract(y, copy, z);
       }
   }
   else
   {
       sign = x[0];    /* save the sign */
       for(i = ARRAYSIZE - 1; i > 0; -i)
           z[i] = x[i] + y[i];
       z[0] = sign;
       normalize(z);
   }
}

/* Subtract y from x, put result in z */
subtract(x, y, z)
int *x;
int *y;
int *z;
{
   int sign, i;
   int copy[ARRAYSIZE];

   if(x[0] != y[0])    /* if signs are different */
   {
       negative(y, copy);  /* change sign of y */
       sign = x[0];        /* save sign of x */
       add(x, copy, z);
       z[0] = sign;
       return;
   }
   else if(y[0] == MINUS)  /* (-x) - (-y) */
   {
       /* if(abs(y) < abs(x)) sign = MINUS */
       sign = lessthan(y, x);
   }
   else
   {
       /* if(x < y) sign = MINUS */
       sign = lessthan(x, y);
   }
   /* Subtract based on the absolute values */
   if (lessthan(x, y))
   {
       for(i = ARRAYSIZE - 1; i > 0; -i)
           z[i] = y[i] - x[i];
   }
   else
   {
       for(i = ARRAYSIZE - 1; i > 0; -i)
           z[i] = x[i] - y[i];
   }
   z[0] = sign;
   normalize(z);
}

/*
 * Multiply factorial-base x by integer y, result in z.
 * Utility routine called by multiply(), divide(),
 * atofact(), fractoa(), and facttoa(), and always
 * with a positive value for y
 */
multfbyi(x, y, z)
int *x;
int y;
int *z;
{
   int i;

   for(i = ARRAYSIZE - 1; i > 0; -i)
       z[i] = x[i] * y;
   normalize(z);
}

/*
 * Divide factorial-base x by integer y, result in z.
 * Utility routine called by multiply(), divide(),
 * and facttoa(), and always with a positive y
 */
divfbyi(x, y, z)
int *x;
int y;
int *z;
{
   int i, j, carry, part;

   carry = 0;
   /*
    * Work the integer part first,
    * from the largest subscript to smallest
    */
   for(i = MAXFBINTEGER, j = i + 1; i >= 1; -i, -j)
   {
       part = x[i] + carry * j;
       carry = part % y;
       z[i] = part / y;
   }
   /*
    * Now, work the fractional part,
    * from the smallest subscript to largest
    */
   for(i = HIGHGUARD + 1, j = 1;
       j <= MAXFBFRACTION; i++, j++)
   {
       part = x[i] + carry * j;
       carry = part % y;
       z[i] = part / y;
   }
   /*
    * propogate any carry into z[LOWGUARD]
    * to mark underflow and loss of precision
    */
   z[i] = carry * j;
   normalize(z);
}

/*
 * Multiply factorial-base x by factorial-base y,
 * assign result to z. Uses the identity
 * number * (integer + fraction)
 * == (number * integer) + (number * fraction)
 */
multiply(x, y, z)
int *x;
int *y;
int *z;
{
   int i, j, k, sign;
   int partial [ARRAYSIZE];
   int temp[ARRAYSIZE];
   int copy[ARRAYSIZE];

   if(x[0] != y[0])    /* if signs different */
       sign = MINUS;
   else
       sign = PLUS;

   assignto(partial, zero);   /* Initialize result */
   /*
    * Work the integer portion first,
    * from smallest subscript to largest
    */
   absolute(x, copy);  /* copy = abs(x) */
   /* first, find largest subscript k where y[k] != 0 */
   for(k = MAXFBINTEGER; (k > 0) && (y[k] == 0); -k)
       ;
   for(i = 1; i <= k; i++)
   {
       /* first shift copy by factorial position */
       multfbyi(copy, i, copy);
       if(y[i])    /* don't bother multiplying by 0 */
       {
           /* multiply by factorial digit */
           multfbyi(copy, y[i], temp);
           add(partial, temp, partial);
       }
   }

   /* now work fraction part */
   assignto(copy, x); /* reset copy */
   /* find largest subscript k where y[k] != 0 */
   for(k = ARRAYSIZE - 1; (k > HIGHGUARD + 1) && (y[k] == 0); -k)
       ;
   for(i = HIGHGUARD + 2, j = 2; i <= k; i++, j++)
   {
       /* first shift copy by factorial position */
       divfbyi(copy, j, copy);
       if(y[i])    /* don't bother multiplying by zero */
       {
           /* multiply by factorial digit */
           multfbyi(copy, y, temp);
           add(partial, temp, partial);
       }
   }
   partial[0] = sign;
   assignto(z, partial);
}
/*
 * Divide factorial-base x by factorial-base y, store
 * result in z. Uses blackboard style long division.
 */
divide(x, y, z)
int *x;
int *y;
int *z;
{
   int i, j, sign;
   int estimate;
   int copyx[ARRAYSIZE];
   int copyy[ARRAYSIZE];
   int temp[ARRAYSIZE];
   int partial[ARRAYSIZE];

   if(x[0] != y[0])    /* if signs different */
       sign = MINUS;
   else
       sign = PLUS;
   absolute(x, copyx);
   absolute(y, copyy);
   assignto(partial, copyx);
   assignto(temp, copyy);
   /*
    * First, estimate the integer part of result by
    * driving y to 1.xxx. Division is VERY slow, so
    * the extra time spent to identify special cases
    * is well worth it.
    */
   if(equalto(temp, zero))
   {
       /* division by zero fault */
       /* not handled here */
       return;
   }
   else if(lessthan(partial, temp))
   {
       assignto(partial, zero);    /* integer part 0 */
   }
   else if(lessthan(one, temp))
   {
       /*
        * This could cause a spurious UNDERFLOW
        * message even though the final result
        * would be exact, so we set a flag to
        * suppress the warning.
        */
       nowarning = TRUE;
       while(lessthan(one, temp))
       {
           divfbyi(partial, 2, partial);
           divfbyi(temp, 2, temp);
       }
       multfbyi(partial, 2, partial);
       nowarning = FALSE;  /* reset flag */
   }
   else if(lessthan(temp, one))
   {
       while(lessthan(temp, one))
       {
           multfbyi(partial, 2, partial);
           multfbyi(temp, 2, temp);
       }
   }
   else     /* division by 1 or -1 */
   {
       assignto(z, x);
       z[0] = sign;
       return;
   }
   /* Now, delete fractional part of estimate */
   for(i = HIGHGUARD + 1; i < ARRAYSIZE; i++)
       partial[i] = 0;
   multiply(copyy, partial, temp);
   while(greaterthan(temp, copyx))
   {
       subtract(partial, one partial);
       subtract(temp, copyy, temp);
   }
   subtract(copyx, temp, copyx);
   multfbyi(copyx, 2, copyx);
   /* partial now holds integer part of result */

   /*
    * Now, proceed by long division to divide by
    * the fractional part - using the subscript
    * (less 1) as the estimate at each position
    */
   for(i = HIGHGUARD + 2, j = 2;
       !equalto(copyx, zero) && (i < ARRAYSIZE);
       i++)
   {
       estimate = j - 1;
       do {
           multfbyi(copyy, estimate-, temp);
       } while(greaterthan(temp, copyx));
       subtract(copyx, temp, copyx);
       partial[i] = ++estimate;
       multfbyi(copyx, ++j, copyx);
   }
   normalize(partial);
   if(sign == MINUS)
       partial[0] = MINUS;
   assignto(z, partial);
}

/*
 * ASCII to factorial-base conversion
 * Just like ASCII to binary conversion!
 */
atofact(s, z)
char s[];
int *z;
}
   int i, j, sign;
   int ipart[ARRAYSIZE];
   int fpart[ARRAYSIZE];

   assignto(ipart, zero);
   assignto(fpart, zero);
   i = 0;
   sign = PLUS;
   while(s[i] == ' ' || s[i] == '\t')
       i++;
   if(s[i] == '-' || s[i] == '+')
   {
       if(s[i++] == '-' )
           sign = MINUS;
   }
   for(; s[i] >= '0' && s[i] <= '9'; i++)
   {
       multfbyi(ipart, 10, ipart);
       ipart[1] = s[i] - '0';
       normalize(ipart);
   }
   if(s[i] == '.')
   {
       i++;
       j = 0;
       for( ; s[i] >= '0' && s[i] <= '9'; i++)
       {
           multfbyi(fpart, 10, fpart);
           ++j;
           fpart[1] = s[i] - '0';
           normalize(fpart);
       }
       while(j-)
           divfbyi(fpart, 10, fpart);
       add(ipart, fpart, ipart);
   }
   ipart[0] = sign;
   assignto(z, ipart);
}

/*
 * Convert the fractional part of x to ASCII.
 * Up to count characters go to stdout.
 */
fractoa(x, count)
int *x;
int count;
{
   int i;
   int temp[ARRAYSIZE];

   assignto(temp, x);
   for(i = 1; i <= MAXFBINTEGER; i++)
       temp[i] = 0;   /* erase integer part */
   temp[0] = PLUS;     /* always positive */
   if(equalto(temp, zero))
       return; /* no fractional part to print out */
   putchar('.');
   while(count- && !equalto(temp, zero))
   {
       multfbyi(temp, 10, temp);
       putchar('0' + 6 * temp[3] + 2 * temp[2] + temp[1]);
       /* Now erase the integer part */
       temp[3] = temp[2] = temp[1] = 0;
   }
}

/*
 * CAUTION - altering the size of outbuff requires
 * some art. If MAXFBINTEGER == 100, it must be
 * large enough to hold the 160 decimal digit integer
 * 9.4259 * 10^159. If MAXFBINTEGER == 180, it must
 * be large enough for the 332 decimal digit integer
 * 3.6362 * 10^331. If you want to deal with really
 * big numbers and increase MAXFBINTEGER, you'll have
 * to give some thought as to how large the conversion
 * buffer is going to have to be.
 */
/* Allow a little slack */
char outbuff[MAXFBINTEGER*2] = { 0 };
int outptr; /* Actually an index and not a "ptr" */

/*
 * Factorial-base to ASCII conversion, integer
 * part end up to count characters of fractional
 * portion go to stdout.
 */
facttoa(x, count)
int *x;
int count;
{
   int i, j, sign;
   int temp[ARRAYSIZE];
   int val[ARRAYSIZE];

   outptr = 0;
   assignto(val, zero);
   assignto(temp, x);
   if((sign = temp[0]) == MINUS)
   {
       temp[0] = PLUS;
       putchar('-');
   }
   for(i = ARRAYSIZE - 1; i > HIGHGUARD; -i)
       temp[i] = 0;    /* erase fractional part */
   while(!equalto(temp, zero))
   {
       divbyi(temp, 10, temp);
       for(i = HIGHGUARD + 2, j = 1; j <= 4; i++, j++)
       {
           val[i] = temp[i];
       }
       multfbyi(val, 10, val);
       outbuff[outptr++] =
           '0' + 6 * val[3] + 2 * val[2] + val[1];
       val[3] = val[2] = val[1] = 0;
       /*Now erase fractional portion of temp */
       temp[i-1] = temp[i-2] = temp[i-3] = temp[i-4] = 0;
   }
   if (outptr == 0) /* if no integer part */
       putchar ('0');
   else
   {
       while(outptr-)
           putchar(outbuff[outptr]);
   }
   fractoa(x, count);  /* to print fractional portion */
   putchar('\n');
}

*/ Remainder of file is RPN calculator & display */
#ifdef DEBUG
/*
 * Print a factorial-base number in factorial-base.
 * "digits" are printed between '<' and '>',
 * output goes to stdout.
 */
facprint(x)
int *x;
{
   int i, j;
   if(x[0] == MINUS)
       printf("-");
   /* Delete any leading zeroes */
   for(i = MAXFBINTEGER; i >= 1; i-)
   {
       if(x[i] != 0)
           break;
   }
   /* Print any integer portion */
   for( ; i >= 2; )
       printf("<%d>", x[i-]);
   /* Make sure to print at least one digit */
   printf("<%d>", x[1]);
   printf(".");
   /*
    * Print fractional part, deleting any trailing
    * zeroes but printing at least one digit
    */
   i = HIGHGUARD + 2; /* start at 2! */
   printf("<%d>", x[i++]);
   for(j = 0; i < ARRAYSIZE; i++)
   {
       if(x[i] == 0)
           j += 1;
       else
       {
           while(j)
           {
printf("<0>");

-j;

           }
           printf("<%d>", x[i]);
       }
   }
   printf("\n");
}

/*
 * A simple but VERY slow reverse Polish calculator.
 * Commands +, -, *, /, D to print in decimal,
 * F to print in factorial, C to clear the stack,
 * S to display the whole stack in decimal,
 * a decimal number to use as an operand,
 * only 1 operator or operand per line!!
 */

#define IPSIZE 256
char input[IPSIZE];
#define STACKSIZE 8
int stack[STACKSIZE][ARRAYSIZE];

main(argc, argv)
int argc;
char *argv[];
{
   int x, i, prompt, depth;

   prompt = (argc > 1) ? TRUE : FALSE;
   depth = 0;
   for( ; ; )
   {
       if(prompt)
           printf(%d> ", depth);
       if(fgets(input, IPSIZE, stdin) == NULL)
           break;
       switch(x = input[0])
       {
       case 'c':   /* clear the stack */
       case 'C':   depth = 0;
                 continue;
       case 'f':   /* print top of stack in factorial */
       case 'F':   if(depth < 1)
                 {
                     printf("empty stack\n");
                     continue;
                 }
                 printf("%d: ", depth - 1);
                 facprint(&stack[depth-1][0]);
                 continue;
       case'd':    /* print top of stack in decimal  */
       case 'D':   if(depth < 1)
                 {
                     printf("empty stack\n");
                     continue;
                 }
                 printf("%d: ", depth - 1);
                 facttoa(&stack[depth-1] [0], 50);
                 continue;
       case's':    /* display contents of stack  */
       case 'S':   if(depth < 1)
                 {
                     printf("empty stack\n");
                     continue;
                 {
                 for(i = 0; i < depth; i++)
                 }
                     printf("%d: ", i);
                     facttoa(&stack[i][0], 50);
                 }
                 continue;
       case '+':   if(depth < 2)
                 }
                     printf("stack will underflow\n");
                     continue;
                 {
                 add (&stack [depth-2] [0],
                     &stack [depth- 1] [0],
                     &stack[depth-2] [0] );
                 -depth;
                 continue;
       case '/':   if(depth < 2)
                 {
                     printf ("stack will underflow\n");
                     continue;
                 }
                 if (equal to(&stack [depth-1][0], zero))
                 {
                     printf ("division by zero\n");
                     /*  allow the 0 to be discarded  */
                 }
                 else
                 {
                     divide (&stack [depth-2][0],
                         &stack [depth-1][0],
                         &stack[depth-2][0]);
                 }
                 -depth;
                 continue;
       case '*':   if(depth < 2)
                 {
                     printf ("stack will underflow\n");
                     continue;
                 }
                 multiply(&stack [depth-2][0],
                     &stack [depth-1] [0],
                     &stack[depth-2] [0]);
                 -depth;
                 continue;
       case '-':   if(input[1] == '\n')
                 {
                     if(depth < 2)
                     {
                         printf(
                          "stack will underflow\n");
                         continue;
                     }
                     subtract(&stack[depth-2][0],
                         &stack[depth-1][0],
                         &stack(depth-2][0]);
                     -depth;
                     continue;
                 }    /* else a negative number */
                 else
                     break;
       default:    if(x != '-' && x != '.' &&
                     (x < 'o' || x > '9'))
                 {
                     printf("invalid entry\n");
                     continue;
                 }
                 break;
                 /* to convert and stack a number */
       }
       if(depth >= STACKSIZE)
       {
           printf("stack will overflow\n");
           continue;
       }
       atofact(input, &stack[depth++][0]);
       continue;
   }
}
#endif

/* End of File */