Initial revision
This commit is contained in:
Ben Gras
266
lib/float/div_ext.c
Executable file
266
lib/float/div_ext.c
Executable file
@@ -0,0 +1,266 @@
|
||||
/*
|
||||
(c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
|
||||
See the copyright notice in the ACK home directory, in the file "Copyright".
|
||||
*/
|
||||
|
||||
/* $Header$ */
|
||||
|
||||
/*
|
||||
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;
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user