path: root/src/benchmarks/cfrac/pdivmod.c

#include "pdefs.h"
#include "precision.h"

#ifdef DEBUG
#include <stdio.h>
#endif

#ifdef ASM_16BIT
#include "asm16bit.h"
#endif

/*
 *   Divide u (dividend) by v (divisor); If non-null, qp and rp are set to
 *   quotient and remainder.  The result returned will be *qp, unless qp is
 *   NULL, then *rp will be returned if non-null, otherwise pUndef is returned.
 *
 *  Produce:
 * 
 *	 q (quotient)	= u div v	(v != 0)
 *			  truncation is toward zero
 *	
 *	 r (remainder)	= u mod v
 *			= u - u div v * v      (v != 0)
 *			= u		       (v == 0)
 *			  ( e.g.  u == q*v + r )
 *			remainder has same sign and dividend
 *
 * Note: this has opposite convention than the C standard div fuction,
 *	 but the same convention of the typical C "/" operator
 *	 It is also inconvienient for the mod function.	 
 */
/*
 *	This algorithm is taken almost verbatum from Knuth Vol 2.
 *	Please note the following trivial(?) array index
 *	transformations (since MSD to LSD order is reversed):
 *	
 *	q[0..m] to Q[0..m]   thus   q[i] == Q[m-i]
 *	r[1..n]	   R[0..n-1]	    r[i] == R[n+1-i]
 *	u[0..m+n]  w[0..m+n]	    u[i] == w[m+n-i]
 *	v[1..n]	   x[0..n-1]	    v[i] == x[n-i]
 *	
 *	let N == n - 1 so that n == N + 1 thus:
 *	
 *	q[0..m] to Q[0..m]   thus   q[i] == Q[m-i]
 *	r[1..n]	   R[0..N]	    r[i] == R[N+2-i]
 *	u[0..m+n]  w[0..m+N+1]	    u[i] == w[m+N+1-i]
 *	v[1..n]	   x[0..N]	    v[i] == x[N+1-i]
 */

/*
 * Note: Be very observent of the usage of uPtr, and vPtr.
 *       They are used to point to u, v, w, q or r as necessary.
 */
