div.c

#define ITERATIONS_INV 2
void ufp_inv(unsigned long *val) {
  unsigned long guess[PRECISION], dblbuf[PRECISION*2],
    sglbuf[PRECISION], two[PRECISION+1];
  int i, j;

  two[0] = HI_BIT;
  for(i=1; i<PRECISION+1; i++)
    two[i] = 0;
  guess[0] = invseed_x(val[0]);
  for(i=1; i<PRECISION; i++)
    guess[i] = 0;

  for(i=0; i<ITERATIONS_INV; i++) {
    mulu_x(val, guess, dblbuf);
    lsr_x(PRECISION+1, dblbuf);
    neg_x(PRECISION+1, dblbuf);
    add_x(PRECISION+1, dblbuf, two);
    lsl_x(PRECISION+1, dblbuf);

    for(j=0; j<PRECISION; j++)
      sglbuf[j] = dblbuf[j];

    mulu_x(guess, sglbuf, dblbuf);
    lsl_x(PRECISION+1, dblbuf);

    for(j=0; j<PRECISION; j++)
      guess[j] = dblbuf[j];
  } for(i=0; i<PRECISION; i++)
    val[i] = guess[i];
}

void inv(num_p num) {
  int i;

  if (num->fmt == FMT_INT) {
    convert_flt(num);
  } else if (num->fmt != FMT_FLT) {
    me_raise(ME_NOTIMPLEMENTED);
    return;
  } if (num->data[0] == 0) {
    me_raise(ME_DIVBYZERO);
    return;
  } if (num->data[0] == HI_BIT) {
    for(i=1; i<PRECISION; i++)
      if (num->data[i] != 0)
	break;
    if (i==PRECISION) {
      num->exponent = -num->exponent;
      return;
    }
  } ufp_inv(num->data);
  num->exponent = -num->exponent-1;
  /* val->sigfigs should remain constant
   * along with val->sign, yet it maxes out
   * at 64.
   */
  if (num->sigfigs >= BT_PRECISION)
    num->sigfigs = BT_PRECISION-1;
}

void div(num_p num, num_p dnm) {
  num_t dcpy, ncpy, dinv, tmp;
  long  seed;

  copy_num(&dcpy, dnm);
  copy_num(&dinv, dnm);
  copy_num(&ncpy, num);

  if (num_iszero(&dinv)) {
    me_raise(ME_DIVBYZERO);
    return;
  }

  if (num->fmt != FMT_FLT)
    convert_flt(num);
  if (dinv.fmt != FMT_FLT)
    convert_flt(&dinv);
  
  inv(&dinv);
  mul(num, &dinv);

  /* Deal with integer divide case */
  if ((ncpy.fmt == FMT_INT)&&(dcpy.fmt == FMT_INT)) {
    copy_num(&tmp, num);
    trunc(&tmp);
    mul(&tmp, &dcpy);

    if (cmp(&tmp, &ncpy)==0)
      trunc(num);
    else {
      quick_int(&dinv, 1);
      copy_num(&tmp, num);
      trunc(&tmp);
      add(&tmp, &dinv);
      mul(&tmp, &dcpy);
      if (cmp(&tmp, &ncpy)==0) {
	trunc(num);
	add(num, &dinv);
      } else {
	copy_num(&tmp, num);
	trunc(&tmp);
	sub(&tmp, &dinv);
	mul(&tmp, &dcpy);
	if (cmp(&tmp, &ncpy) == 0) {
	  trunc(num);
	  sub(num, &dinv);
	}
      } /* End subtract one case */
    } /* End non-simple case */
  } /* End integer divide */
}