Minix2.0/src/lib/float/div_ext.c
/*
(c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
See the copyright notice in the ACK home directory, in the file "Copyright".
*/
/* $Header: div_ext.c,v 1.10 93/01/05 12:04:47 ceriel Exp $ */
/*
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;
}
}