cleaned up a bit, and made add_ext and sub_ext proper entry points,

without preconditions on the arguments

mach/proto/fp/add_ext.c

@@ -14,10 +14,39 @@
 add_ext(e1,e2)
 register EXTEND	*e1,*e2;
 {
+	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);
 }

mach/proto/fp/adf4.c

@@ -27,27 +27,7 @@ _float	s1,s2;
 	}
 	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);
-	}
 	compact(&e1,(_double *)&s1,sizeof(SINGLE));
 	return(s1);
 }

mach/proto/fp/adf8.c

@@ -16,7 +16,6 @@ adf8(s2,s1)
 _double	s1,s2;
 {
 	EXTEND	e1,e2;
-	int	swap;
 
 	if (s1.__double[0] == 0 && s1.__double[1] == 0) {
 		s1 = s2;
@@ -28,28 +27,7 @@ _double	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);
-	}
 	compact(&e1,&s1,sizeof(_double));
 	return(s1);
 }

mach/proto/fp/nrm_ext.c

@@ -17,13 +17,6 @@
 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	*/
@@ -32,24 +25,25 @@ EXTEND	*e1;
 		/* 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)	{
+	if (e1->m1 == 0L)	{
-		*mant_1++ = e1->m2;
+		e1->m1 = e1->m2;
-		*mant_1-- = 0L;
+		e1->m2 = 0L;
 		e1->exp -= 32;
 	}
-	while ((*mant_1 & NORMBIT) == 0) {
+	if ((e1->m1 & NORMBIT) == 0) {
-		e1->exp--;
+		unsigned long l = ((unsigned long)NORMBIT >> 1);
-		*mant_1 <<= 1;
+		int cnt = -1;
-		if ((*mant_2 & CARRYBIT) == 0)
+
-				;	/* empty statement */
+		while (! (l & e1->m1)) {
-		else	{
+			l >>= 1;
-			*mant_1 += 1;
+			cnt--;
 		}
-		*mant_2 <<= 1;
+		e1->exp += cnt;
+		b64_sft(&(e1->m1), cnt);
 	}
 }

mach/proto/fp/sft_ext.c

@@ -8,6 +8,7 @@
 /*
 	SHIFT TWO EXTENDED NUMBERS INTO PROPER
 	ALIGNMENT FOR ADDITION (exponents are equal)
+	Numbers should not be zero on entry.
 */
 
 #include "FP_types.h"
@@ -17,7 +18,6 @@ EXTEND	*e1,*e2;
 {
 	register	EXTEND	*s;
 	register	int	diff;
-			long	tmp;
 
 	diff = e1->exp - e2->exp;
 
@@ -34,27 +34,5 @@ EXTEND	*e1,*e2;
 		s = e2;
 
 	s->exp += diff;
-
+	b64_sft(&(s->m1), 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;
-	}
 }

mach/proto/fp/shifter.c

@@ -11,26 +11,46 @@ b64_sft(e1,n)
 B64	*e1;
 int	n;
 {
+	if (n > 0) {
+		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 = (e1->l_32 >> n) | (e1->h_32 << (32 - n));
+			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 > 0) {
+		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 = (e1->h_32 << n) | (e1->l_32 >> (32 - n));
+			e1->h_32 <<= n;
+			if (e1->l_32 != 0) {
+				e1->h_32 |= (e1->l_32 >> (32 - n));
 				e1->l_32 <<= n;
 			}
+		}
+	}
 }
 
 b64_lsft(e1)

mach/proto/fp/sub_ext.c

@@ -8,21 +8,44 @@
 /*
 	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->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);
 }

mach/proto/fp/zrf_ext.c

@@ -14,13 +14,8 @@
 zrf_ext(e)
 EXTEND	*e;
 {
-	register short *ipt;
+	e->m1 = 0;
-	register int  i;
+	e->m2 = 0;
-
+	e->exp = 0;
-	/*	local CAST conversion	*/
+	e->sign = 0;
-	ipt = (short *) e;
-
-	i = sizeof(EXTEND)/sizeof(short);
-	while (i--)
-		*ipt++ = 0;
 }