[Toybox] Long division algorithm for the toybox project by Graff

[CM Graff]

// Fixed-point arbitrary precision long hand division algorithm written by
// Christopher M. Graff for Rob Landley's `toybox' project.
//
// A commonly known multiplication algorithm is included to help test the
// effectiveness of the fixed-point system.
//
// The division algorithm needs a large hardware type such as int in order to
// hold carried values. The multiplication algorithm uses "partial carries"
// so it will work with extremely small types such as char. As far as I know,
// no "partial carry" system has been invented for long division.

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

typedef struct {    // Toym fxdpnt type
  int *number;      // The actual number
  int sign;         // Sign
  size_t fp;        // Length left of radix
  size_t len;       // Length of number (count of digits / limbs)
  size_t allocated; // Length of allocated memory
  size_t chunk;     // Allocation chunk size
} fxdpnt;

fxdpnt *toym_div(fxdpnt *, fxdpnt *, fxdpnt *, int base);
fxdpnt *toym_mul(fxdpnt *, fxdpnt *, fxdpnt *, int base);
fxdpnt *toym_expand(fxdpnt *, size_t);
fxdpnt *toym_str2fxdpnt(const char *);
fxdpnt *toym_add_precision(fxdpnt *, size_t);
fxdpnt *toym_alloc(size_t);
void *toym_malloc(size_t);
void *toym_realloc(void *, size_t);
void toym_flipsign(fxdpnt *);
void toym_setsign(fxdpnt *, fxdpnt *, fxdpnt *);
void toym_init(fxdpnt *);
void toym_print(fxdpnt *);
void toym_free(fxdpnt *);
void toym_error(char *);
void toym_copyarr(int *, int *, size_t);
void toym_setarr(int *, int, size_t);
int toym_highbase(int);

// helper functions that could be replaced or hard coded
int toym_highbase(int a)
{
  // Handle high bases
  static int glyphs[65] = { '0', '1', '2', '3', '4', '5', '6', '7', '8',
        '9', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h',
        'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q',
        'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z',
        'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I',
        'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R',
        'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '+',
        '/', '='};
  if ( a <= 64 )
    return glyphs[a];
  else // just use the ascii values for bases that are very high
    return a;
}

fxdpnt *toym_alloc(size_t len)
{
  // Allocate the basic requirements of a toym `fxdpnt'
  fxdpnt *ret = toym_malloc(sizeof(fxdpnt));
  ret->number = toym_malloc(sizeof(int) * len);
  ret->sign = '+';
  ret->fp = 0;
  ret->allocated = len;
  ret->len = 0;
  ret->chunk = 4; // set to 4 to force worst case tests, change to >255
  return ret;
}

void toym_init(fxdpnt *flt)
{
  flt->sign = '+';
  flt->len = flt->fp = 0;
}

void toym_flipsign(fxdpnt *flt)
{
  if (flt->sign == '+')
    flt->sign = '-';
  else if (flt->sign == '-')
    flt->sign = '+';
}

void toym_print(fxdpnt *flt)
{
  size_t i = 0;
  if (flt->sign == '-')
    putchar(flt->sign);
  for (i = 0; i < flt->len ; ++i){
    if (flt->fp == i)
      putchar('.');
    putchar(toym_highbase((flt->number[i])));
  }
  putchar('\n');
  fflush(stdout);
}

fxdpnt *toym_str2fxdpnt(const char *str)
{
  // Convert a string to a toym `fxdpnt'
  size_t i = 0;
  int flt_set = 0, sign_set = 0;

  fxdpnt *ret = toym_expand(NULL, 1);

  for (i = 0; str[i] != '\0'; ++i){
    if (str[i] == '.'){
      flt_set = 1;
      ret->fp = i - sign_set;
    }
    else if (str[i] == '+'){
      sign_set = 1;
      ret->sign = '+';
    }
    else if (str[i] == '-'){
      sign_set = 1;
      ret->sign = '-';
    }
    else{
      ret = toym_expand(ret, ret->len + 1);
      ret->number[ret->len++] = str[i] - '0';
    }
  }

  if (flt_set == 0)
    ret->fp = ret->len;

  return ret;
}

void toym_free(fxdpnt *flt)
{
  if (flt->number)
    free(flt->number);
  free(flt);
}

fxdpnt *toym_expand(fxdpnt *flt, size_t request)
{
  // Enlarge or create a fxdpnt
  if (flt == NULL){
    flt = toym_alloc(request);
  } else if (request > flt->allocated){
    flt->allocated = (request + flt->chunk);
    flt->number = toym_realloc(flt->number, flt->allocated * sizeof(int));
  }
  return flt;
}