precision pdivmod(u, v, qp, rp)
   precision u, v, *qp, *rp;
{
   register digitPtr	uPtr, vPtr, qPtr, LastPtr;

   register accumulator temp;		       /* 0 <= temp < base^2 */
   register digit	carry;			   /* 0 <= carry < 2 */
   register digit	hi;			   /* 0 <= hi < base */

   register posit	n, m;
   digit		d;			    /* 0 <= d < base */
   digit		qd;			   /* 0 <= qd < base */
#ifdef DEBUG
   int			i;
#endif

   precision		q, r, w;   /* quotient, remainder, temporary */

   n = v->size;					  /* size of v and r */

   (void) pparm(u);
   (void) pparm(v);
   if (u->size < n) {
      q = pUndef;
      r = pUndef;
      pset(&q, pzero);
      pset(&r, u);
      goto done;
   }

   m = u->size - n;

   uPtr = u->value + m + n;
   vPtr = v->value + n;

   q = palloc(m + 1);
   if (q == pUndef) return q;

   q->sign = (u->sign != v->sign);		  /* can generate -0 */

   r = palloc(n);
   if (r == pUndef) {
      pdestroy(q);
      return r;
   }
   r->sign = u->sign;	    
/* 
 * watch out! does this function return: q=floor(a/b) or trunc(a/b)?
 * it's currently the latter, but every mathmaticion I have talked to
 * prefers the former so that a % b returns between 0 to b-1.  The
 * problem is that this is slower and disagrees with C common practice.
 */
   qPtr = q->value + m + 1;

   if (n == 1) {
      d = *--vPtr;			     /* d is only digit of v */
      if (d == 0) {				  /* divide by zero? */
	 q = pnew(errorp(PDOMAIN, "pdivmod", "divide by zero"));
      } else {				      /* single digit divide */
#ifndef ASM_16BIT
	 vPtr = r->value + n;
	 hi = 0;			  /* hi is current remainder */
	 do {
	     temp    = mulBase(hi);	/* 0 <= temp <= (base-1)^2   */
	     temp   += *--uPtr;		/* 0 <= temp <= base(base-1) */
	     hi	     = uModDiv(temp, d, --qPtr); /* 0 <= hi < base */
	 } while (uPtr > u->value);
	 *--vPtr = hi;
#else
	qPtr -= m + 1;
	*(r->value) = memdivw1(qPtr, u->value, m + 1, d);
#endif
      }					  
   } else {					/* muti digit divide */
      /*
       * normalize:   multiply u and v by d so hi digit of v > b/2 
       */
      d = BASE / (*--vPtr+1);			  /* high digit of v */

      w = palloc(n);					/* size of v */
      if (w == pUndef) return w;

#ifndef ASM_16BIT
      vPtr = v->value;
      uPtr = w->value;		    /* very confusing.	just a temp */
      LastPtr = vPtr + n;
      hi = 0;
      do {				    /* single digit multiply */
	 temp	 = uMul(*vPtr++, d);   /* 0<= temp <= base(base-1)/2 */
	 temp	+= hi;			/* 0 <= temp <= (base^2-1)/2 */
	 hi	 = divBase(temp);	     /* 0 <= hi	  < base / 2 */
	 *uPtr++ = modBase(temp);	     /* 0 <= hi	  < base / 2 */
      } while (vPtr < LastPtr);			  /* on exit hi == 0 */
#else
      hi = memmulw1(w->value, v->value, n, d);
#endif

      pset(&v, w);
      pdestroy(w);

      w = palloc(m + n + 1);
      if (w == pUndef) return w;

#ifndef ASM_16BIT
      uPtr = u->value;
      vPtr = w->value;		     /* very confusing.	 just a temp */
      LastPtr = uPtr + m + n;
      do {				    /* single digit multiply */
	 temp	 = uMul(*uPtr++, d);	      
	 temp	+= hi;	      
	 hi	 = divBase(temp);
	 *vPtr++ = modBase(temp);
      } while (uPtr < LastPtr);
      *vPtr = hi;				 /* note extra digit */
#else
      hi = memmulw1(w->value, u->value, m + n, d);
      w->value[m + n] = hi;
#endif

      pset(&u, w);
      pdestroy(w);

#ifdef DEBUG
      printf("m = %d    n = %d\nd = %d\n", m, n, d);
      printf("norm u = "); pshow(u);
      printf("norm v = "); pshow(v);
#endif

      uPtr = u->value + m + 1;		  /* current least significant digit */
      do {
	 --uPtr;
#ifdef DEBUG
	 printf("   u = ");
	 for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
	 putchar('\n');
	 printf("   v = "); 
	 for (i = 1; i < 3; i++) printf("%.*x ", sizeof(digit) * 2, 
	    v->value[n-i]);
	 putchar('\n');
#endif
#ifndef ASM_16BIT
	 vPtr    = v->value + n;
	 LastPtr = uPtr + n;
	 if (*LastPtr == *--vPtr) {			 /* guess next digit */
	    qd = BASE - 1;
	 } else {
	    temp  = mulBase(*LastPtr);
	    temp += *--LastPtr;				/* 0 <= temp< base^2 */
	    temp  = uModDiv(temp, *vPtr, &qd);
	    --vPtr;
	    --LastPtr;
	    while (uMul(*vPtr, qd) > mulBase(temp) + *LastPtr) {
	       --qd;
	       temp += vPtr[1];
	       if (temp >= BASE) break;  /* if so, vPtr*qd <= temp*base */
	    }		
	    LastPtr += 2;
	 }
	 /*
	  * Single digit Multiply then Subtract
	  */
	 vPtr  = v->value;
	 carry = 1;			 /* noborrow bit */
	 hi    = 0;			 /* hi digit of multiply */
	 do {
	    /* multiply */
	    temp    = uMul(qd, *vPtr++); /* 0 <= temp <= (base-1)^2   */
	    temp   += hi;		 /* 0 <= temp <= base(base-1) */
	    hi	    = divBase(temp);
	    temp    = modBase(temp);
	    /* subtract */
	    temp    = (BASE-1) - temp;	 /* 0 <= temp < base */
	    temp   += *uPtr + carry;	 /* 0 <= temp < 2*base */
	    carry   = divBase(temp);
	    *uPtr++ = modBase(temp);	 /* 0 <= carry < 2 */
	 } while (uPtr < LastPtr);
	 temp  = (BASE-1) - hi;
	 temp += *uPtr + carry;
	 carry = divBase(temp);
	 *uPtr = modBase(temp);
	 uPtr -= n;
#else
#if 0
	 carry = !memmulsubw(uPtr, v->value, n, qd); 	/* 1 if noborrow */
#endif
	 carry = !memdivw(uPtr, v->value, n, &qd); 	/* 1 if noborrow */
#endif
#ifdef DEBUG
	 printf("   qhat = %.*x\n", sizeof(digit) * 2, qd);
	 printf("   new u = "); 
	 for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
	 putchar('\n');
#endif
	 if (carry == 0) {		 /* Test remainder, add back */
	    vPtr    = v->value;
	    LastPtr = uPtr + n;
	    do {
	       temp    = *uPtr + *vPtr++;
	       temp   += carry;
	       carry   = divBase(temp);
	       *uPtr++ = modBase(temp);
	    } while (uPtr < LastPtr);
	    *uPtr += carry - BASE;	   /* real strange but works */
	    uPtr -= n;
	    --qd;
#ifdef DEBUG
	    printf("   decrementing q...adding back\n");
	    printf("   fixed u = "); 
	    for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
	    putchar('\n');
	    printf("   newq = %.*x\n", sizeof(digit) * 2, qd);
#endif
	 }
	 *--qPtr = qd;			/* one leading zero possible */
#ifdef DEBUG
	 putchar('\n'); 
#endif
      } while (uPtr > u->value);

      /*
       * Un-normalize to get remainder
       */
#ifndef ASM_16BIT
      uPtr = u->value + n;		/* skip hi digit (it's zero) */
      vPtr = r->value + n;
      hi   = 0;				/* hi is current remainder   */
      do {				/* single digit divide	     */
	  temp	  = mulBase(hi);	/* 0<=temp < base^2-(base-1) */
	  temp	 += *--uPtr;		/* 0 <= temp < base^2	     */
	  hi	  = uModDiv(temp, d, --vPtr);
      } while (uPtr > u->value);	/* carry will be zero	     */
#else 
      carry = memdivw1(r->value, u->value, n, d); 	 /* always 0 */
#endif
      pnorm(r);		      /* remainder may have many leading 0's */
   }

   if (m > 0 && qPtr[m] == 0) {
      --(q->size);					 /* normalize */
   }
   if (q->size == 1 && *qPtr == 0) q->sign = false;

done:

   pdestroy(u);
   pdestroy(v);

   if (rp == (precision *) -1) {
      if (qp != pNull) pset(qp, q);
      pdestroy(q);
      return presult(r);
   } else if (qp == (precision *) -1) {
      if (rp != pNull) pset(rp, r);
      pdestroy(r);
      return presult(q);
   }
   if (qp != pNull) pset(qp, q);
   if (rp != pNull) pset(rp, r);
   pdestroy(q);
   pdestroy(r);
   return pUndef;
}