Initial revision
This commit is contained in:
ceriel
parent
811612634a
commit
c20a2155fa
18
modules/src/flt_arith/.distr
Normal file
18
modules/src/flt_arith/.distr
Normal file
|
|
@ -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
|
||||||
|
|
||||||
68
modules/src/flt_arith/Makefile
Normal file
68
modules/src/flt_arith/Makefile
Normal file
|
|
@ -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
|
||||||
30
modules/src/flt_arith/b64_add.c
Normal file
30
modules/src/flt_arith/b64_add.c
Normal file
|
|
@ -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 */
|
||||||
|
}
|
||||||
75
modules/src/flt_arith/b64_sft.c
Normal file
75
modules/src/flt_arith/b64_sft.c
Normal file
|
|
@ -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;
|
||||||
|
}
|
||||||
82
modules/src/flt_arith/flt_add.c
Normal file
82
modules/src/flt_arith/flt_add.c
Normal file
|
|
@ -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);
|
||||||
|
}
|
||||||
48
modules/src/flt_arith/flt_ar2flt.c
Normal file
48
modules/src/flt_arith/flt_ar2flt.c
Normal file
|
|
@ -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 <em_arith.h>
|
||||||
|
|
||||||
|
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);
|
||||||
|
}
|
||||||
206
modules/src/flt_arith/flt_arith.3
Normal file
206
modules/src/flt_arith/flt_arith.3
Normal file
|
|
@ -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 <flt_arith.h>
|
||||||
|
.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 <em_arith.h>
|
||||||
|
.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
|
||||||
30
modules/src/flt_arith/flt_arith.h
Normal file
30
modules/src/flt_arith/flt_arith.h
Normal file
|
|
@ -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__
|
||||||
26
modules/src/flt_arith/flt_chk.c
Normal file
26
modules/src/flt_arith/flt_chk.c
Normal file
|
|
@ -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;
|
||||||
|
}
|
||||||
|
}
|
||||||
26
modules/src/flt_arith/flt_cmp.c
Normal file
26
modules/src/flt_arith/flt_cmp.c
Normal file
|
|
@ -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;
|
||||||
|
}
|
||||||
131
modules/src/flt_arith/flt_div.c
Normal file
131
modules/src/flt_arith/flt_div.c
Normal file
|
|
@ -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);
|
||||||
|
}
|
||||||
55
modules/src/flt_arith/flt_flt2ar.c
Normal file
55
modules/src/flt_arith/flt_flt2ar.c
Normal file
|
|
@ -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 <em_arith.h>
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
32
modules/src/flt_arith/flt_modf.c
Normal file
32
modules/src/flt_arith/flt_modf.c
Normal file
|
|
@ -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);
|
||||||
|
}
|
||||||
83
modules/src/flt_arith/flt_mul.c
Normal file
83
modules/src/flt_arith/flt_mul.c
Normal file
|
|
@ -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);
|
||||||
|
}
|
||||||
34
modules/src/flt_arith/flt_nrm.c
Normal file
34
modules/src/flt_arith/flt_nrm.c
Normal file
|
|
@ -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);
|
||||||
|
}
|
||||||
|
}
|
||||||
443
modules/src/flt_arith/flt_str2fl.c
Normal file
443
modules/src/flt_arith/flt_str2fl.c
Normal file
|
|
@ -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<x) {
|
||||||
|
b=a;n+=n;a*=a
|
||||||
|
}
|
||||||
|
n/=2
|
||||||
|
a=b
|
||||||
|
while(b<x) {
|
||||||
|
a=b;b*=c;n+=32
|
||||||
|
}
|
||||||
|
n-=32
|
||||||
|
b=a
|
||||||
|
while(a<x) {
|
||||||
|
b=a;a+=a;n+=1
|
||||||
|
}
|
||||||
|
n-=1
|
||||||
|
x*=16^16
|
||||||
|
b=x%a
|
||||||
|
x/=a
|
||||||
|
if(a<=(2*b)) x+=1
|
||||||
|
obase=10
|
||||||
|
n
|
||||||
|
obase=16
|
||||||
|
return(x)
|
||||||
|
}
|
||||||
|
for (i=1;i<28;i++) {
|
||||||
|
t(10^i)
|
||||||
|
}
|
||||||
|
0
|
||||||
|
for (i=1;i<20;i++) {
|
||||||
|
t(10^(28*i))
|
||||||
|
}
|
||||||
|
0
|
||||||
|
define r(x){
|
||||||
|
auto a, b, c
|
||||||
|
a=2;b=1;c=2^32;n=1
|
||||||
|
while(a<x) {
|
||||||
|
b=a;n+=n;a*=a
|
||||||
|
}
|
||||||
|
n/=2
|
||||||
|
a=b
|
||||||
|
while(b<x) {
|
||||||
|
a=b;b*=c;n+=32
|
||||||
|
}
|
||||||
|
n-=32
|
||||||
|
b=a
|
||||||
|
while(a<x) {
|
||||||
|
b=a;a+=a;n+=1
|
||||||
|
}
|
||||||
|
a=b
|
||||||
|
a*=16^16
|
||||||
|
b=a%x
|
||||||
|
a/=x
|
||||||
|
if(x<=(2*b)) a+=1
|
||||||
|
obase=10
|
||||||
|
-n
|
||||||
|
obase=16
|
||||||
|
return(a)
|
||||||
|
}
|
||||||
|
for (i=1;i<28;i++) {
|
||||||
|
r(10^i)
|
||||||
|
}
|
||||||
|
0
|
||||||
|
for (i=1;i<20;i++) {
|
||||||
|
r(10^(28*i))
|
||||||
|
}
|
||||||
|
0
|
||||||
|
|
||||||
|
*/
|
||||||
|
static flt_arith s10pow[] = { /* representation of 10 ** i */
|
||||||
|
{ 0, 0, 0x80000000, 0 },
|
||||||
|
{ 0, 3, 0xA0000000, 0 },
|
||||||
|
{ 0, 6, 0xC8000000, 0 },
|
||||||
|
{ 0, 9, 0xFA000000, 0 },
|
||||||
|
{ 0, 13, 0x9C400000, 0 },
|
||||||
|
{ 0, 16, 0xC3500000, 0 },
|
||||||
|
{ 0, 19, 0xF4240000, 0 },
|
||||||
|
{ 0, 23, 0x98968000, 0 },
|
||||||
|
{ 0, 26, 0xBEBC2000, 0 },
|
||||||
|
{ 0, 29, 0xEE6B2800, 0 },
|
||||||
|
{ 0, 33, 0x9502F900, 0 },
|
||||||
|
{ 0, 36, 0xBA43B740, 0 },
|
||||||
|
{ 0, 39, 0xE8D4A510, 0 },
|
||||||
|
{ 0, 43, 0x9184E72A, 0 },
|
||||||
|
{ 0, 46, 0xB5E620F4, 0x80000000 },
|
||||||
|
{ 0, 49, 0xE35FA931, 0xA0000000 },
|
||||||
|
{ 0, 53, 0x8E1BC9BF, 0x04000000 },
|
||||||
|
{ 0, 56, 0xB1A2BC2E, 0xC5000000 },
|
||||||
|
{ 0, 59, 0xDE0B6B3A, 0x76400000 },
|
||||||
|
{ 0, 63, 0x8AC72304, 0x89E80000 },
|
||||||
|
{ 0, 66, 0xAD78EBC5, 0xAC620000 },
|
||||||
|
{ 0, 69, 0xD8D726B7, 0x177A8000 },
|
||||||
|
{ 0, 73, 0x87867832, 0x6EAC9000 },
|
||||||
|
{ 0, 76, 0xA968163F, 0x0A57B400 },
|
||||||
|
{ 0, 79, 0xD3C21BCE, 0xCCEDA100 },
|
||||||
|
{ 0, 83, 0x84595161, 0x401484A0 },
|
||||||
|
{ 0, 86, 0xA56FA5B9, 0x9019A5C8 },
|
||||||
|
{ 0, 89, 0xCECB8F27, 0xF4200F3A }
|
||||||
|
};
|
||||||
|
static flt_arith big_10pow[] = { /* representation of 10 ** (28*i) */
|
||||||
|
{ 0, 0, 0x80000000, 0 },
|
||||||
|
{ 0, 93, 0x813F3978, 0xF8940984 },
|
||||||
|
{ 0, 186, 0x82818F12, 0x81ED44A0 },
|
||||||
|
{ 0, 279, 0x83C7088E, 0x1AAB65DB },
|
||||||
|
{ 0, 372, 0x850FADC0, 0x9923329E },
|
||||||
|
{ 0, 465, 0x865B8692, 0x5B9BC5C2 },
|
||||||
|
{ 0, 558, 0x87AA9AFF, 0x79042287 },
|
||||||
|
{ 0, 651, 0x88FCF317, 0xF22241E2 },
|
||||||
|
{ 0, 744, 0x8A5296FF, 0xE33CC930 },
|
||||||
|
{ 0, 837, 0x8BAB8EEF, 0xB6409C1A },
|
||||||
|
{ 0, 930, 0x8D07E334, 0x55637EB3 },
|
||||||
|
{ 0, 1023, 0x8E679C2F, 0x5E44FF8F },
|
||||||
|
{ 0, 1116, 0x8FCAC257, 0x558EE4E6 },
|
||||||
|
{ 0, 1209, 0x91315E37, 0xDB165AA9 },
|
||||||
|
{ 0, 1302, 0x929B7871, 0xDE7F22B9 },
|
||||||
|
{ 0, 1395, 0x940919BB, 0xD4620B6D },
|
||||||
|
{ 0, 1488, 0x957A4AE1, 0xEBF7F3D4 },
|
||||||
|
{ 0, 1581, 0x96EF14C6, 0x454AA840 },
|
||||||
|
{ 0, 1674, 0x98678061, 0x27ECE4F5 },
|
||||||
|
{ 0, 1767, 0x99E396C1, 0x3A3ACFF2 }
|
||||||
|
};
|
||||||
|
|
||||||
|
static flt_arith r_10pow[] = { /* representation of 10 ** -i */
|
||||||
|
{ 0, 0, 0x80000000, 0 },
|
||||||
|
{ 0, -4, 0xCCCCCCCC, 0xCCCCCCCD },
|
||||||
|
{ 0, -7, 0xA3D70A3D, 0x70A3D70A },
|
||||||
|
{ 0, -10, 0x83126E97, 0x8D4FDF3B },
|
||||||
|
{ 0, -14, 0xD1B71758, 0xE219652C },
|
||||||
|
{ 0, -17, 0xA7C5AC47, 0x1B478423 },
|
||||||
|
{ 0, -20, 0x8637BD05, 0xAF6C69B6 },
|
||||||
|
{ 0, -24, 0xD6BF94D5, 0xE57A42BC },
|
||||||
|
{ 0, -27, 0xABCC7711, 0x8461CEFD },
|
||||||
|
{ 0, -30, 0x89705F41, 0x36B4A597 },
|
||||||
|
{ 0, -34, 0xDBE6FECE, 0xBDEDD5BF },
|
||||||
|
{ 0, -37, 0xAFEBFF0B, 0xCB24AAFF },
|
||||||
|
{ 0, -40, 0x8CBCCC09, 0x6F5088CC },
|
||||||
|
{ 0, -44, 0xE12E1342, 0x4BB40E13 },
|
||||||
|
{ 0, -47, 0xB424DC35, 0x095CD80F },
|
||||||
|
{ 0, -50, 0x901D7CF7, 0x3AB0ACD9 },
|
||||||
|
{ 0, -54, 0xE69594BE, 0xC44DE15B },
|
||||||
|
{ 0, -57, 0xB877AA32, 0x36A4B449 },
|
||||||
|
{ 0, -60, 0x9392EE8E, 0x921D5D07 },
|
||||||
|
{ 0, -64, 0xEC1E4A7D, 0xB69561A5 },
|
||||||
|
{ 0, -67, 0xBCE50864, 0x92111AEB },
|
||||||
|
{ 0, -70, 0x971DA050, 0x74DA7BEF },
|
||||||
|
{ 0, -74, 0xF1C90080, 0xBAF72CB1 },
|
||||||
|
{ 0, -77, 0xC16D9A00, 0x95928A27 },
|
||||||
|
{ 0, -80, 0x9ABE14CD, 0x44753B53 },
|
||||||
|
{ 0, -84, 0xF79687AE, 0xD3EEC551 },
|
||||||
|
{ 0, -87, 0xC6120625, 0x76589DDB },
|
||||||
|
{ 0, -90, 0x9E74D1B7, 0x91E07E48 }
|
||||||
|
};
|
||||||
|
|
||||||
|
static flt_arith r_big_10pow[] = { /* representation of 10 ** -(28*i) */
|
||||||
|
{ 0, 0, 0x80000000, 0 },
|
||||||
|
{ 0, -94, 0xFD87B5F2, 0x8300CA0E },
|
||||||
|
{ 0, -187, 0xFB158592, 0xBE068D2F },
|
||||||
|
{ 0, -280, 0xF8A95FCF, 0x88747D94 },
|
||||||
|
{ 0, -373, 0xF64335BC, 0xF065D37D },
|
||||||
|
{ 0, -466, 0xF3E2F893, 0xDEC3F126 },
|
||||||
|
{ 0, -559, 0xF18899B1, 0xBC3F8CA2 },
|
||||||
|
{ 0, -652, 0xEF340A98, 0x172AACE5 },
|
||||||
|
{ 0, -745, 0xECE53CEC, 0x4A314EBE },
|
||||||
|
{ 0, -838, 0xEA9C2277, 0x23EE8BCB },
|
||||||
|
{ 0, -931, 0xE858AD24, 0x8F5C22CA },
|
||||||
|
{ 0, -1024, 0xE61ACF03, 0x3D1A45DF },
|
||||||
|
{ 0, -1117, 0xE3E27A44, 0x4D8D98B8 },
|
||||||
|
{ 0, -1210, 0xE1AFA13A, 0xFBD14D6E },
|
||||||
|
{ 0, -1303, 0xDF82365C, 0x497B5454 },
|
||||||
|
{ 0, -1396, 0xDD5A2C3E, 0xAB3097CC },
|
||||||
|
{ 0, -1489, 0xDB377599, 0xB6074245 },
|
||||||
|
{ 0, -1582, 0xD91A0545, 0xCDB51186 },
|
||||||
|
{ 0, -1675, 0xD701CE3B, 0xD387BF48 },
|
||||||
|
{ 0, -1768, 0xD4EEC394, 0xD6258BF8 }
|
||||||
|
};
|
||||||
|
|
||||||
|
static
|
||||||
|
add_exponent(e, exp)
|
||||||
|
register flt_arith *e;
|
||||||
|
{
|
||||||
|
int neg = exp < 0;
|
||||||
|
int divsz, modsz;
|
||||||
|
flt_arith x;
|
||||||
|
|
||||||
|
if (neg) exp = -exp;
|
||||||
|
divsz = exp / (sizeof(s10pow)/sizeof(s10pow[0]));
|
||||||
|
modsz = exp % (sizeof(s10pow)/sizeof(s10pow[0]));
|
||||||
|
if (neg) {
|
||||||
|
flt_mul(e, &r_10pow[modsz], &x);
|
||||||
|
while (divsz >= 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);
|
||||||
|
}
|
||||||
24
modules/src/flt_arith/misc.h
Normal file
24
modules/src/flt_arith/misc.h
Normal file
|
|
@ -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 <flt_arith.h>
|
||||||
|
|
||||||
|
/* 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
|
||||||
21
modules/src/flt_arith/ucmp.c
Normal file
21
modules/src/flt_arith/ucmp.c
Normal file
|
|
@ -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;
|
||||||
|
}
|
||||||
Loading…
Reference in a new issue