123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459 |
- /*
- * Sources of the "FLOATING POINT ARITHMETIC" group instructions
- */
- /* $Id$ */
- #include <em_abs.h>
- #include "nofloat.h"
- #include "global.h"
- #include "log.h"
- #include "mem.h"
- #include "trap.h"
- #include "text.h"
- #include "fra.h"
- #include "warn.h"
- #ifndef NOFLOAT
- extern double fpop();
- #ifdef __STDC__
- #include <float.h>
- #define MAXDOUBLE DBL_MAX
- #else /* not __STDC__ */
- #if defined(vax) || defined(pdp) || defined(__vax) || defined(__pdp)
- #define MAXDOUBLE 1.701411834604692293e+38
- #else
- #define MAXDOUBLE 1.7976931348623157e+308
- #endif
- #endif /* not __STDC__ */
- #define SMALL (1.0/MAXDOUBLE)
- PRIVATE double adf(), sbf(), mlf(), dvf();
- PRIVATE double ttttp();
- PRIVATE double floor(), fabs();
- PRIVATE fef(), fif();
- #endif /* NOFLOAT */
- DoADF(l)
- register size l;
- {
- /* ADF w: Floating add (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoADF(%ld)", l));
- spoilFRA();
- fpush(adf(fpop(l), t), l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoSBF(l)
- register size l;
- {
- /* SBF w: Floating subtract (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoSBF(%ld)", l));
- spoilFRA();
- fpush(sbf(fpop(l), t), l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoMLF(l)
- register size l;
- {
- /* MLF w: Floating multiply (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoMLF(%ld)", l));
- spoilFRA();
- fpush(mlf(fpop(l), t), l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoDVF(l)
- register size l;
- {
- /* DVF w: Floating divide (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoDVF(%ld)", l));
- spoilFRA();
- fpush(dvf(fpop(l), t), l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoNGF(l)
- register size l;
- {
- /* NGF w: Floating negate (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoNGF(%ld)", l));
- spoilFRA();
- fpush(-t, l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoFIF(l)
- register size l;
- {
- /* FIF w: Floating multiply and split integer and fraction part (*) */
- #ifndef NOFLOAT
- double t = fpop(arg_wf(l));
- LOG(("@F6 DoFIF(%ld)", l));
- spoilFRA();
- fif(fpop(l), t, l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- DoFEF(l)
- register size l;
- {
- /* FEF w: Split floating number in exponent and fraction part (*) */
- #ifndef NOFLOAT
- LOG(("@F6 DoFEF(%ld)", l));
- spoilFRA();
- fef(fpop(arg_wf(l)), l);
- #else /* NOFLOAT */
- nofloat();
- #endif /* NOFLOAT */
- }
- #ifndef NOFLOAT
- /* Service routines */
- PRIVATE double adf(f1, f2) /* returns f1 + f2 */
- double f1, f2;
- {
- if (must_test && !(IgnMask&BIT(EFOVFL))) {
- if (f1 > 0.0 && f2 > 0.0) {
- if (MAXDOUBLE - f1 < f2) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- else if (f1 < 0.0 && f2 < 0.0) {
- if (-(MAXDOUBLE + f1) > f2) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- }
- return (f1 + f2);
- }
- PRIVATE double sbf(f1, f2) /* returns f1 - f2 */
- double f1, f2;
- {
- if (must_test && !(IgnMask&BIT(EFOVFL))) {
- if (f2 < 0.0 && f1 > 0.0) {
- if (MAXDOUBLE - f1 < -f2) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- else if (f2 > 0.0 && f1 < 0.0) {
- if (f2 - MAXDOUBLE > f1) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- }
- return (f1 - f2);
- }
- PRIVATE double mlf(f1, f2) /* returns f1 * f2 */
- double f1, f2;
- {
- double ff1 = fabs(f1), ff2 = fabs(f2);
- if (f1 == 0.0 || f2 == 0.0)
- return (0.0);
- if ((ff1 >= 1.0 && ff2 <= 1.0) || (ff2 >= 1.0 && ff1 <= 1.0))
- return (f1 * f2);
- if (must_test && !(IgnMask&BIT(EFUNFL))) {
- if (ff1 < 1.0 && ff2 < 1.0) {
- if (SMALL / ff1 > ff2) {
- trap(EFUNFL);
- return (0.0);
- }
- return (f1 * f2);
- }
- }
- if (must_test && !(IgnMask&BIT(EFOVFL))) {
- if (MAXDOUBLE / ff1 < ff2) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- return (f1 * f2);
- }
- PRIVATE double dvf(f1, f2) /* returns f1 / f2 */
- double f1, f2;
- {
- double ff1 = fabs(f1), ff2 = fabs(f2);
- if (f2 == 0.0) {
- if (!(IgnMask&BIT(EFDIVZ))) {
- trap(EFDIVZ);
- }
- else return (0.0);
- }
- if (f1 == 0.0)
- return (0.0);
- if ((ff2 >= 1.0 && ff1 >= 1.0) || (ff1 <= 1.0 && ff2 <= 1.0))
- return (f1 / f2);
- if (must_test && !(IgnMask&BIT(EFUNFL))) {
- if (ff2 > 1.0 && ff1 < 1.0) {
- if (SMALL / ff2 > ff1) {
- trap(EFUNFL);
- return (0.0);
- }
- return (f1 / f2);
- }
- }
- if (must_test && !(IgnMask&BIT(EFOVFL))) {
- if (MAXDOUBLE * ff2 < ff1) {
- trap(EFOVFL);
- return (0.0);
- }
- }
- return (f1 / f2);
- }
- PRIVATE fif(f1, f2, n)
- double f1, f2;
- size n;
- {
- double f = mlf(f1, f2);
- double fl = floor(fabs(f));
- fpush(fabs(f) - fl, n); /* push fraction */
- fpush((f < 0.0) ? -fl : fl, n); /* push integer-part */
- }
- PRIVATE fef(f, n)
- double f;
- size n;
- {
- register long exponent, sign = (long) (f < 0.0);
- if (f == 0.0) {
- fpush(f, n);
- wpush(0L);
- return;
- }
- for (f = fabs(f), exponent = 0; f >= 1.0; exponent++)
- f /= 2.0;
- for (; f < 0.5; exponent--)
- f *= 2.0;
- fpush((sign) ? -f : f, n); /* push mantissa */
- wpush(exponent); /* push exponent */
- }
- /* floating point service routines, to avoid having to use -lm */
- PRIVATE double fabs(f)
- double f;
- {
- return (f < 0.0 ? -f : f);
- }
- PRIVATE double floor(f)
- double f;
- {
- double res, d;
- register int sign = 1;
-
- /* eliminate the sign */
- if (f < 0) {
- sign = -1, f = -f;
- }
-
- /* get the largest power of 2 <= f */
- d = 1.0;
- while (f - d >= d) {
- d *= 2.0;
- }
-
- /* reconstruct f by deminishing powers of 2 */
- res = 0.0;
- while (d >= 1.0) {
- if (res + d <= f)
- res += d;
- d /= 2.0;
- }
-
- /* undo the sign elimination */
- if (sign == -1) {
- res = -res, f = -f;
- if (res > f)
- res -= 1.0;
- }
-
- return res;
- }
- PRIVATE double ttttp(f, n) /* times ten to the power */
- double f;
- {
- while (n > 0) {
- f = mlf(f, 10.0);
- n--;
- }
- while (n < 0) {
- f = dvf(f, 10.0);
- n++;
- }
- return f;
- }
- /* Str2double is used to initialize the global data area with floats;
- we do not use, e.g., sscanf(), to be able to check the grammar of
- the string and to give warnings.
- */
- double str2double(str)
- char *str;
- {
- register char b;
- register int sign = 1; /* either +1 or -1 */
- register int frac = 0; /* how far in fraction part ? */
- register int ex; /* to store exponent */
- double mantissa = 0.0; /* to store mantissa */
- double d; /* double to be returned */
-
- b = *str++;
- if (b == '-') {
- sign = -1;
- b = *str++;
- }
- else if (b == '+') {
- sign = 1;
- b = *str++;
- }
-
- if ('0' <= b && b <= '9') {
- mantissa = (double) (b-'0');
- }
- else if (b == '.') {
- /* part before dot cannot be empty */
- warning(WBADFLOAT);
- frac = 1;
- }
- else {
- goto BadFloat;
- }
-
- LOG((" q9 str2double : (before while) mantissa = %20.20g", mantissa));
-
- while ((b = *str++) != 'e' && b != 'E' && b != '\0') {
- if (b == '.') {
- if (frac == 0) {
- frac++;
- }
- else { /* there already was a '.' in input */
- goto BadFloat;
- }
- }
- else if ('0' <= b && b <= '9') {
- double bval = b - '0';
-
- if (frac) {
- mantissa =
- adf(mantissa, ttttp(bval, -frac));
- frac++;
- }
- else {
- mantissa =
- adf(mlf(mantissa, 10.0), bval);
- }
- }
- else {
- goto BadFloat;
- }
- LOG((" q9 str2double : (inside while) mantissa = %20.20g",
- mantissa));
- }
- LOG((" q9 str2double : mantissa = %10.10g", mantissa));
- mantissa = sign * mantissa;
- if (b == '\0')
- return (mantissa);
- /* else we have b == 'e' or b== 'E' */
-
- /* Optional sign for exponent */
- b = *str++;
- if (b == '-') {
- sign = -1;
- b = *str++;
- }
- else if (b == '+') {
- sign = 1;
- b = *str++;
- }
- else {
- sign = 1;
- }
-
- ex = 0;
- do {
- if ('0' <= b && b <= '9') {
- ex = 10*ex + (b-'0');
- }
- else {
- goto BadFloat;
- }
- } while ((b = *str++) != '\0');
- LOG((" q9 str2double : exponent = %d", ex));
-
- /* Construct total value of float */
- ex = sign * ex;
- d = ttttp(mantissa, ex);
- return (d);
- BadFloat:
- fatal("Float garbled in loadfile");
- return (0.0);
- }
- #else /* NOFLOAT */
- nofloat() {
- fatal("attempt to execute a floating point instruction on an EM machine without FP");
- }
- #endif /* NOFLOAT */
|