fxdpnt *toym_add_precision(fxdpnt *flt, size_t more)
{
  // Increase the precision of a toym `fxdpnt'
  flt = toym_expand(flt, flt->len + more);
  toym_setarr(flt->number + flt->len, 0, more);
  flt->len += more;
  return flt;
}

void toym_error(char *message)
{
  fprintf(stderr, "%s\n", message);
  exit(1);
}

void *toym_malloc(size_t len)
{
  void *ret;
  if(!(ret = malloc(len)))
    toym_error("malloc failed\n");
  return ret;
}

void *toym_realloc(void *ptr, size_t len)
{
  void *ret;
  if(!(ret = realloc(ptr, len)))
    toym_error("realloc failed\n");
  return ret;
}

void toym_copyarr(int *answer, int *from, size_t len)
{ // this could be replaced with memcpy
  size_t i = 0;
  for( i = 0; i < len ; i++)
    answer[i] = from[i];
}

void toym_setarr(int *answer, int delim, size_t len)
{ // this could be replaced with memset
  size_t i = 0;
  for( i = 0; i < len; i++)
    answer[i] = delim;
}

void toym_setsign(fxdpnt *a, fxdpnt *b, fxdpnt *c)
{
  toym_init(c);
  if (a->sign == '-')
    toym_flipsign(c);
  if (b->sign == '-')
    toym_flipsign(c);
}

fxdpnt *toym_mul(fxdpnt *a, fxdpnt *b, fxdpnt *c, int base)
{
  int i = 0, j = 0, sum = 0, carry = 0;
  size_t k = 0;

  toym_setsign(a, b, c);
  c = toym_expand(c, a->len + b->len);
  toym_setarr(c->number, 0, a->len + b->len);

  for ( i = a->len - 1; i >= 0 ; i--){
    for ( j = b->len - 1, k = i + j + 1, carry = 0; j >= 0 ; j--, k--){
      sum = (a->number[i]) * (b->number[j]) + (c->number[k]) + carry;
      carry = sum / base;
      c->number[k] = (sum % base);
    }
    c->number[k] += carry;
    c->len++;
  }
  c->fp = a->fp + b->fp;
  return c;
}

fxdpnt *toym_div(fxdpnt *a, fxdpnt *b, fxdpnt *c, int base)
{
  size_t i = 0, j = 0, z = 0, diff = 0, off = 0;
  size_t width = a->len + b->len;
  int sum = 0, rec = 0;
  // temporary variables for division "guesses"
  int *mir = toym_malloc(sizeof(int) * width);
  int *tmir = toym_malloc(sizeof(int) * width);

  toym_setsign(a, b, c);

  c = toym_expand(c, a->len + b->len);
  toym_setarr(mir + a->len, 0, width - a->len);
  toym_copyarr(mir, a->number, a->len);
  c->number[z] = 0;

  // count the zeros before a non-zero value
  while (b->number[off] == 0 && off < b->len)
    ++off;

  // precompute the position of the radix for "c"
  if (a->fp < b->fp){
    diff = b->fp - a->fp - 1;
    toym_setarr(c->number, 0, diff);
    c->len = z = diff;
    c->fp = off;
    c->number[z] = 0;
  }
  else if (a->fp + 1> b->fp)
    c->fp = a->fp - b->fp + off + 1;
  else
    c->fp = a->fp + off;

  // this is the actual long division routine which operates on two arrays
  for (;z < a->len + diff;){
    for (rec = 0, i = off, j = z - diff;i < b->len;j++, i++){
      sum = (mir[j]) - (b->number[i]);
      if (sum < 0){
        if (j == z - diff){
          mir[j + 1] += ((mir[j]) * base);
          z++;
          c->len++;
          c->number[z] = 0;
        }
         else{
          mir[j - 1] -= 1;
          mir[j] += base;
        }
        rec = 1;
        break;
      }
      tmir[j] = sum;
    }
    if (rec == 0){
      toym_copyarr(mir, tmir, j);
      c->number[z] += 1;
    }
  }
  free(mir);
  free(tmir);
  return c;
}

int main(int argc, char *argv[])
{
  if (argc < 3 )
    toym_error("Needs 2 numbers");

  size_t scale = 100;
  int base = 10;
  fxdpnt *a, *b, *c;

  a = toym_str2fxdpnt(argv[1]);
  b = toym_str2fxdpnt(argv[2]);
  a = toym_add_precision(a, scale);
  b = toym_add_precision(b, scale);
  c = toym_expand(NULL, 1);

  // multiplication
  c = toym_mul(a, b, c, base);
  toym_print(c);
  // division
  c = toym_div(a, b, c, base);
  toym_print(c);

  toym_free(a);
  toym_free(b);
  toym_free(c);

  return 0;
}