/********************************************************/
/*
THESE STRUCTURES ARE USED TO ADDRESS THE INDIVIDUAL
- PARTS OF THE FLOATING NUMBER REPRESENTATIONS.
+ PARTS OF THE FLOATING POINT NUMBER REPRESENTATIONS.
THREE STRUCTURES ARE DEFINED:
SINGLE: single precision floating format
*/
/********************************************************/
-#ifdef EXT_DEBUG
-#ifndef __FPSTDIO
-#define __FPSTDIO
-#include <stdio.h>
-#endif
-#endif
-
#ifndef __FPTYPES
#define __FPTYPES
-typedef unsigned long _float;
-typedef union {
- double _dbl;
- unsigned long __double[2];
-} _double;
+typedef struct {
+ unsigned long h_32; /* higher 32 bits of 64 */
+ unsigned long l_32; /* lower 32 bits of 64 */
+} B64;
-typedef union { /* address parts of float */
- /* field to extract exponent */
- short sgl_exp[sizeof(float)/sizeof(short)];
- /* same as fract below */
- _float f[sizeof(float)/sizeof(_float)];
- unsigned long fract;
-} SINGLE;
+typedef unsigned long SINGLE;
-typedef union {
- /* field to extract exponent */
- short dbl_exp[sizeof(double)/sizeof(short)];
- /* for normal syntax use */
- _double _f8[sizeof(double)/sizeof(_double)];
- /* a float in a double - for returns */
- _float _f4[sizeof(double)/sizeof(_float)];
- /* to be deleted eventually */
- struct { /* address parts of float */
- SINGLE p1; /* part one */
- unsigned long p2; /* part two */
- } _s;
-} DOUBLE;
+typedef struct {
+ unsigned long d[2];
+} DOUBLE;
typedef struct { /* expanded float format */
short sign;
short exp;
- unsigned long m1;
- unsigned long m2; /* includes guard byte */
+ B64 mantissa;
+#define m1 mantissa.h_32
+#define m2 mantissa.l_32
} EXTEND;
+
+struct fef4_returns {
+ int e;
+ SINGLE f;
+};
+
+struct fef8_returns {
+ int e;
+ DOUBLE f;
+};
+
+struct fif4_returns {
+ SINGLE ipart;
+ SINGLE fpart;
+};
+
+struct fif8_returns {
+ DOUBLE ipart;
+ DOUBLE fpart;
+};
+
+#if __STDC__
+#define _PROTOTYPE(function, params) function params
+#else
+#define _PROTOTYPE(function, params) function()
+#endif
+_PROTOTYPE( void add_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void mul_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void div_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sub_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sft_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void nrm_ext, (EXTEND *e1));
+_PROTOTYPE( void zrf_ext, (EXTEND *e1));
+_PROTOTYPE( void extend, (unsigned long *from, EXTEND *to, int size));
+_PROTOTYPE( void compact, (EXTEND *from, unsigned long *to, int size));
+_PROTOTYPE( void _fptrp, (int));
+_PROTOTYPE( SINGLE adf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( DOUBLE adf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( SINGLE sbf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( DOUBLE sbf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( SINGLE dvf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( DOUBLE dvf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( SINGLE mlf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( DOUBLE mlf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( SINGLE ngf4, (SINGLE f));
+_PROTOTYPE( DOUBLE ngf8, (DOUBLE f));
+_PROTOTYPE( void zrf4, (SINGLE *l));
+_PROTOTYPE( void zrf8, (DOUBLE *z));
+_PROTOTYPE( SINGLE cff4, (DOUBLE src));
+_PROTOTYPE( DOUBLE cff8, (SINGLE src));
+_PROTOTYPE( SINGLE cif4, (int ss, long src));
+_PROTOTYPE( DOUBLE cif8, (int ss, long src));
+_PROTOTYPE( SINGLE cuf4, (int ss, long src));
+_PROTOTYPE( DOUBLE cuf8, (int ss, long src));
+_PROTOTYPE( long cfu, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( long cfi, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( int cmf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( int cmf8, (DOUBLE d1, DOUBLE d2));
+_PROTOTYPE( void fef4, (struct fef4_returns *r, SINGLE s1));
+_PROTOTYPE( void fef8, (struct fef8_returns *r, DOUBLE s1));
+_PROTOTYPE( void fif4, (struct fif4_returns *p, SINGLE x, SINGLE y));
+_PROTOTYPE( void fif8, (struct fif8_returns *p, DOUBLE x, DOUBLE y));
+
+_PROTOTYPE( void b64_sft, (B64 *, int));
+_PROTOTYPE( void b64_lsft, (B64 *));
+_PROTOTYPE( void b64_rsft, (B64 *));
+_PROTOTYPE( int b64_add, (B64 *, B64 *));
#endif
#include "FP_types.h"
+void
add_ext(e1,e2)
register EXTEND *e1,*e2;
{
}
}
else {
- if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
- b64_rsft(&e1->m1); /* shift mantissa one bit RIGHT */
+ if (b64_add(&e1->mantissa,&e2->mantissa)) { /* addition carry */
+ b64_rsft(&e1->mantissa); /* shift mantissa one bit RIGHT */
e1->m1 |= 0x80000000L; /* set max bit */
e1->exp++; /* increase the exponent */
}
# include <stdio.h>
# endif
-# include "adder.h"
+# include "FP_types.h"
# define UNKNOWN -1
# define TRUE 1
# define FALSE 0
/*
* add 64 bits
*/
+int
b64_add(e1,e2)
/*
* pointers to 64 bit 'registers'
# endif
if ((carry) && (++e1->h_32 == 0))
return(TRUE); /* had a 64 bit overflow */
- else
- return(overflow); /* return status from higher add */
+ return(overflow); /* return status from higher add */
}
#include "FP_types.h"
-_float
+SINGLE
adf4(s2,s1)
-_float s1,s2;
+SINGLE s1,s2;
{
EXTEND e1,e2;
int swap = 0;
- if (s1 == (_float) 0) {
+ if (s1 == (SINGLE) 0) {
s1 = s2;
return s1;
}
- if (s2 == (_float) 0) {
+ if (s2 == (SINGLE) 0) {
return s1;
}
- extend((_double *)&s1,&e1,sizeof(SINGLE));
- extend((_double *)&s2,&e2,sizeof(SINGLE));
+ extend(&s1,&e1,sizeof(SINGLE));
+ extend(&s2,&e2,sizeof(SINGLE));
add_ext(&e1,&e2);
- compact(&e1,(_double *)&s1,sizeof(SINGLE));
+ compact(&e1,&s1,sizeof(SINGLE));
return s1;
}
#include "FP_types.h"
-_double
+DOUBLE
adf8(s2,s1)
-_double s1,s2;
+DOUBLE s1,s2;
{
EXTEND e1,e2;
- if (s1.__double[0] == 0 && s1.__double[1] == 0) {
+ if (s1.d[0] == 0 && s1.d[1] == 0) {
s1 = s2;
return s1;
}
- if (s2.__double[0] == 0 && s2.__double[1] == 0) {
+ if (s2.d[0] == 0 && s2.d[1] == 0) {
return s1;
}
- extend(&s1,&e1,sizeof(_double));
- extend(&s2,&e2,sizeof(_double));
+ extend(&s1.d[0],&e1,sizeof(DOUBLE));
+ extend(&s2.d[0],&e2,sizeof(DOUBLE));
add_ext(&e1,&e2);
- compact(&e1,&s1,sizeof(_double));
+ compact(&e1,&s1.d[0],sizeof(DOUBLE));
return s1;
}
/* $Header$ */
/*
- CONVERT DOUBLE TO FLOAT (CFF 8 4)
+ CONVERT DOUBLE TO SINGLE (CFF 8 4)
This routine works quite simply. A floating point
of size 08 is converted to extended format.
#include "FP_types.h"
-_float
+SINGLE
cff4(src)
-_double src; /* the source itself - THIS TIME it's DOUBLE */
+DOUBLE src; /* the source itself - THIS TIME it's DOUBLE */
{
EXTEND buf;
- extend(&src,&buf,8); /* no matter what */
- compact(&buf,(_double *) &(src.__double[1]),4);
- return *(_float *)&(src.__double[1]);
+ extend(&src.d[0],&buf,sizeof(DOUBLE)); /* no matter what */
+ compact(&buf,&(src.d[1]),sizeof(SINGLE));
+ return *(SINGLE *)&(src.d[1]);
}
/* $Header$ */
/*
- CONVERT FLOAT TO DOUBLE (CFF 4 8)
+ CONVERT SINGLE TO DOUBLE (CFF 4 8)
This routine works quite simply. A floating point
of size 04 is converted to extended format.
#include "FP_types.h"
-_double
+DOUBLE
cff8(src)
-_float src;
+SINGLE src;
{
EXTEND buf;
- extend((_double *) &src,&buf,4); /* no matter what */
- compact(&buf,(_double *) &src,8);
- return *(_double *) &src;
+ extend(&src,&buf,sizeof(SINGLE)); /* no matter what */
+ compact(&buf, &src,sizeof(DOUBLE));
+ return *(DOUBLE *) ((void *) &src);
}
cfi(ds,ss,src)
int ds; /* destination size (2 or 4) */
int ss; /* source size (4 or 8) */
-_double src; /* assume worst case */
+DOUBLE src; /* assume worst case */
{
EXTEND buf;
long new;
short max_exp;
- extend(&src,&buf,ss); /* get extended format */
+ extend(&src.d[0],&buf,ss); /* get extended format */
if (buf.exp < 0) { /* no conversion needed */
- src.__double[ss == 8] = 0L;
+ src.d[ss == 8] = 0L;
return(0L);
}
max_exp = (ds << 3) - 2; /* signed numbers */
if (buf.sign)
new = -new;
done:
- src.__double[ss == 8] = new;
+ src.d[ss == 8] = new;
return(new);
}
cfu(ds,ss,src)
int ds; /* destination size (2 or 4) */
int ss; /* source size (4 or 8) */
-_double src; /* assume worst case */
+DOUBLE src; /* assume worst case */
{
EXTEND buf;
long new;
short newint, max_exp;
- extend(&src,&buf,ss); /* get extended format */
+ extend(&src.d[0],&buf,ss); /* get extended format */
if (buf.exp < 0) { /* no conversion needed */
- src.__double[ss == 8] = 0L;
+ src.d[ss == 8] = 0L;
return(0L);
}
max_exp = (ds << 3) - 1;
}
new = buf.m1 >> (31-buf.exp);
done:
- src.__double[ss == 8] = new;
+ src.d[ss == 8] = new;
return(new);
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT (CIF n 4)
+ CONVERT INTEGER TO SINGLE (CIF n 4)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
#include "FP_types.h"
-_float
+SINGLE
cif4(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
EXTEND buf;
short *ipt;
long i_src;
- _float *result;
+ SINGLE *result;
zrf_ext(&buf);
if (ss == sizeof(long)) {
buf.exp = 31;
i_src = src;
- result = (_float *) &src;
+ result = (SINGLE *) &src;
}
else {
ipt = (short *) &src;
i_src = (long) *ipt;
buf.exp = 15;
- result = (_float *) &ss;
+ result = (SINGLE *) &ss;
}
if (i_src == 0) {
- *result = (_float) 0L;
+ *result = (SINGLE) 0L;
return(0L);
}
/* ESTABLISHED THAT src != 0 */
if (ss != sizeof(long))
buf.m1 <<= 16;
nrm_ext(&buf); /* adjust mantissa field */
- compact(&buf,(_double *) result,4); /* put on stack */
+ compact(&buf, result,sizeof(SINGLE)); /* put on stack */
return(*result);
}
#include "FP_types.h"
-_double
+DOUBLE
cif8(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
{
EXTEND buf;
- _double *result; /* for return value */
+ DOUBLE *result; /* for return value */
short *ipt;
long i_src;
- result = (_double *) &ss; /* always */
+ result = (DOUBLE *) ((void *) &ss); /* always */
zrf_ext(&buf);
if (ss == sizeof(long)) {
buf.exp = 31;
if (ss != sizeof(long))
buf.m1 <<= 16;
nrm_ext(&buf);
- compact(&buf,result,8);
+ compact(&buf,&result->d[0],8);
return(*result);
}
int
cmf4(f1,f2)
-_float f1,f2;
+SINGLE f1,f2;
{
/*
* return ((f1 < f2) ? 1 : (f1 - f2))
#include "FP_types.h"
#include "get_put.h"
+int
cmf8(d1,d2)
-_double d1,d2;
+DOUBLE d1,d2;
{
#define SIGN(x) (((x) < 0) ? -1 : 1)
/*
# include "FP_types.h"
# include "get_put.h"
+void
compact(f,to,size)
EXTEND *f;
-_double *to;
+unsigned long *to;
int size;
{
int error = 0;
- if (size == sizeof(_double)) {
+ if (size == sizeof(DOUBLE)) {
/*
* COMPACT EXTENDED INTO DOUBLE
*/
- DOUBLE *DBL;
+ DOUBLE *DBL = (DOUBLE *) (void *) to;
if ((f->m1|(f->m2 & DBL_ZERO)) == 0L) {
- zrf8(to);
+ zrf8(DBL);
return;
}
f->exp += DBL_BIAS; /* restore proper bias */
return;
}
else if (f->exp < DBL_MIN) {
- b64_rsft(&(f->m1));
+ b64_rsft(&(f->mantissa));
if (f->exp < 0) {
- b64_sft(&(f->m1), -f->exp);
+ b64_sft(&(f->mantissa), -f->exp);
f->exp = 0;
}
/* underflow ??? */
}
/* local CAST conversion */
- DBL = (DOUBLE *) to;
/* because of special format shift only 10 bits */
/* bit shift mantissa 10 bits */
/* first align within words, then do store operation */
- DBL->_s.p1.fract = f->m1 >> DBL_RUNPACK; /* plus 22 == 32 */
- DBL->_s.p2 = f->m2 >> DBL_RUNPACK; /* plus 22 == 32 */
- DBL->_s.p2 |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
+ DBL->d[0] = f->m1 >> DBL_RUNPACK; /* plus 22 == 32 */
+ DBL->d[1] = f->m2 >> DBL_RUNPACK; /* plus 22 == 32 */
+ DBL->d[1] |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
/* if not exact then round to nearest */
/* on a tie, round to even */
if (((f->m2 & DBL_EXACT) > DBL_ROUNDUP)
|| ((f->m2 & DBL_EXACT) == DBL_ROUNDUP
&& (f->m2 & (DBL_ROUNDUP << 1)))) {
- DBL->_s.p2++; /* rounding up */
- if (DBL->_s.p2 == 0L) { /* carry out */
- DBL->_s.p1.fract++;
+ DBL->d[1]++; /* rounding up */
+ if (DBL->d[1] == 0L) { /* carry out */
+ DBL->d[0]++;
- if (f->exp == 0 && (DBL->_s.p1.fract & ~DBL_MASK)) {
+ if (f->exp == 0 && (DBL->d[0] & ~DBL_MASK)) {
f->exp++;
}
- if (DBL->_s.p1.fract & DBL_CARRYOUT) { /* carry out */
- if (DBL->_s.p1.fract & 01)
- DBL->_s.p2 = CARRYBIT;
- DBL->_s.p1.fract >>= 1;
+ if (DBL->d[0] & DBL_CARRYOUT) { /* carry out */
+ if (DBL->d[0] & 01)
+ DBL->d[1] = CARRYBIT;
+ DBL->d[0] >>= 1;
f->exp++;
}
}
* 2) shift and store exponent
*/
- DBL->_s.p1.fract &= DBL_MASK;
- DBL->_s.p1.fract |=
+ DBL->d[0] &= DBL_MASK;
+ DBL->d[0] |=
((long) (f->exp << DBL_EXPSHIFT) << EXP_STORE);
if (f->sign)
- DBL->_s.p1.fract |= CARRYBIT;
+ DBL->d[0] |= CARRYBIT;
/*
* STORE MANTISSA
*/
#if FL_MSL_AT_LOW_ADDRESS
- put4(DBL->_s.p1.fract, (char *) &DBL->_s.p1.fract);
- put4(DBL->_s.p2, (char *) &DBL->_s.p2);
+ put4(DBL->d[0], (char *) &DBL->d[0]);
+ put4(DBL->d[1], (char *) &DBL->d[1]);
#else
{ unsigned long l;
- put4(DBL->_s.p2, (char *) &l);
- put4(DBL->_s.p1.fract, (char *) &DBL->_s.p2);
- DBL->_s.p1.fract = l;
+ put4(DBL->d[1], (char *) &l);
+ put4(DBL->d[0], (char *) &DBL->d[1]);
+ DBL->d[0] = l;
}
#endif
}
SINGLE *SGL;
/* local CAST conversion */
- SGL = (SINGLE *) to;
+ SGL = (SINGLE *) (void *) to;
if ((f->m1 & SGL_ZERO) == 0L) {
- SGL->fract = 0L;
+ *SGL = 0L;
return;
}
f->exp += SGL_BIAS; /* restore bias */
return;
}
else if (f->exp < SGL_MIN) {
- b64_rsft(&(f->m1));
+ b64_rsft(&(f->mantissa));
if (f->exp < 0) {
- b64_sft(&(f->m1), -f->exp);
+ b64_sft(&(f->mantissa), -f->exp);
f->exp = 0;
}
/* underflow ??? */
}
/* shift mantissa and store */
- SGL->fract = (f->m1 >> SGL_RUNPACK);
+ *SGL = (f->m1 >> SGL_RUNPACK);
/* check for rounding to nearest */
/* on a tie, round to even */
if (((f->m1 & SGL_EXACT) > SGL_ROUNDUP)
|| ((f->m1 & SGL_EXACT) == SGL_ROUNDUP
&& (f->m1 & (SGL_ROUNDUP << 1)))) {
- SGL->fract++;
- if (f->exp == 0 && (SGL->fract & ~SGL_MASK)) {
+ (*SGL)++;
+ if (f->exp == 0 && (*SGL & ~SGL_MASK)) {
f->exp++;
}
/* check normal */
- if (SGL->fract & SGL_CARRYOUT) {
- SGL->fract >>= 1;
+ if (*SGL & SGL_CARRYOUT) {
+ *SGL >>= 1;
f->exp++;
}
if (f->exp > SGL_MAX)
* 2) shift and store exponent
*/
- SGL->fract &= SGL_MASK; /* B23-B31 are 0 */
- SGL->fract |=
- ((long) (f->exp << SGL_EXPSHIFT) << EXP_STORE);
+ *SGL &= SGL_MASK; /* B23-B31 are 0 */
+ *SGL |= ((long) (f->exp << SGL_EXPSHIFT) << EXP_STORE);
if (f->sign)
- SGL->fract |= CARRYBIT;
+ *SGL |= CARRYBIT;
/*
* STORE MANTISSA
*/
- put4(SGL->fract, (char *) &SGL->fract);
+ put4(*SGL, (char *) &SGL);
}
}
/* $Header$ */
/*
- CONVERT INTEGER TO FLOAT (CUF n 4)
+ CONVERT INTEGER TO SINGLE (CUF n 4)
THIS ROUTINE WORKS BY FILLING AN EXTENDED
WITH THE INTEGER VALUE IN EXTENDED FORMAT
#include "FP_types.h"
-_float
+SINGLE
cuf4(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
{
EXTEND buf;
short *ipt;
- _float *result;
+ SINGLE *result;
long i_src;
zrf_ext(&buf);
if (ss == sizeof(long)) {
buf.exp = 31;
i_src = src;
- result = (_float *) &src;
+ result = (SINGLE *) &src;
}
else {
ipt = (short *) &src;
i_src = (long) *ipt;
buf.exp = 15;
- result = (_float *) &ss;
+ result = (SINGLE *) ((void *) &ss);
}
if (i_src == 0) {
- *result = (_float) 0L;
- return (_float) 0L;
+ *result = (SINGLE) 0L;
+ return (SINGLE) 0L;
}
/* ESTABLISHED THAT src != 0 */
/* adjust mantissa field */
nrm_ext(&buf);
- compact(&buf,(_double *) result,4);
+ compact(&buf,result,4);
return *result;
}
#include "FP_types.h"
-_double
+DOUBLE
cuf8(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
buf.exp = 15;
}
if (i_src == 0) {
- zrf8(&ss);
+ zrf8((DOUBLE *)((void *)&ss));
return;
}
/* ESTABLISHED THAT src != 0 */
/* adjust mantissa field */
nrm_ext(&buf);
- compact(&buf,(_double *) &ss,8);
- return *((_double *) &ss);
+ compact(&buf,(unsigned long *) (void *)&ss,8);
+ return *((DOUBLE *) (void *)&ss);
}
*/
/********************************************************/
+void
div_ext(e1,e2)
EXTEND *e1,*e2;
{
- short error = 0;
- unsigned long result[2];
+ short error = 0;
+ B64 result;
register unsigned long *lp;
#ifndef USE_DIVIDE
- short count;
+ short count;
#else
- unsigned short u[9], v[5];
- register int j;
- register unsigned short *u_p = u;
- int maxv = 4;
+ 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) {
* that m1 is quaranteed to be larger if its
* maximum bit is set
*/
- b64_rsft(&e1->m1); /* 64 bit shift right */
- b64_rsft(&e2->m1); /* 64 bit shift right */
+ b64_rsft(&e1->mantissa); /* 64 bit shift right */
+ b64_rsft(&e2->mantissa); /* 64 bit shift right */
e1->exp++;
e2->exp++;
#endif
/* init control variables */
count = 64;
- result[0] = 0L;
- result[1] = 0L;
+ result.h_32 = 0L;
+ result.l_32 = 0L;
/* partial product division loop */
/* first left shift result 1 bit */
/* this is ALWAYS done */
- b64_lsft(result);
+ b64_lsft(&result);
/* compare dividend and divisor */
/* if dividend >= divisor add a bit */
; /* null statement */
/* i.e., don't add or subtract */
else {
- result[1]++; /* ADD */
+ result.l_32++; /* ADD */
if (e2->m2 > e1->m2)
e1->m1 -= 1; /* carry in */
e1->m1 -= e2->m1; /* do SUBTRACTION */
error = ((*lp | *(lp+1)) != 0L) ? 1 : 0;
if (error) { /* more work */
/* assume max bit == 0 (see above) */
- b64_lsft(&e1->m1);
+ b64_lsft(&e1->mantissa);
continue;
}
else
} /* end of divide by subtraction loop */
if (count > 0) {
- lp = result;
+ lp = &result.h_32;
if (count > 31) { /* move to higher word */
*lp = *(lp+1);
count -= 32;
}
if (*lp)
*lp <<= count; /* shift rest of way */
- lp++; /* == &result[1] */
+ lp++; /* == &result.l_32 */
if (*lp) {
- result[0] |= (*lp >> 32-count);
+ result.h_32 |= (*lp >> 32-count);
*lp <<= count;
}
}
#else /* USE_DIVIDE */
u[4] = (e1->m2 & 1) << 15;
- b64_rsft(&(e1->m1));
+ b64_rsft(&(e1->mantissa));
u[0] = e1->m1 >> 16;
u[1] = e1->m1;
u[2] = e1->m2 >> 16;
v[3] = e2->m2 >> 16;
v[4] = e2->m2;
while (! v[maxv]) maxv--;
- result[0] = 0;
- result[1] = 0;
- lp = result;
+ result.h_32 = 0;
+ result.l_32 = 0;
+ lp = &result.h_32;
/*
* Use an algorithm of Knuth (The art of programming, Seminumerical
*/
INEXACT();
#endif
- e1->m1 = result[0];
- e1->m2 = result[1];
+ e1->mantissa = result;
nrm_ext(e1);
if (e1->exp < EXT_MIN) {
/* $Header$ */
/*
- DIVIDE TWO FLOATS - SINGLE Precision (dvf 4)
+ DIVIDE TWO SINGLES - SINGLE Precision (dvf 4)
*/
#include "FP_types.h"
-_float
+SINGLE
dvf4(s2,s1)
-_float s1,s2;
+SINGLE s1,s2;
{
EXTEND e1,e2;
- extend((_double *)&s1,&e1,sizeof(_float));
- extend((_double *)&s2,&e2,sizeof(_float));
+ extend(&s1,&e1,sizeof(SINGLE));
+ extend(&s2,&e2,sizeof(SINGLE));
/* do a divide */
div_ext(&e1,&e2);
- compact(&e1,(_double *)&s1,sizeof(_float));
+ compact(&e1,&s1,sizeof(SINGLE));
return s1;
}
#include "FP_types.h"
-_double
+DOUBLE
dvf8(s2,s1)
-_double s1,s2;
+DOUBLE s1,s2;
{
EXTEND e1,e2;
- extend(&s1,&e1,sizeof(_double));
- extend(&s2,&e2,sizeof(_double));
+ extend(&s1.d[0],&e1,sizeof(DOUBLE));
+ extend(&s2.d[0],&e2,sizeof(DOUBLE));
/* do a divide */
div_ext(&e1,&e2);
- compact(&e1,&s1,sizeof(_double));
+ compact(&e1,&s1.d[0],sizeof(DOUBLE));
return s1;
}
#include "get_put.h"
/********************************************************/
+void
extend(from,to,size)
-_double *from;
+unsigned long *from;
EXTEND *to;
int size;
{
#include "FP_types.h"
-struct fef4_returns {
- int e;
- _float f;
-};
-
+void
fef4(r,s1)
-_float s1;
+SINGLE s1;
struct fef4_returns *r;
{
EXTEND buf;
to itself (see table)
*/
- extend((_double *) &s1,&buf,sizeof(_float));
+ extend(&s1,&buf,sizeof(SINGLE));
if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
p->e = 0;
}
p->e = buf.exp+1;
buf.exp = -1;
}
- compact(&buf,(_double *) &p->f,sizeof(_float));
+ compact(&buf,&p->f,sizeof(SINGLE));
}
#include "FP_types.h"
-struct fef8_returns {
- int e;
- _double f;
-};
-
+void
fef8(r, s1)
-_double s1;
+DOUBLE s1;
struct fef8_returns *r;
{
EXTEND buf;
to itself (see table)
*/
- extend(&s1,&buf,sizeof(_double));
+ extend(&s1.d[0],&buf,sizeof(DOUBLE));
if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
p->e = 0;
}
p->e = buf.exp + 1;
buf.exp = -1;
}
- compact(&buf,&p->f,sizeof(_double));
+ compact(&buf,&p->f.d[0],sizeof(DOUBLE));
}
#include "FP_types.h"
#include "FP_shift.h"
-_float sbf4();
-
-struct fif4_returns {
- _float ipart;
- _float fpart;
-};
-
+void
fif4(p,x,y)
-_float x,y;
+SINGLE x,y;
struct fif4_returns *p;
{
EXTEND e1,e2;
- extend((_double *)&y,&e1,sizeof(_float));
- extend((_double *)&x,&e2,sizeof(_float));
+ extend(&y,&e1,sizeof(SINGLE));
+ extend(&x,&e2,sizeof(SINGLE));
/* do a multiply */
mul_ext(&e1,&e2);
e2 = e1;
- compact(&e2, (_double *)&y, sizeof(_float));
+ compact(&e2,&y,sizeof(SINGLE));
if (e1.exp < 0) {
p->ipart = 0;
p->fpart = y;
p->fpart = 0;
return;
}
- b64_sft(&e1.m1, 63 - e1.exp);
- b64_sft(&e1.m1, e1.exp - 63); /* "loose" low order bits */
- compact(&e1,(_double *) &(p->ipart), sizeof(SINGLE));
+ b64_sft(&e1.mantissa, 63 - e1.exp);
+ b64_sft(&e1.mantissa, e1.exp - 63); /* "loose" low order bits */
+ compact(&e1,&(p->ipart),sizeof(SINGLE));
p->fpart = sbf4(p->ipart, y);
}
#include "FP_types.h"
#include "FP_shift.h"
-_double sbf8();
-
-struct fif8_returns {
- _double ipart;
- _double fpart;
-};
-
+void
fif8(p,x,y)
-_double x,y;
+DOUBLE x,y;
struct fif8_returns *p;
{
EXTEND e1,e2;
- extend(&y,&e1,sizeof(_double));
- extend(&x,&e2,sizeof(_double));
+ extend(&y.d[0],&e1,sizeof(DOUBLE));
+ extend(&x.d[0],&e2,sizeof(DOUBLE));
/* do a multiply */
mul_ext(&e1,&e2);
e2 = e1;
- compact(&e2, &y, sizeof(_double));
+ compact(&e2, &y.d[0], sizeof(DOUBLE));
if (e1.exp < 0) {
- p->ipart.__double[0] = 0;
- p->ipart.__double[1] = 0;
+ p->ipart.d[0] = 0;
+ p->ipart.d[1] = 0;
p->fpart = y;
return;
}
if (e1.exp > 62 - DBL_M1LEFT) {
p->ipart = y;
- p->fpart.__double[0] = 0;
- p->fpart.__double[1] = 0;
+ p->fpart.d[0] = 0;
+ p->fpart.d[1] = 0;
return;
}
- b64_sft(&e1.m1, 63 - e1.exp);
- b64_sft(&e1.m1, e1.exp - 63); /* "loose" low order bits */
- compact(&e1, &(p->ipart), sizeof(DOUBLE));
+ b64_sft(&e1.mantissa, 63 - e1.exp);
+ b64_sft(&e1.mantissa, e1.exp - 63); /* "loose" low order bits */
+ compact(&e1, &(p->ipart.d[0]), sizeof(DOUBLE));
p->fpart = sbf8(p->ipart, y);
}
#include "FP_types.h"
-_float
+SINGLE
mlf4(s2,s1)
-_float s1,s2;
+SINGLE s1,s2;
{
EXTEND e1,e2;
- extend((_double *)&s1,&e1,sizeof(_float));
- extend((_double *)&s2,&e2,sizeof(_float));
+ extend(&s1,&e1,sizeof(SINGLE));
+ extend(&s2,&e2,sizeof(SINGLE));
/* do a multiply */
mul_ext(&e1,&e2);
- compact(&e1,(_double *)&s1,sizeof(_float));
+ compact(&e1,&s1,sizeof(SINGLE));
return(s1);
}
/* $Header$ */
/*
- * Multiply Single Precision Float (MLF 8)
+ * Multiply Double Precision Float (MLF 8)
*/
#include "FP_types.h"
-_double
+DOUBLE
mlf8(s2,s1)
-_double s1,s2;
+DOUBLE s1,s2;
{
EXTEND e1,e2;
- extend(&s1,&e1,sizeof(_double));
- extend(&s2,&e2,sizeof(_double));
+ extend(&s1.d[0],&e1,sizeof(DOUBLE));
+ extend(&s2.d[0],&e2,sizeof(DOUBLE));
/* do a multiply */
mul_ext(&e1,&e2);
- compact(&e1,&s1,sizeof(_double));
+ compact(&e1,&s1.d[0],sizeof(DOUBLE));
return(s1);
}
ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
*/
-# include "adder.h"
# include "FP_bias.h"
# include "FP_trap.h"
# include "FP_types.h"
# include "FP_shift.h"
+void
mul_ext(e1,e2)
EXTEND *e1,*e2;
{
#include "get_put.h"
#define OFF ((FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
-_float
+SINGLE
ngf4(f)
-_float f;
+SINGLE f;
{
unsigned char *p;
- if (f != (_float) 0) {
+ if (f != (SINGLE) 0) {
p = (unsigned char *) &f + OFF;
*p ^= 0x80;
}
#define OFF ((FL_MSL_AT_LOW_ADDRESS ? 0 : 4) + (FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
-_double
+DOUBLE
ngf8(f)
-_double f;
+DOUBLE f;
{
unsigned char *p;
- if (f.__double[0] != 0 || f.__double[1] != 0) {
+ if (f.d[0] != 0 || f.d[1] != 0) {
p = (unsigned char *) &f + OFF;
*p ^= 0x80;
}
#include "FP_shift.h"
#include "FP_types.h"
+void
nrm_ext(e1)
EXTEND *e1;
{
cnt--;
}
e1->exp += cnt;
- b64_sft(&(e1->m1), cnt);
+ b64_sft(&(e1->mantissa), cnt);
}
}
/*
- (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
See the copyright notice in the ACK home directory, in the file "Copyright".
*/
#include "FP_types.h"
-extern _float adf4(), ngf4();
-
-_float
+SINGLE
sbf4(s2,s1)
-_float s1,s2;
+SINGLE s1,s2;
{
- _float *result = &s1; /* s1 may not be in a register! */
+ SINGLE *result = &s1; /* s1 may not be in a register! */
- if (s2 == (_float) 0) {
+ if (s2 == (SINGLE) 0) {
return s1;
}
s2 = ngf4(s2);
#include "FP_types.h"
-extern _double adf8(), ngf8();
-
-_double
+DOUBLE
sbf8(s2,s1)
-_double s1,s2;
+DOUBLE s1,s2;
{
- _double *result = &s1; /* s1 may not be in a register! */
+ DOUBLE *result = &s1; /* s1 may not be in a register! */
- if (s2.__double[0] == 0 && s2.__double[1] == 0) {
+ if (s2.d[0] == 0 && s2.d[1] == 0) {
return s1;
}
s2 = ngf8(s2);
#include "FP_types.h"
+void
sft_ext(e1,e2)
EXTEND *e1,*e2;
{
s = e2;
s->exp += diff;
- b64_sft(&(s->m1), diff);
+ b64_sft(&(s->mantissa), diff);
}
/* $Header$ */
-# include "adder.h"
+# include "FP_types.h"
+void
b64_sft(e1,n)
B64 *e1;
int n;
}
}
+void
b64_lsft(e1)
B64 *e1;
{
e1->l_32 <<= 1;
}
+void
b64_rsft(e1)
B64 *e1;
{
#include "FP_types.h"
+void
sub_ext(e1,e2)
EXTEND *e1,*e2;
{
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 */
+ if (b64_add(&e1->mantissa,&e2->mantissa)) { /* addition carry */
+ b64_rsft(&e1->mantissa); /* shift mantissa one bit RIGHT */
e1->m1 |= 0x80000000L; /* set max bit */
e1->exp++; /* increase the exponent */
}
return a zero float (ZRF 4)
*/
+#include "FP_types.h"
+
+void
zrf4(l)
-long *l;
+SINGLE *l;
{
*l = 0L;
}
#include "FP_types.h"
+void
zrf8(z)
-_double *z;
+DOUBLE *z;
{
- z->__double[0] = 0L;
- z->__double[1] = 0L;
+ z->d[0] = 0L;
+ z->d[1] = 0L;
}
#include "FP_types.h"
+void
zrf_ext(e)
EXTEND *e;
{
public / ack.git / commitdiff