Updated for ANSI C, cleaned up a bit

This commit is contained in:
ceriel 1993-01-05 12:06:58 +00:00
parent 55131b091f
commit 005f32298f
38 changed files with 295 additions and 262 deletions

View file

@ -13,7 +13,7 @@
/********************************************************/
/*
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
@ -22,48 +22,92 @@
*/
/********************************************************/
#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

View file

@ -11,6 +11,7 @@
#include "FP_types.h"
void
add_ext(e1,e2)
register EXTEND *e1,*e2;
{
@ -45,8 +46,8 @@ 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 */
}

View file

@ -13,7 +13,7 @@
# include <stdio.h>
# endif
# include "adder.h"
# include "FP_types.h"
# define UNKNOWN -1
# define TRUE 1
# define FALSE 0
@ -22,6 +22,7 @@
/*
* add 64 bits
*/
int
b64_add(e1,e2)
/*
* pointers to 64 bit 'registers'
@ -45,6 +46,5 @@ register B64 *e1,*e2;
# 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 */
}

View file

@ -11,23 +11,23 @@
#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;
}

View file

@ -11,23 +11,23 @@
#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;
}

View file

@ -6,7 +6,7 @@
/* $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.
@ -17,13 +17,13 @@
#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]);
}

View file

@ -6,7 +6,7 @@
/* $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.
@ -17,13 +17,13 @@
#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);
}

View file

@ -21,15 +21,15 @@ long
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 */
@ -47,6 +47,6 @@ _double src; /* assume worst case */
if (buf.sign)
new = -new;
done:
src.__double[ss == 8] = new;
src.d[ss == 8] = new;
return(new);
}

View file

@ -20,15 +20,15 @@ long
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;
@ -38,6 +38,6 @@ _double src; /* assume worst case */
}
new = buf.m1 >> (31-buf.exp);
done:
src.__double[ss == 8] = new;
src.d[ss == 8] = new;
return(new);
}

View file

@ -6,7 +6,7 @@
/* $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
@ -16,7 +16,7 @@
#include "FP_types.h"
_float
SINGLE
cif4(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
@ -24,22 +24,22 @@ 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 */
@ -52,6 +52,6 @@ long src; /* largest possible integer to convert */
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);
}

View file

@ -16,17 +16,17 @@
#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;
@ -51,6 +51,6 @@ long src; /* largest possible integer to convert */
if (ss != sizeof(long))
buf.m1 <<= 16;
nrm_ext(&buf);
compact(&buf,result,8);
compact(&buf,&result->d[0],8);
return(*result);
}

View file

@ -14,7 +14,7 @@
int
cmf4(f1,f2)
_float f1,f2;
SINGLE f1,f2;
{
/*
* return ((f1 < f2) ? 1 : (f1 - f2))

View file

@ -12,8 +12,9 @@
#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)
/*

View file

@ -15,21 +15,22 @@
# 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 */
@ -42,25 +43,24 @@ dbl_over: trap(EFOVFL);
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 */
@ -72,17 +72,17 @@ dbl_over: trap(EFOVFL);
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++;
}
}
@ -101,24 +101,24 @@ dbl_over: trap(EFOVFL);
* 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
}
@ -129,9 +129,9 @@ dbl_over: trap(EFOVFL);
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 */
@ -144,16 +144,16 @@ sgl_over: trap(EFOVFL);
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 */
@ -165,13 +165,13 @@ sgl_over: trap(EFOVFL);
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)
@ -188,16 +188,15 @@ sgl_over: trap(EFOVFL);
* 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);
}
}

View file

@ -6,7 +6,7 @@
/* $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
@ -16,31 +16,31 @@
#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 */
@ -53,6 +53,6 @@ long src; /* largest possible integer to convert */
/* adjust mantissa field */
nrm_ext(&buf);
compact(&buf,(_double *) result,4);
compact(&buf,result,4);
return *result;
}

View file

