123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266 |
- /*
- (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
- See the copyright notice in the ACK home directory, in the file "Copyright".
- */
- /* $Id$ */
- /*
- DIVIDE EXTENDED FORMAT
- */
- #include "FP_bias.h"
- #include "FP_trap.h"
- #include "FP_types.h"
- /*
- November 15, 1984
- This is a routine to do the work.
- There are two versions:
- One is based on the partial products method
- and makes no use possible machine instructions
- to divide (hardware dividers).
- The other is used when USE_DIVIDE is defined. It is much faster on
- machines with fast 4 byte operations.
- */
- /********************************************************/
- void
- div_ext(e1,e2)
- EXTEND *e1,*e2;
- {
- short error = 0;
- B64 result;
- register unsigned long *lp;
- #ifndef USE_DIVIDE
- short count;
- #else
- unsigned short u[9], v[5];
- register int j;
- register unsigned short *u_p = u;
- int maxv = 4;
- #endif
- if ((e2->m1 | e2->m2) == 0) {
- /*
- * Exception 8.2 - Divide by zero
- */
- trap(EFDIVZ);
- e1->m1 = e1->m2 = 0L;
- e1->exp = EXT_MAX;
- return;
- }
- if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */
- e1->exp = 0; /* make sure */
- return;
- }
- #ifndef USE_DIVIDE
- /*
- * numbers are right shifted one bit to make sure
- * that m1 is quaranteed to be larger if its
- * maximum bit is set
- */
- b64_rsft(&e1->mantissa); /* 64 bit shift right */
- b64_rsft(&e2->mantissa); /* 64 bit shift right */
- e1->exp++;
- e2->exp++;
- #endif
- /* check for underflow, divide by zero, etc */
- e1->sign ^= e2->sign;
- e1->exp -= e2->exp;
- #ifndef USE_DIVIDE
- /* do division of mantissas */
- /* uses partial product method */
- /* init control variables */
- count = 64;
- result.h_32 = 0L;
- result.l_32 = 0L;
- /* partial product division loop */
- while (count--) {
- /* first left shift result 1 bit */
- /* this is ALWAYS done */
- b64_lsft(&result);
- /* compare dividend and divisor */
- /* if dividend >= divisor add a bit */
- /* and subtract divisior from dividend */
- if ( (e1->m1 < e2->m1) ||
- ((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
- ; /* null statement */
- /* i.e., don't add or subtract */
- else {
- result.l_32++; /* ADD */
- if (e2->m2 > e1->m2)
- e1->m1 -= 1; /* carry in */
- e1->m1 -= e2->m1; /* do SUBTRACTION */
- e1->m2 -= e2->m2; /* SUBTRACTION */
- }
- /* shift dividend left one bit OR */
- /* IF it equals ZERO we can break out */
- /* of the loop, but still must shift */
- /* the quotient the remaining count bits */
- /* NB save the results of this test in error */
- /* if not zero, then the result is inexact. */
- /* this would be reported in IEEE standard */
- /* lp points to dividend */
- lp = &e1->m1;
- error = ((*lp | *(lp+1)) != 0L) ? 1 : 0;
- if (error) { /* more work */
- /* assume max bit == 0 (see above) */
- b64_lsft(&e1->mantissa);
- continue;
- }
- else
- break; /* leave loop */
- } /* end of divide by subtraction loop */
- if (count > 0) {
- lp = &result.h_32;
- if (count > 31) { /* move to higher word */
- *lp = *(lp+1);
- count -= 32;
- *(lp+1) = 0L; /* clear low word */
- }
- if (*lp)
- *lp <<= count; /* shift rest of way */
- lp++; /* == &result.l_32 */
- if (*lp) {
- result.h_32 |= (*lp >> 32-count);
- *lp <<= count;
- }
- }
- #else /* USE_DIVIDE */
- u[4] = (e1->m2 & 1) << 15;
- b64_rsft(&(e1->mantissa));
- u[0] = e1->m1 >> 16;
- u[1] = e1->m1;
- u[2] = e1->m2 >> 16;
- u[3] = e1->m2;
- u[5] = 0; u[6] = 0; u[7] = 0;
- v[1] = e2->m1 >> 16;
- v[2] = e2->m1;
- v[3] = e2->m2 >> 16;
- v[4] = e2->m2;
- while (! v[maxv]) maxv--;
- result.h_32 = 0;
- result.l_32 = 0;
- lp = &result.h_32;
- /*
- * Use an algorithm of Knuth (The art of programming, Seminumerical
- * algorithms), to divide u by v. u and v are both seen as numbers
- * with base 65536.
- */
- for (j = 0; j <= 3; j++, u_p++) {
- unsigned long q_est, temp;
- if (j == 2) lp++;
- if (u_p[0] == 0 && u_p[1] < v[1]) continue;
- temp = ((unsigned long)u_p[0] << 16) + u_p[1];
- if (u_p[0] >= v[1]) {
- q_est = 0x0000FFFFL;
- }
- else {
- q_est = temp / v[1];
- }
- temp -= q_est * v[1];
- while (temp < 0x10000 && v[2]*q_est > ((temp<<16)+u_p[2])) {
- q_est--;
- temp += v[1];
- }
- /* Now, according to Knuth, we have an estimate of the
- quotient, that is either correct or one too big, but
- almost always correct.
- */
- if (q_est != 0) {
- int i;
- unsigned long k = 0;
- int borrow = 0;
- for (i = maxv; i > 0; i--) {
- unsigned long tmp = q_est * v[i] + k + borrow;
- unsigned short md = tmp;
- borrow = (md > u_p[i]);
- u_p[i] -= md;
- k = tmp >> 16;
- }
- k += borrow;
- borrow = u_p[0] < k;
- u_p[0] -= k;
- if (borrow) {
- /* So, this does not happen often; the estimate
- was one too big; correct this
- */
- *lp |= (j & 1) ? (q_est - 1) : ((q_est-1)<<16);
- borrow = 0;
- for (i = maxv; i > 0; i--) {
- unsigned long tmp
- = v[i]+(unsigned long)u_p[i]+borrow;
-
- u_p[i] = tmp;
- borrow = tmp >> 16;
- }
- u_p[0] += borrow;
- }
- else *lp |= (j & 1) ? q_est : (q_est<<16);
- }
- }
- #ifdef EXCEPTION_INEXACT
- u_p = &u[0];
- for (j = 7; j >= 0; j--) {
- if (*u_p++) {
- error = 1;
- break;
- }
- }
- #endif
- #endif
- #ifdef EXCEPTION_INEXACT
- if (error) {
- /*
- * report here exception 8.5 - Inexact
- * from Draft 8.0 of IEEE P754:
- * In the absence of an invalid operation exception,
- * if the rounded result of an operation is not exact or if
- * it overflows without a trap, then the inexact exception
- * shall be assigned. The rounded or overflowed result
- * shall be delivered to the destination.
- */
- INEXACT();
- #endif
- e1->mantissa = result;
- nrm_ext(e1);
- if (e1->exp < EXT_MIN) {
- /*
- * Exception 8.4 - Underflow
- */
- trap(EFUNFL); /* underflow */
- e1->exp = EXT_MIN;
- e1->m1 = e1->m2 = 0L;
- return;
- }
- if (e1->exp >= EXT_MAX) {
- /*
- * Exception 8.3 - Overflow
- */
- trap(EFOVFL); /* overflow */
- e1->exp = EXT_MAX;
- e1->m1 = e1->m2 = 0L;
- return;
- }
- }
|