/* new double-float divide routine, by jrd */

#include "flonum.h"

#ifdef L_divdf3
double __divdf3(num, denom)
double num, denom;
{
  double local_num = num;
  double local_denom = denom;
  struct bitdouble * num_bdp = (struct bitdouble *)(&local_num);
  struct bitdouble * den_bdp = (struct bitdouble *)(&local_denom);
  short num_exp = num_bdp->exp,
	den_exp = den_bdp->exp;
  short result_sign = 0;
/*  register */ unsigned long num_m1, num_m2, num_m3, num_m4;
  register unsigned long den_m1, den_m2, den_m3, den_m4;
  unsigned long result_mant[3];
  unsigned long result_mask;
  short result_idx;

  if ((num_exp == 0) || (den_exp == 0))		/* zzz should really cause trap */
  	return(0.0);

  /* deal with signs */
  result_sign = result_sign + num_bdp->sign - den_bdp->sign;
	
  /* unpack the numbers */
  num_m1 = num_bdp->mant1 | 0x00100000;		/* hidden bit */
  num_m2 = num_bdp->mant2;
  num_m4 = num_m3 = 0;
  den_m1 = den_bdp->mant1 | 0x00100000;		/* hidden bit */
  den_m2 = den_bdp->mant2;
  den_m4 = den_m3 = 0;
  
#if 0
/* buy a little extra accuracy by shifting num and denum up 10 bits */
  for (result_idx /* loop counter */ = 0 ; result_idx < 10 ; result_idx++)
	{
	ASL3(num_m1, num_m2, num_m3);
	ASL3(den_m1, den_m2, den_m3);
	}
#endif

/* hot wire result mask and index */
  result_mask = 0x00100000;			/* start at top mantissa bit */
  result_idx = 0;				/* of first word */
  result_mant[0] = 0;
  result_mant[1] = 0;
  result_mant[2] = 0;

/* if denom is greater than num here, shift denom right one and dec num expt */
  if (den_m1 < num_m1)
  	goto kludge1;				/* C is assembly language,
						   remember? */
  if (den_m1 > num_m1)
  	goto kludge0;
  if (den_m2 <= num_m2)				/* first word eq, try 2nd */
  	goto kludge1;

kludge0:

  num_exp--;
  ASR4(den_m1, den_m2, den_m3, den_m4);

kludge1:
  
  for ( ; !((result_idx == 2) && (result_mask & 0x40000000)) ; )
	{
/* if num >= den, subtract den from num and set bit in result */
	if (num_m1 > den_m1) goto kludge2;
	if (num_m1 < den_m1) goto kludge3;
	if (num_m2 > den_m2) goto kludge2;
	if (num_m2 < den_m2) goto kludge3;
	if (num_m3 > den_m3) goto kludge2;
	if (num_m3 < den_m3) goto kludge3;
	if (num_m4 < den_m4) goto kludge3;

kludge2:
	result_mant[result_idx] |= result_mask;
	SUB4(den_m1, den_m2, den_m3, den_m4, num_m1, num_m2, num_m3, num_m4);
kludge3:
	ASR4(den_m1, den_m2, den_m3, den_m4);
	result_mask >>= 1;
	if (result_mask == 0)		/* next word? */
		{
		result_mask = 0x80000000;
		result_idx++;
		}
	}

/* stick in last half-bit if necessary */
  if (result_mant[2])
	{
	den_m1 = 0;		/* handy register */
	den_m2 = 1;
/* busted
	ADD2(den_m1, den_m2, result_mant[0], result_mant[1]);
*/
	result_mant[1]++;
	if (result_mant[1] == 0)
		result_mant[0]++;

	if (result_mant[0] & 0x00200000)	/* overflow? */
		{
		num_exp--;
		}
	}
  /* compute the resultant exponent */
  num_exp = num_exp - den_exp + 0x3FF;

  /* reconstruct the result in local_num */
  num_bdp->sign = result_sign;
  num_bdp->exp = num_exp;
  num_bdp->mant1 = result_mant[0] & 0xFFFFF;
  num_bdp->mant2 = result_mant[1];
  
  /* done! */
  return(local_num);
}
#endif