From c20a2155fafa59eb57c4a26e1ee0ff13e94ab7ec Mon Sep 17 00:00:00 2001 From: ceriel Date: Mon, 10 Jul 1989 11:17:19 +0000 Subject: [PATCH] Initial revision --- modules/src/flt_arith/.distr | 18 ++ modules/src/flt_arith/Makefile | 68 +++++ modules/src/flt_arith/b64_add.c | 30 ++ modules/src/flt_arith/b64_sft.c | 75 +++++ modules/src/flt_arith/flt_add.c | 82 ++++++ modules/src/flt_arith/flt_ar2flt.c | 48 ++++ modules/src/flt_arith/flt_arith.3 | 206 ++++++++++++++ modules/src/flt_arith/flt_arith.h | 30 ++ modules/src/flt_arith/flt_chk.c | 26 ++ modules/src/flt_arith/flt_cmp.c | 26 ++ modules/src/flt_arith/flt_div.c | 131 +++++++++ modules/src/flt_arith/flt_flt2ar.c | 55 ++++ modules/src/flt_arith/flt_modf.c | 32 +++ modules/src/flt_arith/flt_mul.c | 83 ++++++ modules/src/flt_arith/flt_nrm.c | 34 +++ modules/src/flt_arith/flt_str2fl.c | 443 +++++++++++++++++++++++++++++ modules/src/flt_arith/misc.h | 24 ++ modules/src/flt_arith/ucmp.c | 21 ++ 18 files changed, 1432 insertions(+) create mode 100644 modules/src/flt_arith/.distr create mode 100644 modules/src/flt_arith/Makefile create mode 100644 modules/src/flt_arith/b64_add.c create mode 100644 modules/src/flt_arith/b64_sft.c create mode 100644 modules/src/flt_arith/flt_add.c create mode 100644 modules/src/flt_arith/flt_ar2flt.c create mode 100644 modules/src/flt_arith/flt_arith.3 create mode 100644 modules/src/flt_arith/flt_arith.h create mode 100644 modules/src/flt_arith/flt_chk.c create mode 100644 modules/src/flt_arith/flt_cmp.c create mode 100644 modules/src/flt_arith/flt_div.c create mode 100644 modules/src/flt_arith/flt_flt2ar.c create mode 100644 modules/src/flt_arith/flt_modf.c create mode 100644 modules/src/flt_arith/flt_mul.c create mode 100644 modules/src/flt_arith/flt_nrm.c create mode 100644 modules/src/flt_arith/flt_str2fl.c create mode 100644 modules/src/flt_arith/misc.h create mode 100644 modules/src/flt_arith/ucmp.c diff --git a/modules/src/flt_arith/.distr b/modules/src/flt_arith/.distr new file mode 100644 index 000000000..4cd9e5206 --- /dev/null +++ b/modules/src/flt_arith/.distr @@ -0,0 +1,18 @@ +Makefile +b64_add.c +b64_sft.c +flt_add.c +flt_arith.h +flt_arith2flt.c +flt_chk.c +flt_cmp.c +flt_div.c +flt_flt2arith.c +flt_modf.c +flt_mul.c +flt_nrm.c +flt_str2flt.c +misc.h +ucmp.c +flt_arith.3 + diff --git a/modules/src/flt_arith/Makefile b/modules/src/flt_arith/Makefile new file mode 100644 index 000000000..ae65e3c50 --- /dev/null +++ b/modules/src/flt_arith/Makefile @@ -0,0 +1,68 @@ +# $Header$ +EMHOME = ../../.. +MODDIR=$(EMHOME)/modules +INSTALL = $(MODDIR)/install +COMPARE = $(MODDIR)/compare +INCLUDES = -I. -I$(MODDIR)/h +CFLAGS = -O $(INCLUDES) $(COPT) +AR = ar +SUF = o +LIBSUF = a + +LIBFLT = libflt.$(LIBSUF) + +SRC = b64_add.c flt_arith2flt.c flt_div.c flt_nrm.c b64_sft.c flt_chk.c \ + flt_flt2arith.c flt_str2flt.c flt_add.c flt_cmp.c flt_mul.c ucmp.c \ + flt_modf.c +OBJ = b64_add.$(SUF) flt_arith2flt.$(SUF) flt_div.$(SUF) flt_nrm.$(SUF) \ + b64_sft.$(SUF) flt_chk.$(SUF) flt_flt2arith.$(SUF) flt_str2flt.$(SUF) \ + flt_add.$(SUF) flt_cmp.$(SUF) flt_mul.$(SUF) ucmp.$(SUF) \ + flt_modf.$(SUF) + +.SUFFIXES: .$(SUF) +.c.$(SUF): + $(CC) -c $(CFLAGS) $*.c + +all: $(LIBFLT) + +$(LIBFLT): $(OBJ) + rm -f $(LIBFLT) + $(AR) r $(LIBFLT) $(OBJ) + -sh -c 'ranlib $(LIBFLT)' + +install: all + $(INSTALL) lib/$(LIBFLT) + $(INSTALL) h/flt_arith.h + $(INSTALL) man/flt_arith.3 + +cmp: all + -$(COMPARE) lib/$(LIBFLT) + -$(COMPARE) h/flt_arith.h + -$(COMPARE) man/flt_arith.3 + +pr: + @pr Makefile $(SRC) + +opr: + make pr | opr + +clean: + rm -f *.$(SUF) $(LIBFLT) + +lintlib: + lint $(INCLUDES) -Cflt $(SRC) + mv llib-lflt.ln $(MODDIR)/lib + +b64_add.$(SUF): misc.h flt_arith.h +flt_arith2flt.$(SUF): misc.h flt_arith.h +flt_div.$(SUF): misc.h flt_arith.h +flt_nrm.$(SUF): misc.h flt_arith.h +b64_sft.$(SUF): misc.h flt_arith.h +flt_chk.$(SUF): misc.h flt_arith.h +flt_flt2arith.$(SUF): misc.h flt_arith.h +flt_str2flt.$(SUF): misc.h flt_arith.h +flt_add.$(SUF): misc.h flt_arith.h +flt_cmp.$(SUF): misc.h flt_arith.h +flt_mul.$(SUF): misc.h flt_arith.h +flt_modf.$(SUF): misc.h flt_arith.h +ucmp.$(SUF): misc.h flt_arith.h diff --git a/modules/src/flt_arith/b64_add.c b/modules/src/flt_arith/b64_add.c new file mode 100644 index 000000000..c1b2ca817 --- /dev/null +++ b/modules/src/flt_arith/b64_add.c @@ -0,0 +1,30 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +# include "misc.h" + +int +b64_add(e1,e2) + register struct _mantissa *e1,*e2; +{ + int overflow; + int carry; + + /* add higher pair of 32 bits */ + overflow = ucmp(0xFFFFFFFF - e1->flt_h_32, e2->flt_h_32) < 0; + e1->flt_h_32 += e2->flt_h_32; + + /* add lower pair of 32 bits */ + carry = ucmp(0xFFFFFFFF - e1->flt_l_32, e2->flt_l_32) < 0; + e1->flt_l_32 += e2->flt_l_32; + + if ((carry) && ((++e1->flt_h_32 &~0xFFFFFFFF) || e1->flt_h_32 == 0)) { + e1->flt_h_32 = 0; + return(1); /* had a 64 bit overflow */ + } + return(overflow); /* return status from higher add */ +} diff --git a/modules/src/flt_arith/b64_sft.c b/modules/src/flt_arith/b64_sft.c new file mode 100644 index 000000000..e66c85acc --- /dev/null +++ b/modules/src/flt_arith/b64_sft.c @@ -0,0 +1,75 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +b64_sft(e,n) + register struct _mantissa *e; + register int n; +{ + if (n > 0) { + if (n > 63) { + e->flt_l_32 = 0; + e->flt_h_32 = 0; + return; + } + if (n >= 32) { + e->flt_l_32 = e->flt_h_32; + e->flt_h_32 = 0; + n -= 32; + } + if (n > 0) { + e->flt_l_32 = (e->flt_l_32 >> 1) & 0x7FFFFFFF; + e->flt_l_32 >>= (n - 1); + if (e->flt_h_32 != 0) { + e->flt_l_32 |= (e->flt_h_32 << (32 - n)) & 0xFFFFFFFF; + e->flt_h_32 = (e->flt_h_32 >> 1) & 0x7FFFFFFF; + e->flt_h_32 >>= (n - 1); + } + } + return; + } + n = -n; + if (n > 0) { + if (n > 63) { + e->flt_l_32 = 0; + e->flt_h_32 = 0; + return; + } + if (n >= 32) { + e->flt_h_32 = e->flt_l_32; + e->flt_l_32 = 0; + n -= 32; + } + if (n > 0) { + e->flt_h_32 = (e->flt_h_32 << n) & 0xFFFFFFFF; + if (e->flt_l_32 != 0) { + long l = (e->flt_l_32 >> 1) & 0x7FFFFFFF; + e->flt_h_32 |= (l >> (31 - n)); + e->flt_l_32 = (e->flt_l_32 << n) & 0xFFFFFFFF; + } + } + } +} + +b64_lsft(e) + register struct _mantissa *e; +{ + /* shift left 1 bit */ + e->flt_h_32 = (e->flt_h_32 << 1) & 0xFFFFFFFF; + if (e->flt_l_32 & 0x80000000L) e->flt_h_32 |= 1; + e->flt_l_32 = (e->flt_l_32 << 1) & 0xFFFFFFFF; +} + +b64_rsft(e) + register struct _mantissa *e; +{ + /* shift right 1 bit */ + e->flt_l_32 = (e->flt_l_32 >> 1) & 0x7FFFFFFF; + if (e->flt_h_32 & 1) e->flt_l_32 |= 0x80000000L; + e->flt_h_32 = (e->flt_h_32 >> 1) & 0x7FFFFFFF; +} diff --git a/modules/src/flt_arith/flt_add.c b/modules/src/flt_arith/flt_add.c new file mode 100644 index 000000000..4a38d63ef --- /dev/null +++ b/modules/src/flt_arith/flt_add.c @@ -0,0 +1,82 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +flt_add(e1,e2,e3) + register flt_arith *e1,*e2,*e3; +{ + /* Add two extended numbers e1 and e2, and put the result + in e3 + */ + flt_arith ce2; + int diff; + + flt_status = 0; + if ((e2->m1 | e2->m2) == 0L) { + *e3 = *e1; + return; + } + if ((e1->m1 | e1->m2) == 0L) { + *e3 = *e2; + return; + } + ce2 = *e2; + *e3 = *e1; + e1 = &ce2; + + /* adjust mantissas to equal power */ + diff = e3->flt_exp - e1->flt_exp; + if (diff < 0) { + diff = -diff; + e3->flt_exp += diff; + b64_sft(&(e3->flt_mantissa), diff); + } + else if (diff > 0) { + e1->flt_exp += diff; + b64_sft(&(e1->flt_mantissa), diff); + } + if (e1->flt_sign != e3->flt_sign) { + /* e3 + e1 = e3 - (-e1) */ + int tmp = ucmp(e1->m1, e3->m1); + int tmp2 = ucmp(e1->m2, e3->m2); + if (tmp > 0 || (tmp == 0 && tmp2 > 0)) { + /* abs(e1) > abs(e3) */ + if (tmp2 < 0) { + e1->m1 -= 1; /* carry in */ + } + e1->m1 -= e3->m1; + e1->m2 -= e3->m2; + *e3 = *e1; + } + else { + if (tmp2 > 0) + e3->m1 -= 1; /* carry in */ + e3->m1 -= e1->m1; + e3->m2 -= e1->m2; + } + } + else { + if (b64_add(&e3->flt_mantissa,&e1->flt_mantissa)) {/* addition carry */ + b64_sft(&e3->flt_mantissa,1);/* shift mantissa one bit RIGHT */ + e3->m1 |= 0x80000000L; /* set max bit */ + e3->flt_exp++; /* increase the exponent */ + } + } + flt_nrm(e3); + flt_chk(e3); +} + +flt_sub(e1,e2,e3) + flt_arith *e1,*e2,*e3; +{ + flt_arith ce2; + + ce2 = *e2; + ce2.flt_sign = ! ce2.flt_sign; + flt_add(e1,&ce2,e3); +} diff --git a/modules/src/flt_arith/flt_ar2flt.c b/modules/src/flt_arith/flt_ar2flt.c new file mode 100644 index 000000000..c13f97f43 --- /dev/null +++ b/modules/src/flt_arith/flt_ar2flt.c @@ -0,0 +1,48 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" +#include + +flt_arith2flt(n, e) + register arith n; + register flt_arith *e; +{ + /* Convert the arith "n" to a flt_arith "e". + */ + register int i; + + if (n == 0) { + e->flt_exp = 0; + e->flt_sign = 0; + e->m1 = e->m2 = 0; + return; + } + if (n < 0) { + e->flt_sign = 1; + n = -n; + } + else e->flt_sign = 0; + e->flt_exp = 63; + e->m1 = e->m2 = 0; + if (n < 0) { + n = 0x40000000; + while ((n << 1) > 0) n <<= 1; + e->flt_exp++; + } + for (i = 64; i > 0 && n != 0; i--) { + b64_rsft(&(e->flt_mantissa)); + e->m1 |= (n & 1) << 31; + n >>= 1; + } + + if (i > 0) { + b64_sft(&(e->flt_mantissa), i); + } + flt_status = 0; + flt_nrm(e); +} diff --git a/modules/src/flt_arith/flt_arith.3 b/modules/src/flt_arith/flt_arith.3 new file mode 100644 index 000000000..1307a5f60 --- /dev/null +++ b/modules/src/flt_arith/flt_arith.3 @@ -0,0 +1,206 @@ +.TH FLT_ARITH 3ACK " Fri June 30 1989" +.ad +.SH NAME +flt_arith \- high precision floating point arithmetic +.SH SYNOPSIS +.nf +.B #include +.PP +.B flt_add(e1, e2, e3) +.B flt_arith *e1, *e2, *e3; +.PP +.B flt_mul(e1, e2, e3) +.B flt_arith *e1, *e2, *e3; +.PP +.B flt_sub(e1, e2, e3) +.B flt_arith *e1, *e2, *e3; +.PP +.B flt_div(e1, e2, e3) +.B flt_arith *e1, *e2, *e3; +.PP +.B flt_modf(e1, intpart, fractpart) +.B flt_arith *e1, *intpart, *fractpart; +.PP +.B int flt_cmp(e1, e2) +.B flt_arith *e1, *e2; +.PP +.B int flt_str2flt(s, e) +.B char *s; +.B flt_arith *e; +.PP +.B flt_flt2str(e, buf, bufsize) +.B flt_arith *e; +.B char *buf; +.B int bufsize; +.PP +.B int flt_status; +.PP +.B #include +.B flt_arith2flt(n, e) +.B arith n; +.B flt_arith *e; +.PP +.B arith flt_flt2arith(e) +.B flt_arith *e; +.SH DESCRIPTION +This set of routines emulates floating point arithmetic, in a high +precision. It is intended primarily for compilers that need to evaluate +floating point expressions at compile-time. It could be argued that this +should be done in the floating point arithmetic of the target machine, +but EM does not define its floating point arithmetic. +.PP +.B flt_add +adds the numbers indicated by +.I e1 +and +.I e2 +and stores the result indirectly through +.IR e3 . +.PP +.B flt_mul +multiplies the numbers indicated by +.I e1 +and +.I e2 +and stores the result indirectly through +.IR e3 . +.PP +.B flt_sub +subtracts the number indicated by +.I e2 +from the one indicated by +.I e1 +and stores the result indirectly through +.IR e3 . +.PP +.B flt_div +divides the number indicated by +.I e1 +by the one indicated by +.I e2 +and stores the result indirectly through +.IR e3 . +.PP +.B flt_modf +splits the number indicated by +.I e +in an integer and a fraction part, and stores the integer part through +.I intpart +and the fraction part through +.IR fractpart . +So, adding the numbers indicated by +.I intpart +and +.I fractpart +results (in the absence of rounding error) in the number +indicated by +.IR e . +Also, the absolute value of the number indicated by +.I intpart +is less than or equal to the absolute value of the number indicated by +.IR e . +The absolute value of the number indicated by +.I fractpart +is less than 1. +.PP +.B flt_cmp +compares the numbers indicated by +.I e1 +and +.I e2 +and returns -1 if +.I e1 +< +.IR e2 , +0 if +.I e1 += +.IR e2 , +and 1 if +.I e1 +> +.IR e2 . +.PP +.B flt_str2flt +converts the string indicated by +.I s +to a floating point number, and stores this number through +.IR e. +The string should contain a floating point constant, which consists of +an integer part, a decimal point, a fraction part, an \f(CWe\fP or an +\f(CWE\fP, and an optionally signed integer exponent. The integer and +fraction parts both consist of a sequence of digits. They may not both be +missing. The decimal point, the \f(CWe\fP and the exponent may be +missing. +.PP +.B flt_flt2str +converts the number indicated by +.I e +into a string, in a scientific notation acceptable for EM. The result is +stored in +.IR buf . +At most +.I bufsize +characters are stored. +.PP +.B flt_arith2flt +converts the number +.I n +to the floating point format use in this package and returns the result +in +.IR e . +.PP +.B flt_flt2arith +truncates the number indicated by +.I e +to the largest integer value smaller than or equal to the number indicated by +.IR e . +It returns this value. +.PP +Before each operation, the +.I flt_status +variable is reset to 0. After an operation, it can be checked for one +of the following values: +.IP FLT_OVFL +.br +an overflow occurred. The result is a large value with the correct sign. +This can occur with the routines +.IR flt_add , +.IR flt_sub , +.IR flt_div , +.IR flt_mul , +.IR flt_flt2arith , +and +.IR flt_str2flt . +.IP FLT_UNFL +.br +an underflow occurred. The result is 0. +This can occur with the routines +.IR flt_div , +.IR flt_mul , +.IR flt_sub , +.IR flt_add , +and +.IR flt_str2flt . +.IP FLT_DIV0 +.br +divide by 0. The result is a large value with the sign of the dividend. +This can only occur with the routine +.IR flt_div . +.IP FLT_NOFLT +.br +indicates that the string did not represent a floating point number. The +result is 0. +This can only occur with the routine +.IR flt_str2flt . +.IP FLT_BTSM +.br +indicates that the buffer is too small. The contents of the buffer is +undefined. This can only occur with the routine +.IR flt_flt2str . +.SH FILES +~em/modules/h/flt_arith.h +.br +~em/modules/h/em_arith.h +.br +~em/modules/lib/libflt.a diff --git a/modules/src/flt_arith/flt_arith.h b/modules/src/flt_arith/flt_arith.h new file mode 100644 index 000000000..f33dd615f --- /dev/null +++ b/modules/src/flt_arith/flt_arith.h @@ -0,0 +1,30 @@ +/* + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + */ +/* $Header$ */ + +#ifndef __FLT_INCLUDED__ +#define __FLT_INCLUDED__ + +struct _mantissa { + long flt_h_32; + long flt_l_32; +}; + +struct _EXTEND { + short flt_sign; + short flt_exp; + struct _mantissa flt_mantissa; +}; + +typedef struct _EXTEND flt_arith; + +extern int flt_status; +#define FLT_OVFL 001 +#define FLT_UNFL 002 +#define FLT_DIV0 004 +#define FLT_NOFLT 010 +#define FLT_BTSM 020 + +#endif __FLT_INCLUDED__ diff --git a/modules/src/flt_arith/flt_chk.c b/modules/src/flt_arith/flt_chk.c new file mode 100644 index 000000000..2a027d014 --- /dev/null +++ b/modules/src/flt_arith/flt_chk.c @@ -0,0 +1,26 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +int flt_status = 0; + +flt_chk(e) + register flt_arith *e; +{ + if (e->flt_exp >= EXT_MAX) { + flt_status = FLT_OVFL; + e->flt_exp = EXT_MAX; + e->m1 = e->m2 = 0; + } + if (e->flt_exp <= EXT_MIN) { + flt_status = FLT_UNFL; + e->flt_exp = 0; + e->m1 = e->m2 = 0; + e->flt_sign = 0; + } +} diff --git a/modules/src/flt_arith/flt_cmp.c b/modules/src/flt_arith/flt_cmp.c new file mode 100644 index 000000000..fd21a9b40 --- /dev/null +++ b/modules/src/flt_arith/flt_cmp.c @@ -0,0 +1,26 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +int +flt_cmp(e1, e2) + register flt_arith *e1, *e2; +{ + int sign = e1->flt_sign ? -1 : 1; + int tmp; + + if (e1->flt_sign > e2->flt_sign) return -1; + if (e1->flt_sign < e2->flt_sign) return 1; + if (e1->flt_exp < e2->flt_exp) return -sign; + if (e1->flt_exp > e2->flt_exp) return sign; + if ((tmp = ucmp(e1->m1, e2->m1)) < 0) return -sign; + if (tmp > 0) return sign; + if ((tmp = ucmp(e1->m2, e2->m2)) < 0) return -sign; + if (tmp > 0) return sign; + return 0; +} diff --git a/modules/src/flt_arith/flt_div.c b/modules/src/flt_arith/flt_div.c new file mode 100644 index 000000000..72d51613d --- /dev/null +++ b/modules/src/flt_arith/flt_div.c @@ -0,0 +1,131 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +flt_div(e1,e2,e3) + register flt_arith *e1,*e2,*e3; +{ + short error = 0; + long result[2]; + register long *lp; + unsigned short u[9], v[5]; + register int j; + register unsigned short *u_p = u; + int maxv = 4; + + flt_status = 0; + e3->flt_sign = e1->flt_sign ^ e2->flt_sign; + if ((e2->m1 | e2->m2) == 0) { + flt_status = FLT_DIV0; + e3->flt_exp = EXT_MAX; + e3->m1 = e3->m2 = 0L; + return; + } + if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */ + e3->flt_exp = 0; /* make sure */ + e3->m1 = e3->m2 = 0L; + e3->flt_sign = 0; + return; + } + e3->flt_exp = e1->flt_exp - e2->flt_exp + 2; + + u[4] = (e1->m2 & 1) << 15; + b64_rsft(&(e1->flt_mantissa)); + u[0] = (e1->m1 >> 16) & 0xFFFF; + u[1] = e1->m1 & 0xFFFF; + u[2] = (e1->m2 >> 16) & 0xFFFF; + u[3] = e1->m2 & 0xFFFF; + u[5] = 0; u[6] = 0; u[7] = 0; + v[1] = (e2->m1 >> 16) & 0xFFFF; + v[2] = e2->m1 & 0xFFFF; + v[3] = (e2->m2 >> 16) & 0xFFFF; + v[4] = e2->m2 & 0xFFFF; + while (! v[maxv]) maxv--; + result[0] = 0; + result[1] = 0; + lp = result; + + /* + * Use an algorithm of Knuth (The art of programming, Seminumerical + * algorithms), to divide u by v. u and v are both seen as numbers + * with base 65536. + */ + for (j = 0; j <= 3; j++, u_p++) { + long q_est, temp; + + if (j == 2) lp++; + if (u_p[0] == 0 && u_p[1] < v[1]) continue; + temp = ((long)u_p[0] << 16) + u_p[1]; + if (u_p[0] >= v[1]) { + q_est = 0x0000FFFF; + } + else if (v[1] == 1) { + q_est = temp; + } + else if (temp >= 0) { + q_est = temp / v[1]; + } + else { + long rem; + q_est = (0x7FFFFFFF/v[1])+((temp&0x7FFFFFFF)/v[1]); + rem = (0x7FFFFFFF%v[1])+((temp&0x7FFFFFFF)%v[1])+1; + while (rem > v[1]) { + q_est++; + rem -= v[1]; + } + } + temp -= q_est * v[1]; + while (!(temp&0xFFFF0000) && ucmp(v[2]*q_est,((temp<<16)+u_p[2])) > 0) { + q_est--; + temp += v[1]; + } + /* Now, according to Knuth, we have an estimate of the + quotient, that is either correct or one too big, but + almost always correct. + */ + if (q_est != 0) { + int i; + long k = 0; + int borrow = 0; + + for (i = maxv; i > 0; i--) { + long tmp = q_est * v[i] + k + borrow; + unsigned short md = tmp & 0xFFFF; + + borrow = (md > u_p[i]); + u_p[i] -= md; + k = (tmp >> 16) & 0xFFFF; + } + k += borrow; + borrow = u_p[0] < k; + u_p[0] -= k; + + if (borrow) { + /* So, this does not happen often; the estimate + was one too big; correct this + */ + *lp |= (j & 1) ? (q_est - 1) : ((q_est-1)<<16); + borrow = 0; + for (i = maxv; i > 0; i--) { + long tmp + = v[i]+(long)u_p[i]+borrow; + + u_p[i] = tmp & 0xFFFF; + borrow = (tmp >> 16) & 0xFFFF; + } + u_p[0] += borrow; + } + else *lp |= (j & 1) ? q_est : (q_est<<16); + } + } + e3->m1 = result[0]; + e3->m2 = result[1]; + + flt_nrm(e3); + flt_chk(e3); +} diff --git a/modules/src/flt_arith/flt_flt2ar.c b/modules/src/flt_arith/flt_flt2ar.c new file mode 100644 index 000000000..4a91a0beb --- /dev/null +++ b/modules/src/flt_arith/flt_flt2ar.c @@ -0,0 +1,55 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" +#include + +arith +flt_flt2arith(e) + register flt_arith *e; +{ + /* Convert the flt_arith "n" to an arith. + */ + arith n; + struct _mantissa a; + + flt_status = 0; + if (e->flt_exp < 0) { + /* absolute value of result < 1. + Return value only depends on sign: + */ + return -e->flt_sign; + } + + if (e->flt_exp > (8*sizeof(arith)-2)) { + /* probably overflow, but there is one exception: + */ + if (e->flt_sign) { + n = 0x80; + while (n << 8) n <<= 8; + if (e->flt_exp == 8*sizeof(arith)-1 && + e->m2 == 0 && + e->m1 == 0x80000000) { + /* No overflow in this case */ + } + else flt_status = FLT_OVFL; + return n; + } + n = 0x7F; + while ((n << 8) > 0) { + n <<= 8; + n |= 0xFF; + } + return n; + } + + a = e->flt_mantissa; + b64_sft(&a, 63-e->flt_exp); + n = a.flt_l_32 | ((a.flt_h_32 << 16) << 16); + /* not << 32; this could be an undefined operation */ + return n; +} diff --git a/modules/src/flt_arith/flt_modf.c b/modules/src/flt_arith/flt_modf.c new file mode 100644 index 000000000..27071c4a7 --- /dev/null +++ b/modules/src/flt_arith/flt_modf.c @@ -0,0 +1,32 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +flt_modf(e, ipart, fpart) + register flt_arith *e, *ipart, *fpart; +{ + if (e->flt_exp < 0) { + *fpart = *e; + ipart->flt_sign = 0; + ipart->flt_exp = 0; + ipart->m1 = ipart->m2 = 0; + return; + } + if (e->flt_exp >= 63) { + fpart->flt_sign = 0; + fpart->flt_exp = 0; + fpart->m1 = fpart->m2 = 0; + *ipart = *e; + return; + } + *ipart = *e; + /* "loose" low order bits */ + b64_sft(&(ipart->flt_mantissa), 63 - e->flt_exp); + b64_sft(&(ipart->flt_mantissa), e->flt_exp - 63); + flt_sub(e, ipart, fpart); +} diff --git a/modules/src/flt_arith/flt_mul.c b/modules/src/flt_arith/flt_mul.c new file mode 100644 index 000000000..19b0aea52 --- /dev/null +++ b/modules/src/flt_arith/flt_mul.c @@ -0,0 +1,83 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +flt_mul(e1,e2,e3) + register flt_arith *e1,*e2,*e3; +{ + /* Multiply the extended numbers e1 and e2, and put the + result in e3. + */ + register int i,j; /* loop control */ + unsigned short mp[4]; + unsigned short mc[4]; + unsigned short result[8]; /* result */ + + register unsigned short *pres; + + flt_status = 0; + + /* first save the sign (XOR) */ + e3->flt_sign = e1->flt_sign ^ e2->flt_sign; + + /* compute new exponent */ + e3->flt_exp = e1->flt_exp + e2->flt_exp + 1; + + /* 128 bit multiply of mantissas */ + + /* assign unknown long formats */ + /* to known unsigned word formats */ + mp[0] = (e1->m1 >> 16) & 0xFFFF; + mp[1] = e1->m1 & 0xFFFF; + mp[2] = (e1->m2 >> 16) & 0xFFFF; + mp[3] = e1->m2 & 0xFFFF; + mc[0] = (e2->m1 >> 16) & 0xFFFF; + mc[1] = e2->m1 & 0xFFFF; + mc[2] = (e2->m2 >> 16) & 0xFFFF; + mc[3] = e2->m2 & 0xFFFF; + for (i = 8; i--;) { + result[i] = 0; + } + /* + * fill registers with their components + */ + for(i=4, pres = &result[4];i--;pres--) if (mp[i]) { + unsigned short k = 0; + long mpi = mp[i]; + for(j=4;j--;) { + long tmp = (long)pres[j] + k; + if (mc[j]) tmp += mpi * mc[j]; + pres[j] = tmp & 0xFFFF; + k = (tmp >> 16) & 0xFFFF; + } + pres[-1] = k; + } + + if (! (result[0] & 0x8000)) { + e3->flt_exp--; + for (i = 0; i <= 3; i++) { + result[i] <<= 1; + if (result[i+1]&0x8000) result[i] |= 1; + } + result[4] <<= 1; + } + /* + * combine the registers to a total + */ + e3->m1 = ((long)result[0] << 16) + result[1]; + e3->m2 = ((long)result[2] << 16) + result[3]; + if (result[4] & 0x8000) { + if (++e3->m2 == 0) { + if (++e3->m1 == 0) { + e3->m1 = 0x80000000; + e3->flt_exp++; + } + } + } + flt_chk(e3); +} diff --git a/modules/src/flt_arith/flt_nrm.c b/modules/src/flt_arith/flt_nrm.c new file mode 100644 index 000000000..166263076 --- /dev/null +++ b/modules/src/flt_arith/flt_nrm.c @@ -0,0 +1,34 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +flt_nrm(e) + register flt_arith *e; +{ + if ((e->m1 | e->m2) == 0L) + return; + + /* if top word is zero mov low word */ + /* to top word, adjust exponent value */ + if (e->m1 == 0L) { + e->m1 = e->m2; + e->m2 = 0L; + e->flt_exp -= 32; + } + if ((e->m1 & 0x80000000) == 0) { + long l = 0x40000000; + int cnt = -1; + + while (! (l & e->m1)) { + l >>= 1; + cnt--; + } + e->flt_exp += cnt; + b64_sft(&(e->flt_mantissa), cnt); + } +} diff --git a/modules/src/flt_arith/flt_str2fl.c b/modules/src/flt_arith/flt_str2fl.c new file mode 100644 index 000000000..63662322f --- /dev/null +++ b/modules/src/flt_arith/flt_str2fl.c @@ -0,0 +1,443 @@ +/* + (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +/* The following tables can be computed with the following bc(1) + program: + +obase=16 +scale=0 +define t(x){ + auto a, b, c + a=2;b=1;c=2^32;n=1 + while(a= sizeof(r_big_10pow)/sizeof(r_big_10pow[0])) { + flt_mul(&x, &r_big_10pow[sizeof(r_big_10pow)/sizeof(r_big_10pow[0])-1],&x); + divsz -= sizeof(r_big_10pow)/sizeof(r_big_10pow[0])-1; + flt_chk(e); + if (flt_status != 0) return; + } + flt_mul(&x, &r_big_10pow[divsz], e); + } + else { + flt_mul(e, &s10pow[modsz], &x); + while (divsz >= sizeof(big_10pow)/sizeof(big_10pow[0])) { + flt_mul(&x, &big_10pow[sizeof(big_10pow)/sizeof(big_10pow[0])-1],&x); + divsz -= sizeof(big_10pow)/sizeof(big_10pow[0])-1; + flt_chk(e); + if (flt_status != 0) return; + } + flt_mul(&x, &big_10pow[divsz], e); + } +} + +flt_str2flt(s, e) + register char *s; + register flt_arith *e; +{ + register int c; + int dotseen = 0; + int digitseen = 0; + int exp = 0; + + while (isspace(*s)) s++; + + flt_status = 0; + + e->flt_sign = 0; + e->flt_exp = 0; + e->m1 = e->m2 = 0; + + c = *s; + switch(c) { + case '-': + e->flt_sign = 1; + case '+': + s++; + } + while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) { + if (c == '.') continue; + digitseen = 1; + if (e->m1 <= 0x7FFFFFFF/5) { + struct _mantissa a1; + + a1 = e->flt_mantissa; + b64_sft(&(e->flt_mantissa), -3); + b64_sft(&a1, -1); + b64_add(&(e->flt_mantissa), &a1); + a1.flt_h_32 = 0; + a1.flt_l_32 = c - '0'; + b64_add(&(e->flt_mantissa), &a1); + } + else exp++; + if (dotseen) exp--; + } + if (! digitseen) { + flt_status = FLT_NOFLT; + return; + } + + if (c == 'E' || c == 'e') { + int exp1 = 0; + int sign = 1; + + switch(*s) { + case '-': + sign = -1; + case '+': + s++; + } + if (c = *s, isdigit(c)) { + do { + exp1 = 10 * exp1 + (c - '0'); + } while (c = *++s, isdigit(c)); + } + exp += sign * exp1; + } + if (e->m1 == 0 && e->m2 == 0) return; + e->flt_exp = 63; + flt_nrm(e); + add_exponent(e, exp); + flt_chk(e); +} + +#define NDIGITS 128 + +static char * +flt_ecvt(e, ndigit, decpt, sign) + register flt_arith *e; + int ndigit, *decpt, *sign; +{ + /* Like ecvt(), but for extended precision */ + + static char buf[NDIGITS+1]; + register char *p = buf; + register char *pe; + + if (ndigit < 0) ndigit = 0; + if (ndigit > NDIGITS) ndigit = NDIGITS; + pe = &buf[ndigit]; + buf[0] = '\0'; + + *sign = 0; + if (e->flt_sign) { + *sign = 1; + e->flt_sign = 0; + } + + *decpt = 0; + if (e->m1 != 0) { +#define BIGSZ (sizeof(big_10pow)/sizeof(big_10pow[0])) +#define SMALLSZ (sizeof(s10pow)/sizeof(s10pow[0])) + register flt_arith *pp = &big_10pow[1]; + + while (flt_cmp(e, &big_10pow[BIGSZ-1]) >= 0) { + flt_mul(e,&r_big_10pow[BIGSZ-1],e); + *decpt += (BIGSZ-1)*SMALLSZ; + } + while (flt_cmp(e,pp) >= 0) { + pp++; + } + pp--; + flt_mul(e,&r_big_10pow[pp-big_10pow],e); + *decpt += (pp - big_10pow)*SMALLSZ; + pp = &s10pow[1]; + while (pp < &s10pow[SMALLSZ] && cmp_ext(e, pp) >= 0) pp++; + pp--; + flt_mul(e, &r_10pow[pp-s10pow], e); + *decpt += pp - s10pow; + + if (cmp_ext(e, &s10pow[0]) < 0) { + while (flt_cmp(e, &r_big_10pow[BIGSZ-1]) < 0) { + flt_mul(e,&big_10pow[BIGSZ-1],e); + *decpt -= (BIGSZ-1)*SMALLSZ; + } + pp = &r_big_10pow[1]; + while(cmp_ext(e,pp) < 0) pp++; + pp--; + flt_mul(e,&big_10pow[pp-r_big_10pow],e); + *decpt -= (pp - r_big_10pow)*SMALLSZ; + /* here, value >= 10 ** -28 */ + flt_mul(e, &s10pow[1], e); + (*decpt)--; + pp = &r_10pow[0]; + while(cmp_ext(e, pp) < 0) pp++; + flt_mul(e, &s10pow[pp-r_10pow], e); + *decpt -= pp - r_10pow; + } + (*decpt)++; /* because now value in [1.0, 10.0) */ + } + while (p <= pe) { + if (e->flt_exp >= 0) { + flt_arith x; + + x.m2 = 0; x.flt_exp = e->flt_exp; + x.flt_sign = 1; + x.m1 = e->m1>>(31-e->flt_exp); + *p++ = (x.m1) + '0'; + x.m1 = x.m1 << (31-e->flt_exp); + add_ext(e, &x, e); + } + else *p++ = '0'; + flt_mul(e, &s10pow[1], e); + } + if (pe >= buf) { + p = pe; + *p += 5; /* round of at the end */ + while (*p > '9') { + *p = '0'; + if (p > buf) ++*--p; + else { + *p = '1'; + ++*decpt; + } + } + *pe = '\0'; + } + return buf; +} + +flt_flt2str(e, buf, bufsize) + flt_arith *e; + char *buf; +{ + +#define NDIG 25 + int sign, dp; + register int i; + register char *s1; + char Xbuf[NDIG+12]; + register char *s = Xbuf; + + flt_status = 0; + s1 = flt_ecvt(e,NDIG,&dp,&sign); + if (sign) + *s++ = '-'; + *s++ = *s1++; + *s++ = '.'; + for (i = NDIG-1; i > 0; i--) { + if (*s1) *s++ = *s1++; + else *s++ = '0'; + } + if ((e->m1 | e->m2) || e->flt_exp == EXT_MIN || e->flt_exp == EXT_MAX) { + --dp ; + } + if (dp == 0) { + *s = 0; + return; + } + *s++ = 'e'; + if ( dp<0 ) { + *s++ = '-' ; dp= -dp ; + } else { + *s++ = '+' ; + } + s1 = &Xbuf[NDIG+12]; + *--s1 = '\0'; + do { + *--s1 = dp % 10 + '0'; + dp /= 10; + } while (dp != 0); + while (*s1) *s++ = *s1++; + *s++ = '\0'; + if (s - Xbuf > bufsize) { + flt_status = FLT_BTSM; + return; + } + s = Xbuf; + s1 = buf; + do { + *s1++ = *s++; + } while (*s); +} diff --git a/modules/src/flt_arith/misc.h b/modules/src/flt_arith/misc.h new file mode 100644 index 000000000..bb496d5b0 --- /dev/null +++ b/modules/src/flt_arith/misc.h @@ -0,0 +1,24 @@ +/* + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + */ +/* $Header$ */ + +#include + +/* some short-hands ... */ +#define m1 flt_mantissa.flt_h_32 +#define m2 flt_mantissa.flt_l_32 + +/* some constants */ +#define EXT_MAX 16384 /* max exponent */ +#define EXT_MIN (-16384) /* min exponent */ + +/* hiding of names: */ +#define ucmp flt__ucmp +#define flt_nrm flt__nrm +#define flt_chk flt__chk +#define b64_add flt__64add +#define b64_sft flt__64sft +#define b64_rsft flt__64rsft +#define b64_lsft flt__64lsft diff --git a/modules/src/flt_arith/ucmp.c b/modules/src/flt_arith/ucmp.c new file mode 100644 index 000000000..016c9ad16 --- /dev/null +++ b/modules/src/flt_arith/ucmp.c @@ -0,0 +1,21 @@ +/* + (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. + See the copyright notice in the ACK home directory, in the file "Copyright". +*/ + +/* $Header$ */ + +#include "misc.h" + +int +ucmp(l1,l2) + long l1,l2; +{ + if (l1 == l2) return 0; + if (l2 >= 0) { + if (l1 > l2) return 1; + return -1; + } + if (l1 > 0 || l1 < l2) return -1; + return 1; +}