@ -16,7 +16,7 @@
#include "FP_types.h"
_double
DOUBLE
cuf8(ss,src)
int ss; /* source size */
long src; /* largest possible integer to convert */
@ -36,7 +36,7 @@ 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 */
@ -50,6 +50,6 @@ long src; /* largest possible integer to convert */
/* 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);
}

View file

@ -26,19 +26,20 @@
*/
/********************************************************/
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) {
@ -60,8 +61,8 @@ EXTEND *e1,*e2;
* 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
@ -75,8 +76,8 @@ EXTEND *e1,*e2;
/* init control variables */
count = 64;
result[0] = 0L;
result[1] = 0L;
result.h_32 = 0L;
result.l_32 = 0L;
/* partial product division loop */
@ -84,7 +85,7 @@ EXTEND *e1,*e2;
/* 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 */
@ -95,7 +96,7 @@ EXTEND *e1,*e2;
; /* 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 */
@ -116,7 +117,7 @@ EXTEND *e1,*e2;
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
@ -124,7 +125,7 @@ EXTEND *e1,*e2;
} /* 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;
@ -132,16 +133,16 @@ EXTEND *e1,*e2;
}
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;
@ -152,9 +153,9 @@ EXTEND *e1,*e2;
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
@ -241,8 +242,7 @@ EXTEND *e1,*e2;
*/
INEXACT();
#endif
e1->m1 = result[0];
e1->m2 = result[1];
e1->mantissa = result;
nrm_ext(e1);
if (e1->exp < EXT_MIN) {

View file

@ -6,22 +6,22 @@
/* $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;
}

View file

@ -11,17 +11,17 @@
#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;
}

View file

@ -32,8 +32,9 @@
#include "get_put.h"
/********************************************************/
void
extend(from,to,size)
_double *from;
unsigned long *from;
EXTEND *to;
int size;
{

View file

@ -11,13 +11,9 @@
#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;
@ -25,7 +21,7 @@ struct fef4_returns *r;
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;
}
@ -33,5 +29,5 @@ struct fef4_returns *r;
p->e = buf.exp+1;
buf.exp = -1;
}
compact(&buf,(_double *) &p->f,sizeof(_float));
compact(&buf,&p->f,sizeof(SINGLE));
}

View file

@ -11,13 +11,9 @@
#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;
@ -25,7 +21,7 @@ struct fef8_returns *r;
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;
}
@ -33,5 +29,5 @@ struct fef8_returns *r;
p->e = buf.exp + 1;
buf.exp = -1;
}
compact(&buf,&p->f,sizeof(_double));
compact(&buf,&p->f.d[0],sizeof(DOUBLE));
}

View file

@ -12,26 +12,20 @@
#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;
@ -42,8 +36,8 @@ struct fif4_returns *p;
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);
}

View file

@ -12,40 +12,34 @@
#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);
}

View file

@ -11,16 +11,16 @@
#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);
}

View file

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

View file

@ -9,12 +9,12 @@
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;
{

View file

@ -14,13 +14,13 @@
#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;
}

View file

@ -15,13 +15,13 @@
#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;
}

View file

@ -14,6 +14,7 @@
#include "FP_shift.h"
#include "FP_types.h"
void
nrm_ext(e1)
EXTEND *e1;
{
@ -44,6 +45,6 @@ EXTEND *e1;
cnt--;
}
e1->exp += cnt;
b64_sft(&(e1->m1), cnt);
b64_sft(&(e1->mantissa), cnt);
}
}

View file

@ -1,5 +1,5 @@
/*
(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".
*/
@ -11,15 +11,13 @@
#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);

View file

@ -11,15 +11,13 @@
#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);

View file

@ -13,6 +13,7 @@
#include "FP_types.h"
void
sft_ext(e1,e2)
EXTEND *e1,*e2;
{
@ -34,5 +35,5 @@ EXTEND *e1,*e2;
s = e2;
s->exp += diff;
b64_sft(&(s->m1), diff);
b64_sft(&(s->mantissa), diff);
}

View file

@ -5,8 +5,9 @@
/* $Header$ */
# include "adder.h"
# include "FP_types.h"
void
b64_sft(e1,n)
B64 *e1;
int n;
@ -53,6 +54,7 @@ int n;
}
}
void
b64_lsft(e1)
B64 *e1;
{
@ -62,6 +64,7 @@ B64 *e1;
e1->l_32 <<= 1;
}
void
b64_rsft(e1)
B64 *e1;
{

View file

@ -11,6 +11,7 @@
#include "FP_types.h"
void
sub_ext(e1,e2)
EXTEND *e1,*e2;
{
@ -25,8 +26,8 @@ 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 */
}

View file

@ -9,8 +9,11 @@
return a zero float (ZRF 4)
*/
#include "FP_types.h"
void
zrf4(l)
long *l;
SINGLE *l;
{
*l = 0L;
}

View file

@ -11,10 +11,11 @@
#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;
}

View file

@ -11,6 +11,7 @@
#include "FP_types.h"
void
zrf_ext(e)
EXTEND *e;
{