ngf4.c
ngf8.c
nrm_ext.c
-prt_dbl.c
-prt_ext.c
sbf4.c
sbf8.c
sft_ext.c
/********************************************************/
/*
+ Type definitions for C Floating Point Package
include file for floating point package
*/
/********************************************************/
EXTEND: double precision extended format
*/
/********************************************************/
+
+#ifdef EXT_DEBUG
+#ifndef __FPSTDIO
+#define __FPSTDIO
+#include <stdio.h>
+#endif
+#endif
+
+#ifndef __FPTYPES
+#define __FPTYPES
typedef unsigned long _float;
typedef union {
unsigned long m1;
unsigned long m2; /* includes guard byte */
} EXTEND;
-#ifdef PRT_EXT
-#include <stdio.h>
#endif
extend.$(SUF) compact.$(SUF)\
add_ext.$(SUF) div_ext.$(SUF) mul_ext.$(SUF) nrm_ext.$(SUF)\
sft_ext.$(SUF) sub_ext.$(SUF) zrf_ext.$(SUF)\
- adder.$(SUF) shifter.$(SUF)\
- prt_dbl.$(SUF)\
- fptrp.$(SUF) prt_ext.$(SUF) # debugging
+ adder.$(SUF) shifter.$(SUF) fptrp.$(SUF)
SLIST = cff4.s cff8.s cfu.s cmf4.s cmf8.s\
cuf4.s cuf8.s\
dvf4.s dvf8.s fef4.s fef8.s\
extend.s compact.s\
add_ext.s div_ext.s mul_ext.s nrm_ext.s\
sft_ext.s sub_ext.s zrf_ext.s\
- adder.s shifter.s\
- prt_dbl.s\
- fptrp.s prt_ext.s # debugging
+ adder.s shifter.s fptrp.s
SRC = FP_bias.h FP_shift.h FP_trap.h FP_types.h adder.h get_put.h\
cff4.c cff8.c cfu.c cmf4.c cmf8.c\
extend.c compact.c\
add_ext.c div_ext.c mul_ext.c nrm_ext.c\
sft_ext.c sub_ext.c zrf_ext.c\
- adder.c shifter.c\
- prt_dbl.c\
- fptrp.e prt_ext.c
+ adder.c shifter.c fptrp.e
all: FP_$(MACH).a
ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/shifter.c
ed - shifter.s <$(CDIR)/FP.script
ack -c -m$(MACH) $(EMFLAGS) shifter.s
-
-prt_dbl.$(SUF): $(CDIR)/prt_dbl.c
- ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/prt_dbl.c
- ed - prt_dbl.s <$(CDIR)/FP.script
- ack -c -m$(MACH) $(EMFLAGS) prt_dbl.s
-
-prt_ext.$(SUF): $(CDIR)/prt_ext.c
- ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/prt_ext.c
- ed - prt_ext.s <$(CDIR)/FP.script
- ack -c -m$(MACH) $(EMFLAGS) prt_ext.s
/* $Header$ */
/*
-#define PRT_EXT
ADD TWO EXTENDED FORMAT NUMBERS
*/
add_ext(e1,e2)
register EXTEND *e1,*e2;
{
-#ifdef PRT_EXT
- prt_ext("before ADD #1:",e1);
- prt_ext("before ADD #2:",e2);
-#endif PRT_EXT
if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
b64_sft(&e1->m1,1); /* shift mantissa one bit RIGHT */
e1->m1 |= 0x80000000L; /* set max bit */
e1->exp++; /* increase the exponent */
}
-#ifdef PRT_EXT
- prt_ext("AFTER ADD :",e1);
-#endif
nrm_ext(e1);
-#ifdef PRT_EXT
- prt_ext("AFTER NRM :",e1);
-#endif
}
* these are the routines the routines to do 32 and 64-bit addition
*/
-# ifdef DEBUG
+# ifdef EXT_DEBUG
# include <stdio.h>
# endif
*/
register B64 *e1,*e2;
{
- register short overflow;
- short carry;
+ register int overflow;
+ int carry;
/* add higher pair of 32 bits */
overflow = b32_add(&e1->h_32,&e2->h_32);
/* add lower pair of 32 bits */
carry = b32_add(&e1->l_32,&e2->l_32);
-# ifdef DEBUG
+# ifdef EXT_DEBUG
printf("\t\t\t\t\tb64_add: overflow (%d); internal carry(%d)\n",
overflow,carry);
fflush(stdout);
b32_add(e1,e2)
register unsigned long *e1,*e2;
{
- register short carry;
+ register int carry;
if (*e1 & *e2 & MAXBIT) /* both max_bits are set */
carry = TRUE; /* so there is a carry */
? UNKNOWN
/* both are clear - no carry */
: FALSE;
-# ifdef DEBUG
+# ifdef EXT_DEBUG
fflush(stdout);
printf("\t\t\t\t\tb32_add: overflow before add(%d) test(%d)\n",
carry,(*e1&MAXBIT)?FALSE:TRUE);
# endif
*e1 += *e2;
-# ifdef DEBUG
+# ifdef EXT_DEBUG
printf("%08X\n",*e1);
fflush(stdout);
# endif
/* $Header$ */
/*
- ADD TWO FLOATS - SINGLE
+ ADD TWO FLOATS - SINGLE (ADF 4)
*/
#include "FP_types.h"
/* if signs differ do subtraction */
/* result in e1 */
if (e1.sign ^ e2.sign) {
-#ifdef DEBUG
- printf("\t\t\tADF calls SUBTRACT\n");
-#endif DEBUG
/* 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;
sft_ext(&e1,&e2);
/* subtract the extended formats */
if (swap) {
-#ifdef DEBUG
- printf("\t\t\t\tSWAP\n");
-#endif DEBUG
sub_ext(&e2,&e1);
e1 = e2;
}
sub_ext(&e1,&e2);
}
else {
-#ifdef DEBUG
- printf("\t\t\tADF calls ADDITION\n");
-#endif DEBUG
/* adjust mantissas to equal powers */
sft_ext(&e1,&e2);
add_ext(&e1,&e2);
/* $Header$ */
/*
- ADD TWO FLOATS - DOUBLE
+ ADD TWO FLOATS - DOUBLE (ADF 8)
*/
#include "FP_types.h"
_double s1,s2;
{
EXTEND e1,e2;
- short swap;
+ int swap;
if (s1.__double[0] == 0 && s1.__double[1] == 0) {
s1 = s2;
extend(&s1,&e1,sizeof(_double));
extend(&s2,&e2,sizeof(_double));
-#ifdef PRT_EXT
- prt_ext("ADF8: e1",&e1);
- prt_ext("ADF8: e2",&e2);
-#endif
/* adjust mantissas to equal powers */
if (e1.sign ^ e2.sign) { /* signs are different */
/* determine which is largest number */
sft_ext(&e1,&e2);
/* subtract the extended formats */
if (swap) { /* &e2 is the largest number */
-#ifdef PRT_EXT
- fprintf(stderr,"ADF8: swaps and subtracts extended\n");
-#endif
sub_ext(&e2,&e1);
e1 = e2;
}
sft_ext(&e1,&e2);
add_ext(&e1,&e2);
}
-#ifdef PRT_EXT
- prt_ext("ADF8: e1 result",&e1);
-#endif
compact(&e1,&s1,sizeof(_double));
return(s1);
}
/* $Header$ */
/*
- CONVERT DOUBLE TO FLOAT
+ CONVERT DOUBLE TO FLOAT (CFF 8 4)
This routine works quite simply. A floating point
of size 08 is converted to extended format.
EXTEND buf;
extend(&src,&buf,8); /* no matter what */
-#ifdef PRT_EXT
- prt_ext("CFF4() entry:",&buf);
- fprintf(stderr,"ds(%d),ss(%d),src(%08X%08X)\n",8,4,src.__double[0],
- src.__double[1]);
-#endif PRT_EXT
compact(&buf,(_double *) &(src.__double[1]),4);
-#ifdef PRT_EXT
- fprintf(stderr,"CFF4() exit : %08X\n",src.__double[1]);
-#endif PRT_EXT
}
/* $Header$ */
/*
- CONVERT FLOAT TO DOUBLE
+ CONVERT FLOAT TO DOUBLE (CFF 4 8)
This routine works quite simply. A floating point
of size 04 is converted to extended format.
EXTEND buf;
extend((_double *) &src,&buf,4); /* no matter what */
-#ifdef PRT_EXT
- prt_ext("CFF8() entry:",&buf);
- fprintf(stderr,"ds(%d),ss(%d),src(%08X)\n",4,8,src);
-#endif
compact(&buf,&src,8);
-#ifdef DEBUG
- fprintf(stderr,"CFF8() exit : %08X",src.__double[0]);
- fprintf(stderr,"%08X\n",src.__double[1]);
-#endif DEBUG
}
/* $Header$ */
/*
- CONVERT FLOAT TO UNSIGNED
+ CONVERT FLOAT TO SIGNED (CFI m n)
N.B. The caller must know what it is getting.
A LONG is always returned. If it is an
short newint, max_exp;
extend(&src,&buf,ss); /* get extended format */
-#ifdef PRT_EXT
- prt_ext("CFI() entry:",&buf);
-#endif PRT_EXT
buf.exp--; /* additional bias correction */
if (buf.exp < 1) { /* no conversion needed */
src.__double[ss == 8] = 0L;
max_exp = (ds << 3) - 1; /* signed numbers */
/* have more limited max_exp */
if (buf.exp > max_exp) {
-#ifdef PRT_EXT
- prt_ext("CFI() INT OVERFLOW", &buf);
-#endif PRT_EXT
trap(EIOVFL); /* integer overflow */
buf.exp %= max_exp; /* truncate */
}
/* $Header$ */
/*
- CONVERT FLOAT TO UNSIGNED
+ CONVERT FLOAT TO UNSIGNED (CFU m n)
N.B. The caller must know what it is getting.
A LONG is always returned. If it is an
short newint, max_exp;
extend(&src,&buf,ss); /* get extended format */
-#ifdef PRT_EXT
- prt_ext("CFU() entry:",&buf);
-#endif PRT_EXT
buf.exp--; /* additional bias correction */
if (buf.exp < 1) { /* no conversion needed */
src.__double[ss == 8] = 0L;
}
max_exp = (ds << 3);
if (buf.exp > max_exp) {
-#ifdef PRT_EXT
- prt_ext("CFU() INT OVERFLOW",&buf);
-#endif PRT_EXT
trap(EIOVFL); /* integer overflow */
buf.exp %= max_exp;
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT
+ CONVERT INTEGER TO FLOAT (CIF n 4)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
buf.exp = 17;
result = (_float *) &ss;
}
-#ifdef PRT_STDERR
- fprintf(stderr,"CIF4(ds(%d),ss(%d),src(%D))\n\n",4,ss,i_src);
-#endif
if (i_src == 0) {
*result = (_float) 0L;
return(0L);
if (ss != sizeof(long))
buf.m1 <<= 16;
nrm_ext(&buf); /* adjust mantissa field */
-#ifdef PRT_STDERR
- fprintf(stderr,"CIF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
-#endif
compact(&buf,(_double *) result,4); /* put on stack */
return(*result);
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT
+ CONVERT INTEGER TO FLOAT (CIF n 8)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
i_src = (long) *ipt;
buf.exp = 17;
}
-#ifdef PRT_STDERR
- fprintf(stderr,"CIF8(ds(%d),ss(%d),src(%D))\n\n",8,ss,i_src);
-#endif
if (i_src == 0) {
zrf8(result);
return(*result);
if (ss != sizeof(long))
buf.m1 <<= 16;
nrm_ext(&buf);
-#ifdef PRT_STDERR
- fprintf(stderr,"CIF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
-#endif
compact(&buf,result,8);
return(*result);
}
/* $Header$ */
/*
- COMPARE DOUBLES
+ COMPARE SINGLES (CMF 4)
*/
#include "FP_types.h"
#include "get_put.h"
-short
cmf4(f1,f2)
_float f1,f2;
{
/* $Header$ */
/*
- COMPARE DOUBLES
+ COMPARE DOUBLES (CMF 8)
*/
#include "FP_types.h"
#include "get_put.h"
-short
cmf8(d1,d2)
_double d1,d2;
{
/* $Header$ */
/*
-#define PRT_EXIT
-#define PRT_TRAP
-#define PRT_ENTRY
COMPACT EXTEND FORMAT INTO FLOAT OF PROPER SIZE
*/
SINGLE *SGL;
int exact;
-#ifdef PRT_ENTRY
- prt_ext("enter compact:",f);
-#endif PRT_ENTRY
- if (size == sizeof(_double))
-/********************************************************/
-/*
- COMPACT EXTENDED INTO DOUBLE
-*/
-/********************************************************/
- {
+ if (size == sizeof(_double)) {
+ /*
+ * COMPACT EXTENDED INTO DOUBLE
+ */
if ((f->m1|(f->m2 & DBL_ZERO)) == 0L) {
zrf8(to);
- goto leave;
+ return;
}
f->exp += DBL_BIAS; /* restore proper bias */
if (f->exp > DBL_MAX) {
dbl_over: trap(EFOVFL);
-#ifdef PRT_TRAP
- prt_ext("FCOMPACT DBL OVERFLOW",f);
-#endif PRT_TRAP
f->exp = DBL_MAX;
f->m1 = f->m2 = 0L;
if (error++)
return;
}
else if (f->exp < DBL_MIN) {
-#ifdef PRT_TRAP
- prt_ext("FCOMPACT DBL UNDERFLOW",f);
-#endif PRT_TRAP
trap(EFUNFL);
f->exp = DBL_MIN;
f->m1 = f->m2 = 0L;
if (f->m2 & DBL_ROUNDUP) {
DBL->_s.p2++; /* rounding up */
if (DBL->_s.p2 == 0L) { /* carry out */
-#ifdef PRT_RNDMSG
- write(2,"rounding up lsb\n",16);
-#endif PRT_RNDMSG
DBL->_s.p1.fract++;
if (DBL->_s.p1.fract & DBL_CARRYOUT) { /* carry out */
-#ifdef PRT_RNDMSG
- write(2,"shift due to rounding\n",22);
-#endif PRT_RNDMSG
if (DBL->_s.p1.fract & 01)
DBL->_s.p2 = CARRYBIT;
DBL->_s.p1.fract >>= 1;
if (f->exp > DBL_MAX)
goto dbl_over;
- /* STORE EXPONENT: */
+ /*
+ * STORE EXPONENT:
+ *
+ * 1) clear leading bits (B4-B15)
+ * 2) shift and store exponent
+ */
- /* 1) clear leading bits (B4-B15) */
DBL->_s.p1.fract &= DBL_MASK;
-
- /* 2) shift and store exponent */
f->exp <<= DBL_EXPSHIFT;
DBL->_s.p1.fract |= ((long) f->exp << EXP_STORE);
}
- else
-/********************************************************/
-/*
- COMPACT EXTENDED INTO FLOAT
-*/
-/********************************************************/
- {
+ else {
+ /*
+ * COMPACT EXTENDED INTO FLOAT
+ */
/* local CAST conversion */
SGL = (SINGLE *) to;
if ((f->m1 & SGL_ZERO) == 0L) {
SGL->fract = 0L;
- goto leave;
+ return;
}
f->exp += SGL_BIAS; /* restore bias */
if (f->exp > SGL_MAX) {
sgl_over: trap(EFOVFL);
-#ifdef PRT_TRAP
- prt_ext("FCOMPACT FLOAT OVERFLOW",f);
-#endif PRT_TRAP
f->exp = SGL_MAX;
f->m1 = f->m2 = 0L;
if (error++)
return;
}
else if (f->exp < SGL_MIN) {
-#ifdef PRT_TRAP
- prt_ext("FCOMPACT FLOAT UNDERFLOW",f);
-#endif PRT_TRAP
trap(EFUNFL);
f->exp = SGL_MIN;
f->m1 = f->m2 = 0L;
/* INEXACT(); */
if (f->m1 & SGL_ROUNDUP) {
SGL->fract++;
-#ifdef PRT_RNDMSG
- write(2,"rounding up lsb\n",16);
-#endif PRT_RNDMSG
/* check normal */
if (SGL->fract & SGL_CARRYOUT) {
SGL->fract >>= 1;
if (f->exp > SGL_MAX)
goto sgl_over;
- /* STORE EXPONENT */
- /* 1) clear leading bit of fraction */
- SGL->fract &= SGL_MASK; /* B23-B31 are 0 */
+ /*
+ * STORE EXPONENT:
+ *
+ * 1) clear leading bit of fraction
+ * 2) shift and store exponent
+ */
- /* 2) shift and store exponent */
+ SGL->fract &= SGL_MASK; /* B23-B31 are 0 */
f->exp <<= SGL_EXPSHIFT;
SGL->fract |= ((long) f->exp << EXP_STORE);
}
-/********************************************************/
-/*
- STORE SIGN BIT
-*/
-/********************************************************/
+ /*
+ * STORE SIGN BIT
+ */
if (f->sign) {
SGL = (SINGLE *) to; /* make sure */
SGL->fract |= CARRYBIT;
}
-/********************************************************/
-/*
- STORE MANTISSA
-/*
-/********************************************************/
+
+ /*
+ * STORE MANTISSA
+ */
if (size == sizeof(_double)) {
put4(DBL->_s.p1.fract, (char *) &DBL->_s.p1.fract);
put4(DBL->_s.p2, (char *) &DBL->_s.p2);
}
- else
+ else
put4(SGL->fract, (char *) &SGL->fract);
-
-leave:
-#ifdef PRT_EXIT
- prt_ext("exit compact:",f);
- prt_dbl((DOUBLE *) to,size); getchar();
-#endif PRT_EXIT
- ; /* end of print statement or null statement */
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT
+ CONVERT INTEGER TO FLOAT (CUF n 4)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
buf.exp = 17;
result = (_float *) &ss;
}
-#ifdef PRT_STDERR
- fprintf(stderr,"CUF4(ds(%d),ss(%d),src(%D))\n\n",4,ss,i_src);
-#endif
if (i_src == 0) {
*result = (_float) 0L;
return (_float) 0L;
/* adjust mantissa field */
nrm_ext(&buf);
-#ifdef PRT_STDERR
- fprintf(stderr,"CUF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
-#endif
compact(&buf,(_double *) result,4);
return *result;
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT
+ CONVERT INTEGER TO FLOAT (CUF n 8)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
i_src = (long) *ipt;
buf.exp = 17;
}
-#ifdef PRT_STDERR
- fprintf(stderr,"CUF8(ds(%d),ss(%d),src(%D))\n\n",8,ss,i_src);
-#endif
if (i_src == 0) {
zrf8(&ss);
return;
/* adjust mantissa field */
nrm_ext(&buf);
-#ifdef PRT_STDERR
- fprintf(stderr,"CUF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
-#endif
compact(&buf,(_double *) &ss,8);
}
/* $Header$ */
/*
-#define PRT_EXT
-#define PRT_ALL
DIVIDE EXTENDED FORMAT
*/
This is a routine to do the work.
It is based on the partial products method
and makes no use possible machine instructions
- to divide (hardware dividers). It is intended
- that it be rewritten to do so, but expedieancy
- requires that something be written NOW - and
- this is it.
+ to divide (hardware dividers).
*/
/********************************************************/
unsigned long result[2];
register unsigned long *lp;
-#ifdef PRT_EXT
- fprintf("stderr:start div_ext:\n");
- prt_ext("dividend:",e1);
- prt_ext("divisor :",e2);
-#endif
if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */
e1->exp = 0; /* make sure */
return;
e1->exp -= e2->exp;
e1->exp += 2; /* bias correction */
if (e1->exp < EXT_MIN) {
- error++;
-#ifdef PRT_EXT
- prt_ext("DIV_EXT UNDERFLOW",e1);
-#endif PRT_EXT
+ /*
+ * Exception 8.4 - Underflow
+ */
trap(EFUNFL); /* underflow */
e1->exp = EXT_MIN;
e1->m1 = e1->m2 = 0L;
+ return;
}
if ((e2->m1 | e2->m2) == 0) {
- error++;
-#ifdef PRT_EXT
- prt_ext("DIV_EXT DIV 0.0",e2);
-#endif PRT_EXT
+ /*
+ * Exception 8.2 - Divide by zero
+ */
trap(EFDIVZ);
e1->m1 = e1->m2 = 0L;
e1->exp = EXT_MAX;
- }
- if (error)
return;
+ }
+ if (e1->exp >= EXT_MAX) {
+ /*
+ * Exception 8.3 - Overflow
+ */
+ trap(EFOVFL); /* overflow */
+ e1->exp = EXT_MAX;
+ e1->m1 = e1->m2 = 0L;
+ return;
+ }
/* do division of mantissas */
/* uses partial product method */
/* compare dividend and divisor */
/* if dividend >= divisor add a bit */
/* and subtract divisior from dividend */
-#ifdef PRT_ALL
- prt_ext("dividend:",e1);
- prt_ext("divisor :",e2);
-#endif
if ( (e1->m1 < e2->m1) ||
((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
e1->m1 -= 1; /* carry in */
e1->m1 -= e2->m1; /* do SUBTRACTION */
e1->m2 -= e2->m2; /* SUBTRACTION */
-#ifdef PRT_ALL
- prt_ext("result :",e1);
-#endif
}
-#ifdef PRT_ALL
- fprintf(stderr,"div_ext %d %08X%08X\n\n",64-count,
- result[0],result[1]);
- fflush(stderr);
-#endif
/* shift dividend left one bit OR */
/* IF it equals ZERO we can break out */
break; /* leave loop */
} /* end of divide by subtraction loop */
- /* DISPLAY RESULTS FOR DEBUGGING */
-#ifdef PRT_ALL
- prt_ext("dividend:",e1);
- prt_ext("divisor :",e2);
- fprintf(stderr,"div_ext %d %08X%08X\n",64-count,
- result[0],result[1]);
-#endif
-
if (count > 0) {
lp = result;
if (count > 31) { /* move to higher word */
*lp <<= count;
}
}
- /*
- if (error)
- INEXACT();
- */
+#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->m1 = result[0];
e1->m2 = result[1];
-#ifdef PRT_EXT
- prt_ext("result :",e1);
-#endif
nrm_ext(e1);
-#ifdef PRT_EXT
- prt_ext("after nrm:",e1);
- /*sleep(4);*/
-#endif
}
/* $Header$ */
/*
- DIVIDE TWO FLOATS - SINGLE Precision
+ DIVIDE TWO FLOATS - SINGLE Precision (dvf 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- DIVIDE TWO FLOATS - DOUBLE Precision
+ DIVIDE TWO FLOATS - DOUBLE Precision (DVF 8)
*/
#include "FP_types.h"
/* $Header$ */
/*
-#define PRT_EXIT
-#define PRT_ENTRY
-#define PRT_DBL
CONVERTS FLOATING POINT TO EXTENDED FORMAT
Two sizes of FLOATING Point are known:
unsigned long tmp;
int leadbit = 0;
-#ifdef PRT_ENTRY
- write(2,"entry extend: ",14);
-#ifdef PRT_DBL
- prt_dbl(from,size);
-#else
- write(2,"\n",1);
-#endif PRT_DBL
-#endif PRT_ENTRY
f = (DOUBLE *) from; /* local cast conversion */
if (f->_s.p1.fract == 0L) {
if (size == sizeof(SINGLE)) {
zero: zrf_ext(to);
- goto ready;
+ return;
}
else if (f->_s.p2 == 0L)
goto zero;
to->m1 |= NORMBIT; /* set bit L */
if (leadbit == 0) /* set or clear Leading Bit */
to->m1 &= ~NORMBIT; /* clear bit L */
-ready:
-#ifdef PRT_EXIT
- prt_ext("exit extend:",to)
-#endif PRT_EXIT
- ; /* end of print statement or null statement */
}
/* $Header$ */
/*
- SEPERATE INTO EXPONENT AND FRACTION
+ SEPERATE INTO EXPONENT AND FRACTION (FEF 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- SEPERATE DOUBLE INTO EXPONENT AND FRACTION
+ SEPERATE DOUBLE INTO EXPONENT AND FRACTION (FEF 8)
*/
#include "FP_types.h"
EXTEND buf;
struct fef8_returns *r = (struct fef8_returns *) &s1;
-#ifdef DEBUG
- printf("FEF8(): ");
-#endif DEBUG
extend(&s1,&buf,sizeof(_double));
r->e = buf.exp - 1;
buf.exp = 1;
compact(&buf,&r->f,sizeof(_double));
-#ifdef DEBUG
- printf("exponent = %3d fraction = 0x%08X%08X: ",
- r->f.__double[0],r->f.__double[1]);
- printf("FEF8()\n");
-#endif DEBUG
}
/* $Header$ */
/*
- MULTIPLY AND DISMEMBER PARTS
+ MULTIPLY AND DISMEMBER PARTS (FIF 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- MULTIPLY AND DISMEMBER PARTS
+ MULTIPLY AND DISMEMBER PARTS (FIF 8)
*/
#include "FP_types.h"
/* $Header$ */
/*
- * Multiply Single Precesion Float
+ * Multiply Single Precesion Float (MLF 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- * Multiply Single Precesion Float
+ * Multiply Single Precision Float (MLF 8)
*/
#include "FP_types.h"
mul_ext(e1,e2)
EXTEND *e1,*e2;
{
- register short k,i,j; /* loop control */
+ register int k,i,j; /* loop control */
long unsigned *reg[7];
long unsigned tmp[4];
short unsigned mp[4]; /* multiplier */
short unsigned mc[4]; /* multipcand */
B64 low64,tmp64; /* 64 bit storage */
-#ifdef PRT_EXT
- prt_ext("before MUL_EXT() e1:",e1);
- prt_ext("before MUL_EXT() e2:",e2);
-#endif
/* first save the sign (XOR) */
e1->sign ^= e2->sign;
/* check for overflow */
if (e1->exp >= EXT_MAX) {
-#ifdef PRT_EXT
- prt_ext("EXT_MUL OVERFLOW",e1);
-#endif
trap(EFOVFL);
/* if caught */
/* return signed infinity */
e1->exp = EXT_MAX;
infinity: e1->m1 = e1->m2 =0L;
-#ifdef PRT_EXT
- prt_ext("after MUL_EXT() e1:",e1);
-#endif
return;
}
/* check for underflow */
if (e1->exp < EXT_MIN) {
-#ifdef PRT_EXT
- prt_ext("EXT_MUL UNDERFLOW",e1);
-#endif
trap(EFUNFL);
e1->exp = EXT_MIN;
goto infinity;
mc[1] = (unsigned short) e2->m1;
mc[2] = e2->m2 >> 16;
mc[3] = (unsigned short) e2->m2;
-# ifdef DEBUG
- for(i=0;i<4;i++)
- printf("%04x",mp[i]);
- putchar('\r');
- putchar('\n');
- for(i=0;i<4;i++)
- printf("%04x",mc[i]);
- putchar('\r');
- putchar('\n');
-# endif
/*
* assign pointers
*/
for(j=4;j--;) if (mc[j]) {
k = i+j;
tmp[0] = (long)mp[i] * (long)mc[j];
-# ifdef PRT_EXT2
- printf("%04x * %04x == %08X ",mp[i],mc[j],tmp[0]);
- printf("index == %d ",k);
- printf("register before add == %08X\n",*reg[k]);
- fflush(stdout);
-# endif
-#ifdef PRT_ADD
- printf("REGISTERS-----\n");
- printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
- *reg[0],*reg[2],*reg[4],*reg[6],
- (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
- (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
-# endif
if (b32_add(reg[k],tmp)) {
for(tmp[0] = 0x10000L;k>0;)
if (b32_add(reg[--k],tmp) == 0)
break;
-#ifdef PRT_ADD
- printf("CARRY---------\n");
- printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
- *reg[0],*reg[2],*reg[4],*reg[6],
- (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
- (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
-#endif
}
}
/*
* combine the registers to a total
*/
-#ifdef PRT_ADD
- printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
- *reg[0],*reg[2],*reg[4],*reg[6],
- (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
- (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
-# endif
tmp64.h_32 = (*reg[1]>>16);
tmp64.l_32 = (*reg[1]<<16) + (*reg[3]>>16);
-# ifdef PRT_ALL
- printf("%08X%08X tmp64\n",tmp64.h_32,tmp64.l_32);
- fflush(stdout);
- printf("%08X%08X e1->m1\n",e1->m1,e1->m2);
- fflush(stdout);
-# endif
b64_add((B64 *)&e1->m1,&tmp64);
-# ifdef PRT_ALL
- printf("b64_add:\n");
- printf("%08X%08X e1->m1\n",e1->m1,e1->m2);
- fflush(stdout);
-# endif
tmp64.l_32 = *reg[5]<<16;
tmp64.h_32 = (*reg[5]>>16) + (*reg[3]<<16);
if (b64_add(&low64,&tmp64))
if (++e1->m2 == 0)
e1->m1++;
-# ifdef PRT_ADD
- printf("%08X %08X %08X %08X\n",e1->m1,e1->m2,low64.h_32,low64.l_32);
- fflush(stdout);
-#endif
-#ifdef PRT_EXT
- prt_ext("after MUL_EXT() e1:",e1);
-#endif PRT_EXT
nrm_ext(e1);
-#ifdef PRT_EXT
- prt_ext("after NRM_EXT() e1:",e1);
- sleep(4);
-#endif PRT_EXT
}
/* $Header$ */
/*
- NEGATE A FLOATING POINT
+ NEGATE A FLOATING POINT (NGF 4)
*/
/********************************************************/
/*
*p ^= 0x80;
}
}
-
/* $Header$ */
/*
- NEGATE A FLOATING POINT
+ NEGATE A FLOATING POINT (NGF 8)
*/
/********************************************************/
/*
register unsigned long *mant_2;
/* local CAST conversion */
-#ifdef PRT_EXT
- prt_ext("before NRM_EXT() e1:",e1);
-#endif PRT_EXT
mant_1 = (unsigned long *) &e1->m1;
- /*
- THIS RESULTS IN A BAD CODE !!!!
- ANOTHER BUG IN EM CODE MAYBE????
mant_2 = (unsigned long *) &e1->m2;
- */
- /* statement that works */
- mant_2 = mant_1 + 1;
/* we assume that the mantissa != 0 */
/* if it is then just return */
*mant_1-- = 0L;
e1->exp -= 32;
}
-#ifdef OLD
- /* check that e1->m1 is not too large */
- if (*mant_1 & CARRYBIT) { /* carry occured */
- e1->exp++; /* increase exponent */
- *mant_2 >>= 1; /* right shift mantissa */
- if ((short) *mant_1 & 01)
- *mant_2 |= CARRYBIT;
- *mant_1 >>= 1;
- }
-#endif
-
-
while ((*mant_1 & NORMBIT) == 0) {
e1->exp--;
*mant_1 <<= 1;
}
*mant_2 <<= 1;
}
-#ifdef PRT_EXT
- prt_ext("after NRM_EXT() e1:",e1);
-#endif PRT_EXT
}
/* $Header$ */
/*
- SUBTRACT TWO FLOATS - SINGLE Precision
+ SUBTRACT TWO FLOATS - SINGLE Precision (SBF 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- SUBTRACT TWO FLOATS - DOUBLE Precision
+ SUBTRACT TWO FLOATS - DOUBLE Precision (SBF 8)
*/
#include "FP_types.h"
/* s2 = -s2; */
char unsigned *p; /* sufficient to access sign bit */
-#ifdef PRT_EXT
- fprintf(stderr,"SBF8 ():\n");
-#endif
if (s2.__double[0] == 0 && s2.__double[1] == 0) {
return s1;
}
s1 = adf8(s2,s1); /* add and return result */
return(s1);
}
-
/* $Header$ */
/*
-#define PRT_EXT
SHIFT TWO EXTENDED NUMBERS INTO PROPER
ALIGNMENT FOR ADDITION (exponents are equal)
*/
EXTEND *e1,*e2;
{
register EXTEND *s;
- register short diff;
+ register int diff;
long tmp;
-#ifdef PRT_EXT
- prt_ext("enter sft_ext e1:",e1);
- prt_ext("enter sft_ext e2:",e2);
-#endif PRT_EXT
diff = e1->exp - e2->exp;
if (!diff)
s->m2 >>= diff;
s->m2 |= tmp;
}
-#ifdef PRT_EXT
- prt_ext("exit sft_ext e1:",e1);
- prt_ext("exit sft_ext e2:",e2);
-#endif PRT_EXT
}
b64_sft(e1,n)
B64 *e1;
-short n;
+int n;
{
if (n>0) do { /* RIGHT shift n bits */
e1->l_32 >>= 1; /* shift 64 bits */
/* $Header$ */
/*
-#define PRT_EXT
- SUBTRACT EXTENDED FORMAT
+ SUBTRACT 2 EXTENDED FORMAT NUMBERS
*/
- /* assumes that e1 >= e2 on entry */
- /* no test is made to check this */
- /* so make sure yourself */
+ /*
+ * 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;
{
-#ifdef PRT_EXT
- prt_ext("before SUB_EXT() e1:",e1);
- prt_ext("before SUB_EXT() e2:",e2);
-#endif PRT_EXT
if (e2->m2 > e1->m2)
e1->m1 -= 1; /* carry in */
e1->m1 -= e2->m1;
e1->m2 -= e2->m2;
-#ifdef PRT_EXT
- prt_ext("after SUB_EXT() e1:",e1);
-#endif PRT_EXT
nrm_ext(e1);
-#ifdef PRT_EXT
- prt_ext("after NRM_EXT() e1:",e1);
-#endif PRT_EXT
}
/* $Header$ */
/*
- return a zero float
+ return a zero float (ZRF 4)
*/
#include "FP_types.h"
/* $Header$ */
/*
- return a zero double
+ return a zero double (ZRF 8)
*/
#include "FP_types.h"
EXTEND *e;
{
register short *ipt;
- register short i;
- register short zero = 0;
+ register int i;
/* local CAST conversion */
ipt = (short *) e;
i = sizeof(EXTEND)/sizeof(short);
while (i--)
- *ipt++ = zero;
+ *ipt++ = 0;
}
public / ack.git / commitdiff