add_ext(e1,e2)
register EXTEND *e1,*e2;
{
- if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
- b64_rsft(&e1->m1); /* shift mantissa one bit RIGHT */
- e1->m1 |= 0x80000000L; /* set max bit */
- e1->exp++; /* increase the exponent */
+ if ((e2->m1 | e2->m2) == 0L) {
+ return;
+ }
+ if ((e1->m1 | e1->m2) == 0L) {
+ *e1 = *e2;
+ return;
+ }
+ sft_ext(e1, e2); /* adjust mantissas to equal powers */
+ if (e1->sign != e2->sign) {
+ /* e1 + e2 = e1 - (-e2) */
+ if (e2->m1 > e1->m1 ||
+ (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+ /* abs(e2) > abs(e1) */
+ if (e1->m2 > e2->m2) {
+ e2->m1 -= 1; /* carry in */
+ }
+ e2->m1 -= e1->m1;
+ e2->m2 -= e1->m2;
+ *e1 = *e2;
+ }
+ else {
+ if (e2->m2 > e1->m2)
+ e1->m1 -= 1; /* carry in */
+ e1->m1 -= e2->m1;
+ e1->m2 -= e2->m2;
+ }
+ }
+ else {
+ if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
+ b64_rsft(&e1->m1); /* shift mantissa one bit RIGHT */
+ e1->m1 |= 0x80000000L; /* set max bit */
+ e1->exp++; /* increase the exponent */
+ }
}
nrm_ext(e1);
}
}
extend((_double *)&s1,&e1,sizeof(SINGLE));
extend((_double *)&s2,&e2,sizeof(SINGLE));
- /* if signs differ do subtraction */
- /* result in e1 */
- if (e1.sign ^ e2.sign) {
- /* set sign of e1 to sign of largest */
- swap = (e2.exp > e1.exp) ? 1 : (e2.exp < e1.exp) ? 0 :
- (e2.m1 > e1.m1 ) ? 1 : 0;
- /* adjust mantissas to equal powers */
- sft_ext(&e1,&e2);
- /* subtract the extended formats */
- if (swap) {
- sub_ext(&e2,&e1);
- e1 = e2;
- }
- else
- sub_ext(&e1,&e2);
- }
- else {
- /* adjust mantissas to equal powers */
- sft_ext(&e1,&e2);
- add_ext(&e1,&e2);
- }
+ add_ext(&e1,&e2);
compact(&e1,(_double *)&s1,sizeof(SINGLE));
return(s1);
}
_double s1,s2;
{
EXTEND e1,e2;
- int swap;
if (s1.__double[0] == 0 && s1.__double[1] == 0) {
s1 = s2;
extend(&s1,&e1,sizeof(_double));
extend(&s2,&e2,sizeof(_double));
- /* adjust mantissas to equal powers */
- if (e1.sign ^ e2.sign) { /* signs are different */
- /* determine which is largest number */
- swap = (e2.exp > e1.exp) ? 1 : (e2.exp < e1.exp) ? 0 :
- (e2.m1 > e1.m1 ) ? 1 : (e2.m1 < e1.m1 ) ? 0 :
- (e2.m2 > e1.m2 ) ? 1 : 0;
- /* adjust mantissas to equal powers */
- sft_ext(&e1,&e2);
- /* subtract the extended formats */
- if (swap) { /* &e2 is the largest number */
- sub_ext(&e2,&e1);
- e1 = e2;
- }
- else {
- sub_ext(&e1,&e2);
- }
- }
- else {
- /* adjust mantissas to equal powers */
- sft_ext(&e1,&e2);
- add_ext(&e1,&e2);
- }
+ add_ext(&e1,&e2);
compact(&e1,&s1,sizeof(_double));
return(s1);
}
nrm_ext(e1)
EXTEND *e1;
{
- register unsigned long *mant_1;
- register unsigned long *mant_2;
-
- /* local CAST conversion */
- mant_1 = (unsigned long *) &e1->m1;
- mant_2 = (unsigned long *) &e1->m2;
-
/* we assume that the mantissa != 0 */
/* if it is then just return */
/* to let it be a problem elsewhere */
/* infinite loop is generated when */
/* mantissa is zero */
- if ((*mant_1 | *mant_2) == 0L)
+ if ((e1->m1 | e1->m2) == 0L)
return;
/* if top word is zero mov low word */
/* to top word, adjust exponent value */
- if (*mant_1 == 0L) {
- *mant_1++ = e1->m2;
- *mant_1-- = 0L;
+ if (e1->m1 == 0L) {
+ e1->m1 = e1->m2;
+ e1->m2 = 0L;
e1->exp -= 32;
}
- while ((*mant_1 & NORMBIT) == 0) {
- e1->exp--;
- *mant_1 <<= 1;
- if ((*mant_2 & CARRYBIT) == 0)
- ; /* empty statement */
- else {
- *mant_1 += 1;
+ if ((e1->m1 & NORMBIT) == 0) {
+ unsigned long l = ((unsigned long)NORMBIT >> 1);
+ int cnt = -1;
+
+ while (! (l & e1->m1)) {
+ l >>= 1;
+ cnt--;
}
- *mant_2 <<= 1;
+ e1->exp += cnt;
+ b64_sft(&(e1->m1), cnt);
}
}
/*
SHIFT TWO EXTENDED NUMBERS INTO PROPER
ALIGNMENT FOR ADDITION (exponents are equal)
+ Numbers should not be zero on entry.
*/
#include "FP_types.h"
{
register EXTEND *s;
register int diff;
- long tmp;
diff = e1->exp - e2->exp;
s = e2;
s->exp += diff;
-
- if (diff > 63) { /* no relative value */
- s->m1 = 0L;
- s->m2 = 0L;
- return;
- }
-
- if (diff > 32) {
- diff -= 32;
- s->m2 = s->m1;
- s->m1 = 0L;
- }
- if (diff) {
- if (s->m1) {
- tmp = s->m1;
- tmp <<= (32-diff);
- s->m1 >>= diff;
- }
- else
- tmp = 0L;
- s->m2 >>= diff;
- s->m2 |= tmp;
- }
+ b64_sft(&(s->m1), diff);
}
B64 *e1;
int n;
{
- if (n >= 32) {
- e1->l_32 = e1->h_32;
- e1->h_32 = 0;
- n -= 32;
- }
if (n > 0) {
- e1->l_32 = (e1->l_32 >> n) | (e1->h_32 << (32 - n));
- e1->h_32 >>= n;
+ if (n > 63) {
+ e1->l_32 = 0;
+ e1->h_32 = 0;
+ return;
+ }
+ if (n >= 32) {
+ e1->l_32 = e1->h_32;
+ e1->h_32 = 0;
+ n -= 32;
+ }
+ if (n > 0) {
+ e1->l_32 >>= n;
+ if (e1->h_32 != 0) {
+ e1->l_32 |= (e1->h_32 << (32 - n));
+ e1->h_32 >>= n;
+ }
+ }
return;
}
n = -n;
- if (n >= 32) {
- e1->h_32 = e1->l_32;
- e1->l_32 = 0;
- n -= 32;
- }
if (n > 0) {
- e1->h_32 = (e1->h_32 << n) | (e1->l_32 >> (32 - n));
- e1->l_32 <<= n;
+ if (n > 63) {
+ e1->l_32 = 0;
+ e1->h_32 = 0;
+ return;
+ }
+ if (n >= 32) {
+ e1->h_32 = e1->l_32;
+ e1->l_32 = 0;
+ n -= 32;
+ }
+ if (n > 0) {
+ e1->h_32 <<= n;
+ if (e1->l_32 != 0) {
+ e1->h_32 |= (e1->l_32 >> (32 - n));
+ e1->l_32 <<= n;
+ }
+ }
}
}
/*
SUBTRACT 2 EXTENDED FORMAT NUMBERS
*/
- /*
- * adf (addition routines) use this rather than
- * add_ext when the signs of the numbers are different.
- * sub_ext requires that e1 >= e2 on entry
- * otherwise nonsense results. If you use this routine
- * make certain this requirement is met.
- */
#include "FP_types.h"
sub_ext(e1,e2)
EXTEND *e1,*e2;
{
- if (e2->m2 > e1->m2)
- e1->m1 -= 1; /* carry in */
- e1->m1 -= e2->m1;
- e1->m2 -= e2->m2;
+ if ((e2->m1 | e2->m2) == 0L) {
+ return;
+ }
+ if ((e1->m1 | e1->m2) == 0L) {
+ *e1 = *e2;
+ e1->sign = e2->sign ? 0 : 1;
+ return;
+ }
+ sft_ext(e1, e2);
+ if (e1->sign != e2->sign) {
+ /* e1 - e2 = e1 + (-e2) */
+ if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
+ b64_rsft(&e1->m1); /* shift mantissa one bit RIGHT */
+ e1->m1 |= 0x80000000L; /* set max bit */
+ e1->exp++; /* increase the exponent */
+ }
+ }
+ else if (e2->m1 > e1->m1 ||
+ (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+ /* abs(e2) > abs(e1) */
+ if (e1->m2 > e2->m2) {
+ e2->m1 -= 1; /* carry in */
+ }
+ e2->m1 -= e1->m1;
+ e2->m2 -= e1->m2;
+ *e1 = *e2;
+ e1->sign = e2->sign ? 0 : 1;
+ }
+ else {
+ if (e2->m2 > e1->m2)
+ e1->m1 -= 1; /* carry in */
+ e1->m1 -= e2->m1;
+ e1->m2 -= e2->m2;
+ }
nrm_ext(e1);
}
zrf_ext(e)
EXTEND *e;
{
- register short *ipt;
- register int i;
-
- /* local CAST conversion */
- ipt = (short *) e;
-
- i = sizeof(EXTEND)/sizeof(short);
- while (i--)
- *ipt++ = 0;
+ e->m1 = 0;
+ e->m2 = 0;
+ e->exp = 0;
+ e->sign = 0;
}
public / ack.git / commitdiff