util/int/do_fpar.c

   1 /*
   2  * Sources of the "FLOATING POINT ARITHMETIC" group instructions
   3  */
   4
   5 /* $Id: do_fpar.c,v 2.6 1994/06/24 10:46:21 ceriel Exp $ */
   6
   7 #include        <em_abs.h>
   8 #include        "nofloat.h"
   9 #include        "global.h"
  10 #include        "log.h"
  11 #include        "mem.h"
  12 #include        "trap.h"
  13 #include        "text.h"
  14 #include        "fra.h"
  15 #include        "warn.h"
  16
  17 #ifndef NOFLOAT
  18
  19 extern double fpop();
  20
  21 #ifdef __STDC__
  22 #include <float.h>
  23 #define MAXDOUBLE DBL_MAX
  24 #else /* not __STDC__ */
  25 #if defined(vax) || defined(pdp) || defined(__vax) || defined(__pdp)
  26 #define MAXDOUBLE       1.701411834604692293e+38
  27 #else
  28 #define MAXDOUBLE       1.7976931348623157e+308
  29 #endif
  30 #endif /* not __STDC__ */
  31 #define SMALL           (1.0/MAXDOUBLE)
  32
  33 PRIVATE double adf(), sbf(), mlf(), dvf();
  34 PRIVATE double ttttp();
  35 PRIVATE double floor(), fabs();
  36 PRIVATE fef(), fif();
  37
  38 #endif  /* NOFLOAT */
  39
  40 DoADF(l)
  41         register size l;
  42 {
  43         /* ADF w: Floating add (*) */
  44 #ifndef NOFLOAT
  45         double t = fpop(arg_wf(l));
  46
  47         LOG(("@F6 DoADF(%ld)", l));
  48         spoilFRA();
  49         fpush(adf(fpop(l), t), l);
  50 #else   /* NOFLOAT */
  51         nofloat();
  52 #endif  /* NOFLOAT */
  53 }
  54
  55 DoSBF(l)
  56         register size l;
  57 {
  58         /* SBF w: Floating subtract (*) */
  59 #ifndef NOFLOAT
  60         double t = fpop(arg_wf(l));
  61
  62         LOG(("@F6 DoSBF(%ld)", l));
  63         spoilFRA();
  64         fpush(sbf(fpop(l), t), l);
  65 #else   /* NOFLOAT */
  66         nofloat();
  67 #endif  /* NOFLOAT */
  68 }
  69
  70 DoMLF(l)
  71         register size l;
  72 {
  73         /* MLF w: Floating multiply (*) */
  74 #ifndef NOFLOAT
  75         double t = fpop(arg_wf(l));
  76
  77         LOG(("@F6 DoMLF(%ld)", l));
  78         spoilFRA();
  79         fpush(mlf(fpop(l), t), l);
  80 #else   /* NOFLOAT */
  81         nofloat();
  82 #endif  /* NOFLOAT */
  83 }
  84
  85 DoDVF(l)
  86         register size l;
  87 {
  88         /* DVF w: Floating divide (*) */
  89 #ifndef NOFLOAT
  90         double t = fpop(arg_wf(l));
  91
  92         LOG(("@F6 DoDVF(%ld)", l));
  93         spoilFRA();
  94         fpush(dvf(fpop(l), t), l);
  95 #else   /* NOFLOAT */
  96         nofloat();
  97 #endif  /* NOFLOAT */
  98 }
  99
 100 DoNGF(l)
 101         register size l;
 102 {
 103         /* NGF w: Floating negate (*) */
 104 #ifndef NOFLOAT
 105         double t = fpop(arg_wf(l));
 106
 107         LOG(("@F6 DoNGF(%ld)", l));
 108         spoilFRA();
 109         fpush(-t, l);
 110 #else   /* NOFLOAT */
 111         nofloat();
 112 #endif  /* NOFLOAT */
 113 }
 114
 115 DoFIF(l)
 116         register size l;
 117 {
 118         /* FIF w: Floating multiply and split integer and fraction part (*) */
 119 #ifndef NOFLOAT
 120         double t = fpop(arg_wf(l));
 121
 122         LOG(("@F6 DoFIF(%ld)", l));
 123         spoilFRA();
 124         fif(fpop(l), t, l);
 125 #else   /* NOFLOAT */
 126         nofloat();
 127 #endif  /* NOFLOAT */
 128 }
 129
 130 DoFEF(l)
 131         register size l;
 132 {
 133         /* FEF w: Split floating number in exponent and fraction part (*) */
 134 #ifndef NOFLOAT
 135         LOG(("@F6 DoFEF(%ld)", l));
 136         spoilFRA();
 137         fef(fpop(arg_wf(l)), l);
 138 #else   /* NOFLOAT */
 139         nofloat();
 140 #endif  /* NOFLOAT */
 141 }
 142
 143 #ifndef NOFLOAT
 144
 145 /* Service routines */
 146
 147 PRIVATE double adf(f1, f2)              /* returns f1 + f2 */
 148         double f1, f2;
 149 {
 150         if (must_test && !(IgnMask&BIT(EFOVFL))) {
 151                 if (f1 > 0.0 && f2 > 0.0) {
 152                         if (MAXDOUBLE - f1 < f2) {
 153                                 trap(EFOVFL);
 154                                 return (0.0);
 155                         }
 156                 }
 157                 else if (f1 < 0.0 && f2 < 0.0) {
 158                         if (-(MAXDOUBLE + f1) > f2) {
 159                                 trap(EFOVFL);
 160                                 return (0.0);
 161                         }
 162                 }
 163         }
 164         return (f1 + f2);
 165 }
 166
 167 PRIVATE double sbf(f1, f2)              /* returns f1 - f2 */
 168         double f1, f2;
 169 {
 170         if (must_test && !(IgnMask&BIT(EFOVFL))) {
 171                 if (f2 < 0.0 && f1 > 0.0) {
 172                         if (MAXDOUBLE - f1 < -f2) {
 173                                 trap(EFOVFL);
 174                                 return (0.0);
 175                         }
 176                 }
 177                 else if (f2 > 0.0 && f1 < 0.0) {
 178                         if (f2 - MAXDOUBLE > f1) {
 179                                 trap(EFOVFL);
 180                                 return (0.0);
 181                         }
 182                 }
 183         }
 184         return (f1 - f2);
 185 }
 186
 187 PRIVATE double mlf(f1, f2)              /* returns f1 * f2 */
 188         double f1, f2;
 189 {
 190         double ff1 = fabs(f1), ff2 = fabs(f2);
 191
 192         if (f1 == 0.0 || f2 == 0.0)
 193                 return (0.0);
 194
 195         if ((ff1 >= 1.0 && ff2 <= 1.0) || (ff2 >= 1.0 && ff1 <= 1.0))
 196                 return (f1 * f2);
 197
 198         if (must_test && !(IgnMask&BIT(EFUNFL))) {
 199                 if (ff1 < 1.0 && ff2 < 1.0) {
 200                         if (SMALL / ff1 > ff2) {
 201                                 trap(EFUNFL);
 202                                 return (0.0);
 203                         }
 204                         return (f1 * f2);
 205                 }
 206         }
 207
 208         if (must_test && !(IgnMask&BIT(EFOVFL))) {
 209                 if (MAXDOUBLE / ff1 < ff2) {
 210                         trap(EFOVFL);
 211                         return (0.0);
 212                 }
 213         }
 214         return (f1 * f2);
 215 }
 216
 217 PRIVATE double dvf(f1, f2)              /* returns f1 / f2 */
 218         double f1, f2;
 219 {
 220         double ff1 = fabs(f1), ff2 = fabs(f2);
 221
 222         if (f2 == 0.0) {
 223                 if (!(IgnMask&BIT(EFDIVZ))) {
 224                         trap(EFDIVZ);
 225                 }
 226                 else    return (0.0);
 227         }
 228
 229         if (f1 == 0.0)
 230                 return (0.0);
 231
 232         if ((ff2 >= 1.0 && ff1 >= 1.0) || (ff1 <= 1.0 && ff2 <= 1.0))
 233                 return (f1 / f2);
 234
 235         if (must_test && !(IgnMask&BIT(EFUNFL))) {
 236                 if (ff2 > 1.0 && ff1 < 1.0) {
 237                         if (SMALL / ff2 > ff1) {
 238                                 trap(EFUNFL);
 239                                 return (0.0);
 240                         }
 241                         return (f1 / f2);
 242                 }
 243         }
 244
 245         if (must_test && !(IgnMask&BIT(EFOVFL))) {
 246                 if (MAXDOUBLE * ff2  < ff1) {
 247                         trap(EFOVFL);
 248                         return (0.0);
 249                 }
 250         }
 251         return (f1 / f2);
 252 }
 253
 254 PRIVATE fif(f1, f2, n)
 255         double f1, f2;
 256         size n;
 257 {
 258         double f = mlf(f1, f2);
 259         double fl = floor(fabs(f));
 260
 261         fpush(fabs(f) - fl, n);         /* push fraction */
 262         fpush((f < 0.0) ? -fl : fl, n); /* push integer-part */
 263 }
 264
 265 PRIVATE fef(f, n)
 266         double f;
 267         size n;
 268 {
 269         register long exponent, sign = (long) (f < 0.0);
 270
 271         if (f == 0.0) {
 272                 fpush(f, n);
 273                 wpush(0L);
 274                 return;
 275         }
 276
 277         for (f = fabs(f), exponent = 0; f >= 1.0; exponent++)
 278                 f /= 2.0;
 279
 280         for (; f < 0.5; exponent--)
 281                 f *= 2.0;
 282
 283         fpush((sign) ? -f : f, n);      /* push mantissa */
 284         wpush(exponent);                /* push exponent */
 285 }
 286
 287 /* floating point service routines, to avoid having to use -lm */
 288
 289 PRIVATE double fabs(f)
 290         double f;
 291 {
 292         return (f < 0.0 ? -f : f);
 293 }
 294
 295 PRIVATE double floor(f)
 296         double f;
 297 {
 298         double res, d;
 299         register int sign = 1;
 300
 301         /* eliminate the sign */
 302         if (f < 0) {
 303                 sign = -1, f = -f;
 304         }
 305
 306         /* get the largest power of 2 <= f */
 307         d = 1.0;
 308         while (f - d >= d) {
 309                 d *= 2.0;
 310         }
 311
 312         /* reconstruct f by deminishing powers of 2 */
 313         res = 0.0;
 314         while (d >= 1.0) {
 315                 if (res + d <= f)
 316                         res += d;
 317                 d /= 2.0;
 318         }
 319
 320         /* undo the sign elimination */
 321         if (sign == -1) {
 322                 res = -res, f = -f;
 323                 if (res > f)
 324                         res -= 1.0;
 325         }
 326
 327         return res;
 328 }
 329
 330 PRIVATE double ttttp(f, n)              /* times ten to the power */
 331         double f;
 332 {
 333         while (n > 0) {
 334                 f = mlf(f, 10.0);
 335                 n--;
 336         }
 337         while (n < 0) {
 338                 f = dvf(f, 10.0);
 339                 n++;
 340         }
 341         return f;
 342 }
 343
 344 /*      Str2double is used to initialize the global data area with floats;
 345         we do not use, e.g., sscanf(), to be able to check the grammar of
 346         the string and to give warnings.
 347 */
 348
 349 double str2double(str)
 350         char *str;
 351 {
 352         register char b;
 353         register int sign = 1;          /* either +1 or -1 */
 354         register int frac = 0;          /* how far in fraction part ? */
 355         register int ex;                /* to store exponent */
 356         double mantissa = 0.0;          /* to store mantissa */
 357         double d;                       /* double to be returned */
 358
 359         b = *str++;
 360         if (b == '-') {
 361                 sign = -1;
 362                 b = *str++;
 363         }
 364         else if (b == '+') {
 365                 sign = 1;
 366                 b = *str++;
 367         }
 368
 369         if ('0' <= b && b <= '9') {
 370                 mantissa = (double) (b-'0');
 371         }
 372         else if (b == '.') {
 373                 /* part before dot cannot be empty */
 374                 warning(WBADFLOAT);
 375                 frac = 1;
 376         }
 377         else {
 378                 goto BadFloat;
 379         }
 380
 381         LOG((" q9 str2double : (before while) mantissa = %20.20g", mantissa));
 382
 383         while ((b = *str++) != 'e' && b != 'E' && b != '\0') {
 384                 if (b == '.') {
 385                         if (frac == 0) {
 386                                 frac++;
 387                         }
 388                         else {  /* there already was a '.' in input */
 389                                 goto BadFloat;
 390                         }
 391                 }
 392                 else if ('0' <= b && b <= '9') {
 393                         double bval = b - '0';
 394
 395                         if (frac) {
 396                                 mantissa =
 397                                         adf(mantissa, ttttp(bval, -frac));
 398                                 frac++;
 399                         }
 400                         else {
 401                                 mantissa =
 402                                         adf(mlf(mantissa, 10.0), bval);
 403                         }
 404                 }
 405                 else {
 406                         goto BadFloat;
 407                 }
 408                 LOG((" q9 str2double : (inside while) mantissa = %20.20g",
 409                                                                 mantissa));
 410         }
 411         LOG((" q9 str2double : mantissa = %10.10g", mantissa));
 412         mantissa = sign * mantissa;
 413         if (b == '\0')
 414                 return (mantissa);
 415         /* else we have b == 'e' or b== 'E' */
 416
 417         /* Optional sign for exponent */
 418         b = *str++;
 419         if (b == '-') {
 420                 sign = -1;
 421                 b = *str++;
 422         }
 423         else if (b == '+') {
 424                 sign = 1;
 425                 b = *str++;
 426         }
 427         else {
 428                 sign = 1;
 429         }
 430
 431         ex = 0;
 432         do {
 433                 if ('0' <= b && b <= '9') {
 434                         ex = 10*ex + (b-'0');
 435                 }
 436                 else {
 437                         goto BadFloat;
 438                 }
 439         } while ((b = *str++) != '\0');
 440         LOG((" q9 str2double : exponent = %d", ex));
 441
 442         /* Construct total value of float */
 443         ex = sign * ex;
 444         d = ttttp(mantissa, ex);
 445         return (d);
 446
 447 BadFloat:
 448         fatal("Float garbled in loadfile");
 449         return (0.0);
 450 }
 451
 452 #else   /* NOFLOAT */
 453
 454 nofloat() {
 455         fatal("attempt to execute a floating point instruction on an EM machine without FP");
 456 }
 457
 458 #endif  /* NOFLOAT */
 459