gcc/real.c

   1 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
   2 and support for XFmode IEEE extended real floating point arithmetic.
   3 Contributed by Stephen L. Moshier (moshier@world.std.com).
   4
   5    Copyright (C) 1993 Free Software Foundation, Inc.
   6
   7 This file is part of GNU CC.
   8
   9 GNU CC is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU General Public License as published by
  11 the Free Software Foundation; either version 2, or (at your option)
  12 any later version.
  13
  14 GNU CC is distributed in the hope that it will be useful,
  15 but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 GNU General Public License for more details.
  18
  19 You should have received a copy of the GNU General Public License
  20 along with GNU CC; see the file COPYING.  If not, write to
  21 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  22
  23 #include <stdio.h>
  24 #include <errno.h>
  25 #include "config.h"
  26 #include "tree.h"
  27
  28 #ifndef errno
  29 extern int errno;
  30 #endif
  31
  32 /* To enable support of XFmode extended real floating point, define
  33 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  34
  35 To support cross compilation between IEEE and VAX floating
  36 point formats, define REAL_ARITHMETIC in the tm.h file.
  37
  38 In either case the machine files (tm.h) must not contain any code
  39 that tries to use host floating point arithmetic to convert
  40 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  41 etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  42 be intelligible to the host computer's native arithmetic.
  43
  44 The emulator defaults to the host's floating point format so that
  45 its decimal conversion functions can be used if desired (see
  46 real.h).
  47
  48 The first part of this file interfaces gcc to ieee.c, which is a
  49 floating point arithmetic suite that was not written with gcc in
  50 mind.  The interface is followed by ieee.c itself and related
  51 items. Avoid changing ieee.c unless you have suitable test
  52 programs available.  A special version of the PARANOIA floating
  53 point arithmetic tester, modified for this purpose, can be found
  54 on usc.edu : /pub/C-numanal/ieeetest.zoo.  Some tutorial
  55 information on ieee.c is given in my book: S. L. Moshier,
  56 _Methods and Programs for Mathematical Functions_, Prentice-Hall
  57 or Simon & Schuster Int'l, 1989.  A library of XFmode elementary
  58 transcendental functions can be obtained by ftp from
  59 research.att.com: netlib/cephes/ldouble.shar.Z  */
  60
  61 /* Type of computer arithmetic.
  62  * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
  63  * The following modification converts gcc macros into the ones
  64  * used by ieee.c.
  65  *
  66  * Note: long double formats differ between IBMPC and MIEEE
  67  * by more than just endian-ness.
  68  */
  69
  70 /* REAL_ARITHMETIC defined means that macros in real.h are
  71    defined to call emulator functions.  */
  72 #ifdef REAL_ARITHMETIC
  73
  74 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  75 /* PDP-11, Pro350, VAX: */
  76 #define DEC 1
  77 #else /* it's not VAX */
  78 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  79 #if WORDS_BIG_ENDIAN
  80 /* Motorola IEEE, high order words come first (Sun workstation): */
  81 #define MIEEE 1
  82 #else /* not big-endian */
  83 /* Intel IEEE, low order words come first:
  84  */
  85 #define IBMPC 1
  86 #endif /*  big-endian */
  87 #else /* it's not IEEE either */
  88 /* UNKnown arithmetic.  We don't support this and can't go on. */
  89 unknown arithmetic type
  90 #define UNK 1
  91 #endif /* not IEEE */
  92 #endif /* not VAX */
  93
  94 #else
  95 /* REAL_ARITHMETIC not defined means that the *host's* data
  96    structure will be used.  It may differ by endian-ness from the
  97    target machine's structure and will get its ends swapped
  98    accordingly (but not here).  Probably only the decimal <-> binary
  99    functions in this file will actually be used in this case.  */
 100 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 101 #define DEC 1
 102 #else /* it's not VAX */
 103 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 104 #ifdef HOST_WORDS_BIG_ENDIAN
 105 #define MIEEE 1
 106 #else /* not big-endian */
 107 #define IBMPC 1
 108 #endif /*  big-endian */
 109 #else /* it's not IEEE either */
 110 unknown arithmetic type
 111 #define UNK 1
 112 #endif /* not IEEE */
 113 #endif /* not VAX */
 114
 115 #endif /* REAL_ARITHMETIC not defined */
 116
 117 /* Define for support of infinity.  */
 118 #ifndef DEC
 119 #define INFINITY
 120 #endif
 121
 122
 123 /* ehead.h
 124  *
 125  * Include file for extended precision arithmetic programs.
 126  */
 127
 128 /* Number of 16 bit words in external e type format */
 129 #define NE 6
 130
 131 /* Number of 16 bit words in internal format */
 132 #define NI (NE+3)
 133
 134 /* Array offset to exponent */
 135 #define E 1
 136
 137 /* Array offset to high guard word */
 138 #define M 2
 139
 140 /* Number of bits of precision */
 141 #define NBITS ((NI-4)*16)
 142
 143 /* Maximum number of decimal digits in ASCII conversion
 144  * = NBITS*log10(2)
 145  */
 146 #define NDEC (NBITS*8/27)
 147
 148 /* The exponent of 1.0 */
 149 #define EXONE (0x3fff)
 150
 151 /* Find a host integer type that is at least 16 bits wide,
 152    and another type at least twice whatever that size is. */
 153
 154 #if HOST_BITS_PER_CHAR >= 16
 155 #define EMUSHORT char
 156 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
 157 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
 158 #else
 159 #if HOST_BITS_PER_SHORT >= 16
 160 #define EMUSHORT short
 161 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
 162 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
 163 #else
 164 #if HOST_BITS_PER_INT >= 16
 165 #define EMUSHORT int
 166 #define EMUSHORT_SIZE HOST_BITS_PER_INT
 167 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
 168 #else
 169 #if HOST_BITS_PER_LONG >= 16
 170 #define EMUSHORT long
 171 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 172 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 173 #else
 174 /*  You will have to modify this program to have a smaller unit size. */
 175 #define EMU_NON_COMPILE
 176 #endif
 177 #endif
 178 #endif
 179 #endif
 180
 181 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
 182 #define EMULONG short
 183 #else
 184 #if HOST_BITS_PER_INT >= EMULONG_SIZE
 185 #define EMULONG int
 186 #else
 187 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
 188 #define EMULONG long
 189 #else
 190 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
 191 #define EMULONG long long int
 192 #else
 193 /*  You will have to modify this program to have a smaller unit size. */
 194 #define EMU_NON_COMPILE
 195 #endif
 196 #endif
 197 #endif
 198 #endif
 199
 200
 201 /* The host interface doesn't work if no 16-bit size exists. */
 202 #if EMUSHORT_SIZE != 16
 203 #define EMU_NON_COMPILE
 204 #endif
 205
 206 /* OK to continue compilation. */
 207 #ifndef EMU_NON_COMPILE
 208
 209 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
 210    In GET_REAL and PUT_REAL, r and e are pointers.
 211    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
 212    in memory, with no holes.  */
 213
 214 #if LONG_DOUBLE_TYPE_SIZE == 96
 215 #define GET_REAL(r,e) bcopy (r, e, 2*NE)
 216 #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
 217 #else /* no XFmode */
 218
 219 #ifdef REAL_ARITHMETIC
 220 /* Emulator uses target format internally
 221    but host stores it in host endian-ness. */
 222
 223 #if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
 224 #define GET_REAL(r,e) e53toe ((r), (e))
 225 #define PUT_REAL(e,r) etoe53 ((e), (r))
 226
 227 #else /* endian-ness differs */
 228 /* emulator uses target endian-ness internally */
 229 #define GET_REAL(r,e)           \
 230 do { EMUSHORT w[4];             \
 231  w[3] = ((EMUSHORT *) r)[0];    \
 232  w[2] = ((EMUSHORT *) r)[1];    \
 233  w[1] = ((EMUSHORT *) r)[2];    \
 234  w[0] = ((EMUSHORT *) r)[3];    \
 235  e53toe (w, (e)); } while (0)
 236
 237 #define PUT_REAL(e,r)           \
 238 do { EMUSHORT w[4];             \
 239  etoe53 ((e), w);               \
 240  *((EMUSHORT *) r) = w[3];      \
 241  *((EMUSHORT *) r + 1) = w[2];  \
 242  *((EMUSHORT *) r + 2) = w[1];  \
 243  *((EMUSHORT *) r + 3) = w[0]; } while (0)
 244
 245 #endif /* endian-ness differs */
 246
 247 #else /* not REAL_ARITHMETIC */
 248
 249 /* emulator uses host format */
 250 #define GET_REAL(r,e) e53toe ((r), (e))
 251 #define PUT_REAL(e,r) etoe53 ((e), (r))
 252
 253 #endif /* not REAL_ARITHMETIC */
 254 #endif /* no XFmode */
 255
 256 void warning ();
 257 extern int extra_warnings;
 258 int ecmp (), enormlz (), eshift (), eisneg (), eisinf ();
 259 void eadd (), esub (), emul (), ediv ();
 260 void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
 261 void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
 262 void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
 263 void eround (), ereal_to_decimal ();
 264 void esqrt (), elog (), eexp (), etanh (), epow ();
 265 void asctoe (), asctoe24 (), asctoe53 (), asctoe64 ();
 266 void etoasc (), e24toasc (), e53toasc (), e64toasc ();
 267 void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
 268 void mtherr ();
 269 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 270 extern unsigned EMUSHORT elog2[], esqrt2[];
 271
 272 /* Pack output array with 32-bit numbers obtained from
 273    array containing 16-bit numbers, swapping ends if required. */
 274 void
 275 endian (e, x, mode)
 276      unsigned EMUSHORT e[];
 277      long x[];
 278      enum machine_mode mode;
 279 {
 280   unsigned long th, t;
 281
 282 #if WORDS_BIG_ENDIAN
 283   switch (mode)
 284     {
 285
 286     case XFmode:
 287
 288       /* Swap halfwords in the third long. */
 289       th = (unsigned long) e[4] & 0xffff;
 290       t = (unsigned long) e[5] & 0xffff;
 291       t |= th << 16;
 292       x[2] = (long) t;
 293       /* fall into the double case */
 294
 295     case DFmode:
 296
 297       /* swap halfwords in the second word */
 298       th = (unsigned long) e[2] & 0xffff;
 299       t = (unsigned long) e[3] & 0xffff;
 300       t |= th << 16;
 301       x[1] = (long) t;
 302       /* fall into the float case */
 303
 304     case SFmode:
 305
 306       /* swap halfwords in the first word */
 307       th = (unsigned long) e[0] & 0xffff;
 308       t = (unsigned long) e[1] & 0xffff;
 309       t |= th << 16;
 310       x[0] = t;
 311       break;
 312
 313     default:
 314       abort ();
 315     }
 316
 317 #else
 318
 319   /* Pack the output array without swapping. */
 320
 321   switch (mode)
 322     {
 323
 324     case XFmode:
 325
 326       /* Pack the third long.
 327          Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
 328          in it.  */
 329       th = (unsigned long) e[5] & 0xffff;
 330       t = (unsigned long) e[4] & 0xffff;
 331       t |= th << 16;
 332       x[2] = (long) t;
 333       /* fall into the double case */
 334
 335     case DFmode:
 336
 337       /* pack the second long */
 338       th = (unsigned long) e[3] & 0xffff;
 339       t = (unsigned long) e[2] & 0xffff;
 340       t |= th << 16;
 341       x[1] = (long) t;
 342       /* fall into the float case */
 343
 344     case SFmode:
 345
 346       /* pack the first long */
 347       th = (unsigned long) e[1] & 0xffff;
 348       t = (unsigned long) e[0] & 0xffff;
 349       t |= th << 16;
 350       x[0] = t;
 351       break;
 352
 353     default:
 354       abort ();
 355     }
 356
 357 #endif
 358 }
 359
 360
 361 /* This is the implementation of the REAL_ARITHMETIC macro.
 362  */
 363 void
 364 earith (value, icode, r1, r2)
 365      REAL_VALUE_TYPE *value;
 366      int icode;
 367      REAL_VALUE_TYPE *r1;
 368      REAL_VALUE_TYPE *r2;
 369 {
 370   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
 371   enum tree_code code;
 372
 373   GET_REAL (r1, d1);
 374   GET_REAL (r2, d2);
 375   code = (enum tree_code) icode;
 376   switch (code)
 377     {
 378     case PLUS_EXPR:
 379       eadd (d2, d1, v);
 380       break;
 381
 382     case MINUS_EXPR:
 383       esub (d2, d1, v);         /* d1 - d2 */
 384       break;
 385
 386     case MULT_EXPR:
 387       emul (d2, d1, v);
 388       break;
 389
 390     case RDIV_EXPR:
 391 #ifndef REAL_INFINITY
 392       if (ecmp (d2, ezero) == 0)
 393         abort ();
 394 #endif
 395       ediv (d2, d1, v); /* d1/d2 */
 396       break;
 397
 398     case MIN_EXPR:              /* min (d1,d2) */
 399       if (ecmp (d1, d2) < 0)
 400         emov (d1, v);
 401       else
 402         emov (d2, v);
 403       break;
 404
 405     case MAX_EXPR:              /* max (d1,d2) */
 406       if (ecmp (d1, d2) > 0)
 407         emov (d1, v);
 408       else
 409         emov (d2, v);
 410       break;
 411     default:
 412       emov (ezero, v);
 413       break;
 414     }
 415 PUT_REAL (v, value);
 416 }
 417
 418
 419 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
 420  * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x))
 421  */
 422 REAL_VALUE_TYPE
 423 etrunci (x)
 424      REAL_VALUE_TYPE x;
 425 {
 426   unsigned EMUSHORT f[NE], g[NE];
 427   REAL_VALUE_TYPE r;
 428   long l;
 429
 430   GET_REAL (&x, g);
 431   eifrac (g, &l, f);
 432   ltoe (&l, g);
 433   PUT_REAL (g, &r);
 434   return (r);
 435 }
 436
 437
 438 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
 439  * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x))
 440  */
 441 REAL_VALUE_TYPE
 442 etruncui (x)
 443      REAL_VALUE_TYPE x;
 444 {
 445   unsigned EMUSHORT f[NE], g[NE];
 446   REAL_VALUE_TYPE r;
 447   unsigned long l;
 448
 449   GET_REAL (&x, g);
 450   euifrac (g, &l, f);
 451   ultoe (&l, g);
 452   PUT_REAL (g, &r);
 453   return (r);
 454 }
 455
 456
 457 /* This is the REAL_VALUE_ATOF function.
 458  * It converts a decimal string to binary, rounding off
 459  * as indicated by the machine_mode argument.  Then it
 460  * promotes the rounded value to REAL_VALUE_TYPE.
 461  */
 462 REAL_VALUE_TYPE
 463 ereal_atof (s, t)
 464      char *s;
 465      enum machine_mode t;
 466 {
 467   unsigned EMUSHORT tem[NE], e[NE];
 468   REAL_VALUE_TYPE r;
 469
 470   switch (t)
 471     {
 472     case SFmode:
 473       asctoe24 (s, tem);
 474       e24toe (tem, e);
 475       break;
 476     case DFmode:
 477       asctoe53 (s, tem);
 478       e53toe (tem, e);
 479       break;
 480     case XFmode:
 481       asctoe64 (s, tem);
 482       e64toe (tem, e);
 483       break;
 484     default:
 485       asctoe (s, e);
 486     }
 487   PUT_REAL (e, &r);
 488   return (r);
 489 }
 490
 491
 492 /* Expansion of REAL_NEGATE.
 493  */
 494 REAL_VALUE_TYPE
 495 ereal_negate (x)
 496      REAL_VALUE_TYPE x;
 497 {
 498   unsigned EMUSHORT e[NE];
 499   REAL_VALUE_TYPE r;
 500
 501   GET_REAL (&x, e);
 502   eneg (e);
 503   PUT_REAL (e, &r);
 504   return (r);
 505 }
 506
 507
 508 /* Round real to int
 509  * implements REAL_VALUE_FIX (x) (eroundi (x))
 510  * The type of rounding is left unspecified by real.h.
 511  * It is implemented here as round to nearest (add .5 and chop).
 512  */
 513 int
 514 eroundi (x)
 515      REAL_VALUE_TYPE x;
 516 {
 517   unsigned EMUSHORT f[NE], g[NE];
 518   EMULONG l;
 519
 520   GET_REAL (&x, f);
 521   eround (f, g);
 522   eifrac (g, &l, f);
 523   return ((int) l);
 524 }
 525
 526 /* Round real to nearest unsigned int
 527  * implements  REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
 528  * Negative input returns zero.
 529  * The type of rounding is left unspecified by real.h.
 530  * It is implemented here as round to nearest (add .5 and chop).
 531  */
 532 unsigned int
 533 eroundui (x)
 534      REAL_VALUE_TYPE x;
 535 {
 536   unsigned EMUSHORT f[NE], g[NE];
 537   unsigned EMULONG l;
 538
 539   GET_REAL (&x, f);
 540   eround (f, g);
 541   euifrac (g, &l, f);
 542   return ((unsigned int)l);
 543 }
 544
 545
 546 /* REAL_VALUE_FROM_INT macro.
 547  */
 548 void
 549 ereal_from_int (d, i, j)
 550      REAL_VALUE_TYPE *d;
 551      long i, j;
 552 {
 553   unsigned EMUSHORT df[NE], dg[NE];
 554   long low, high;
 555   int sign;
 556
 557   sign = 0;
 558   low = i;
 559   if ((high = j) < 0)
 560     {
 561       sign = 1;
 562       /* complement and add 1 */
 563       high = ~high;
 564       if (low)
 565         low = -low;
 566       else
 567         high += 1;
 568     }
 569   eldexp (eone, HOST_BITS_PER_LONG, df);
 570   ultoe (&high, dg);
 571   emul (dg, df, dg);
 572   ultoe (&low, df);
 573   eadd (df, dg, dg);
 574   if (sign)
 575     eneg (dg);
 576   PUT_REAL (dg, d);
 577 }
 578
 579
 580 /* REAL_VALUE_FROM_UNSIGNED_INT macro.
 581  */
 582 void
 583 ereal_from_uint (d, i, j)
 584      REAL_VALUE_TYPE *d;
 585      unsigned long i, j;
 586 {
 587   unsigned EMUSHORT df[NE], dg[NE];
 588   unsigned long low, high;
 589
 590   low = i;
 591   high = j;
 592   eldexp (eone, HOST_BITS_PER_LONG, df);
 593   ultoe (&high, dg);
 594   emul (dg, df, dg);
 595   ultoe (&low, df);
 596   eadd (df, dg, dg);
 597   PUT_REAL (dg, d);
 598 }
 599
 600
 601 /* REAL_VALUE_TO_INT macro
 602  */
 603 void
 604 ereal_to_int (low, high, rr)
 605      long *low, *high;
 606      REAL_VALUE_TYPE rr;
 607 {
 608   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
 609   int s;
 610
 611   GET_REAL (&rr, d);
 612   /* convert positive value */
 613   s = 0;
 614   if (eisneg (d))
 615     {
 616       eneg (d);
 617       s = 1;
 618     }
 619   eldexp (eone, HOST_BITS_PER_LONG, df);
 620   ediv (df, d, dg);             /* dg = d / 2^32 is the high word */
 621   euifrac (dg, high, dh);
 622   emul (df, dh, dg);            /* fractional part is the low word */
 623   euifrac (dg, low, dh);
 624   if (s)
 625     {
 626       /* complement and add 1 */
 627       *high = ~(*high);
 628       if (*low)
 629         *low = -(*low);
 630       else
 631         *high += 1;
 632     }
 633 }
 634
 635
 636 /* REAL_VALUE_LDEXP macro.
 637  */
 638 REAL_VALUE_TYPE
 639 ereal_ldexp (x, n)
 640      REAL_VALUE_TYPE x;
 641      int n;
 642 {
 643   unsigned EMUSHORT e[NE], y[NE];
 644   REAL_VALUE_TYPE r;
 645
 646   GET_REAL (&x, e);
 647   eldexp (e, n, y);
 648   PUT_REAL (y, &r);
 649   return (r);
 650 }
 651
 652 /* These routines are conditionally compiled because functions
 653  * of the same names may be defined in fold-const.c.  */
 654 #ifdef REAL_ARITHMETIC
 655
 656 /* Check for infinity in a REAL_VALUE_TYPE. */
 657 int
 658 target_isinf (x)
 659      REAL_VALUE_TYPE x;
 660 {
 661   unsigned EMUSHORT e[NE];
 662
 663 #ifdef INFINITY
 664   GET_REAL (&x, e);
 665   return (eisinf (e));
 666 #else
 667   return 0;
 668 #endif
 669 }
 670
 671
 672 /* Check whether an IEEE double precision number is a NaN. */
 673
 674 int
 675 target_isnan (x)
 676      REAL_VALUE_TYPE x;
 677 {
 678   return (0);
 679 }
 680
 681
 682 /* Check for a negative IEEE double precision number.
 683  * this means strictly less than zero, not -0.
 684  */
 685
 686 int
 687 target_negative (x)
 688      REAL_VALUE_TYPE x;
 689 {
 690   unsigned EMUSHORT e[NE];
 691
 692   GET_REAL (&x, e);
 693   if (ecmp (e, ezero) < 0)
 694     return (1);
 695   return (0);
 696 }
 697
 698 /* Expansion of REAL_VALUE_TRUNCATE.
 699  * The result is in floating point, rounded to nearest or even.
 700  */
 701 REAL_VALUE_TYPE
 702 real_value_truncate (mode, arg)
 703      enum machine_mode mode;
 704      REAL_VALUE_TYPE arg;
 705 {
 706   unsigned EMUSHORT e[NE], t[NE];
 707   REAL_VALUE_TYPE r;
 708
 709   GET_REAL (&arg, e);
 710   eclear (t);
 711   switch (mode)
 712     {
 713     case XFmode:
 714       etoe64 (e, t);
 715       e64toe (t, t);
 716       break;
 717
 718     case DFmode:
 719       etoe53 (e, t);
 720       e53toe (t, t);
 721       break;
 722
 723     case SFmode:
 724       etoe24 (e, t);
 725       e24toe (t, t);
 726       break;
 727
 728     case SImode:
 729       r = etrunci (e);
 730       return (r);
 731
 732     default:
 733       abort ();
 734     }
 735   PUT_REAL (t, &r);
 736   return (r);
 737 }
 738
 739 #endif /* REAL_ARITHMETIC defined */
 740
 741 /* Target values are arrays of host longs. A long is guaranteed
 742    to be at least 32 bits wide. */
 743 void
 744 etarldouble (r, l)
 745      REAL_VALUE_TYPE r;
 746      long l[];
 747 {
 748   unsigned EMUSHORT e[NE];
 749
 750   GET_REAL (&r, e);
 751   etoe64 (e, e);
 752   endian (e, l, XFmode);
 753 }
 754
 755 void
 756 etardouble (r, l)
 757      REAL_VALUE_TYPE r;
 758      long l[];
 759 {
 760   unsigned EMUSHORT e[NE];
 761
 762   GET_REAL (&r, e);
 763   etoe53 (e, e);
 764   endian (e, l, DFmode);
 765 }
 766
 767 long
 768 etarsingle (r)
 769      REAL_VALUE_TYPE r;
 770 {
 771   unsigned EMUSHORT e[NE];
 772   unsigned long l;
 773
 774   GET_REAL (&r, e);
 775   etoe24 (e, e);
 776   endian (e, &l, SFmode);
 777   return ((long) l);
 778 }
 779
 780 void
 781 ereal_to_decimal (x, s)
 782      REAL_VALUE_TYPE x;
 783      char *s;
 784 {
 785   unsigned EMUSHORT e[NE];
 786
 787   GET_REAL (&x, e);
 788   etoasc (e, s, 20);
 789 }
 790
 791 int
 792 ereal_cmp (x, y)
 793      REAL_VALUE_TYPE x, y;
 794 {
 795   unsigned EMUSHORT ex[NE], ey[NE];
 796
 797   GET_REAL (&x, ex);
 798   GET_REAL (&y, ey);
 799   return (ecmp (ex, ey));
 800 }
 801
 802 int
 803 ereal_isneg (x)
 804      REAL_VALUE_TYPE x;
 805 {
 806   unsigned EMUSHORT ex[NE];
 807
 808   GET_REAL (&x, ex);
 809   return (eisneg (ex));
 810 }
 811
 812 /* End of REAL_ARITHMETIC interface */
 813
 814 /*                                                      ieee.c
 815  *
 816  *    Extended precision IEEE binary floating point arithmetic routines
 817  *
 818  * Numbers are stored in C language as arrays of 16-bit unsigned
 819  * short integers.  The arguments of the routines are pointers to
 820  * the arrays.
 821  *
 822  *
 823  * External e type data structure, simulates Intel 8087 chip
 824  * temporary real format but possibly with a larger significand:
 825  *
 826  *      NE-1 significand words  (least significant word first,
 827  *                               most significant bit is normally set)
 828  *      exponent                (value = EXONE for 1.0,
 829  *                              top bit is the sign)
 830  *
 831  *
 832  * Internal data structure of a number (a "word" is 16 bits):
 833  *
 834  * ei[0]        sign word       (0 for positive, 0xffff for negative)
 835  * ei[1]        biased exponent (value = EXONE for the number 1.0)
 836  * ei[2]        high guard word (always zero after normalization)
 837  * ei[3]
 838  * to ei[NI-2]  significand     (NI-4 significand words,
 839  *                               most significant word first,
 840  *                               most significant bit is set)
 841  * ei[NI-1]     low guard word  (0x8000 bit is rounding place)
 842  *
 843  *
 844  *
 845  *              Routines for external format numbers
 846  *
 847  *      asctoe (string, e)      ASCII string to extended double e type
 848  *      asctoe64 (string, &d)   ASCII string to long double
 849  *      asctoe53 (string, &d)   ASCII string to double
 850  *      asctoe24 (string, &f)   ASCII string to single
 851  *      asctoeg (string, e, prec) ASCII string to specified precision
 852  *      e24toe (&f, e)          IEEE single precision to e type
 853  *      e53toe (&d, e)          IEEE double precision to e type
 854  *      e64toe (&d, e)          IEEE long double precision to e type
 855  *      eabs (e)                        absolute value
 856  *      eadd (a, b, c)          c = b + a
 857  *      eclear (e)              e = 0
 858  *      ecmp (a, b)             compare a to b, return 1, 0, or -1
 859  *      ediv (a, b, c)          c = b / a
 860  *      efloor (a, b)           truncate to integer, toward -infinity
 861  *      efrexp (a, exp, s)      extract exponent and significand
 862  *      eifrac (e, &l, frac)    e to long integer and e type fraction
 863  *      euifrac (e, &l, frac)   e to unsigned long integer and e type fraction
 864  *      einfin (e)              set e to infinity, leaving its sign alone
 865  *      eldexp (a, n, b)        multiply by 2**n
 866  *      emov (a, b)             b = a
 867  *      emul (a, b, c)          c = b * a
 868  *      eneg (e)                        e = -e
 869  *      eround (a, b)           b = nearest integer value to a
 870  *      esub (a, b, c)          c = b - a
 871  *      e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
 872  *      e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
 873  *      e64toasc (&d, str, n)   long double to ASCII string
 874  *      etoasc (e, str, n)      e to ASCII string, n digits after decimal
 875  *      etoe24 (e, &f)          convert e type to IEEE single precision
 876  *      etoe53 (e, &d)          convert e type to IEEE double precision
 877  *      etoe64 (e, &d)          convert e type to IEEE long double precision
 878  *      ltoe (&l, e)            long (32 bit) integer to e type
 879  *      ultoe (&l, e)           unsigned long (32 bit) integer to e type
 880  *      eisneg (e)              1 if sign bit of e != 0, else 0
 881  *      eisinf (e)              1 if e has maximum exponent
 882  *
 883  *
 884  *              Routines for internal format numbers
 885  *
 886  *      eaddm (ai, bi)          add significands, bi = bi + ai
 887  *      ecleaz (ei)             ei = 0
 888  *      ecleazs (ei)            set ei = 0 but leave its sign alone
 889  *      ecmpm (ai, bi)          compare significands, return 1, 0, or -1
 890  *      edivm (ai, bi)          divide  significands, bi = bi / ai
 891  *      emdnorm (ai,l,s,exp)    normalize and round off
 892  *      emovi (a, ai)           convert external a to internal ai
 893  *      emovo (ai, a)           convert internal ai to external a
 894  *      emovz (ai, bi)          bi = ai, low guard word of bi = 0
 895  *      emulm (ai, bi)          multiply significands, bi = bi * ai
 896  *      enormlz (ei)            left-justify the significand
 897  *      eshdn1 (ai)             shift significand and guards down 1 bit
 898  *      eshdn8 (ai)             shift down 8 bits
 899  *      eshdn6 (ai)             shift down 16 bits
 900  *      eshift (ai, n)          shift ai n bits up (or down if n < 0)
 901  *      eshup1 (ai)             shift significand and guards up 1 bit
 902  *      eshup8 (ai)             shift up 8 bits
 903  *      eshup6 (ai)             shift up 16 bits
 904  *      esubm (ai, bi)          subtract significands, bi = bi - ai
 905  *
 906  *
 907  * The result is always normalized and rounded to NI-4 word precision
 908  * after each arithmetic operation.
 909  *
 910  * Exception flags and NaNs are NOT fully supported.
 911  * This arithmetic should never produce a NaN output, but it might
 912  * be confused by a NaN input.
 913  * Define INFINITY in mconf.h for support of infinity; otherwise a
 914  * saturation arithmetic is implemented.
 915  * Denormals are always supported here where appropriate (e.g., not
 916  * for conversion to DEC numbers).
 917  *
 918  */
 919
 920 /*
 921  * Revision history:
 922  *
 923  *  5 Jan 84    PDP-11 assembly language version
 924  *  2 Mar 86    fixed bug in asctoq
 925  *  6 Dec 86    C language version
 926  * 30 Aug 88    100 digit version, improved rounding
 927  * 15 May 92    80-bit long double support
 928  *
 929  * Author:  S. L. Moshier.
 930  */
 931
 932
 933 /*                                                      mconf.h
 934  *
 935  *      Common include file for math routines
 936  *
 937  *
 938  *
 939  * SYNOPSIS:
 940  *
 941  * #include "mconf.h"
 942  *
 943  *
 944  *
 945  * DESCRIPTION:
 946  *
 947  * This file contains definitions for error codes that are
 948  * passed to the common error handling routine mtherr
 949  * (which see).
 950  *
 951  * The file also includes a conditional assembly definition
 952  * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
 953  * IEEE, or UNKnown).
 954  *
 955  * For Digital Equipment PDP-11 and VAX computers, certain
 956  * IBM systems, and others that use numbers with a 56-bit
 957  * significand, the symbol DEC should be defined.  In this
 958  * mode, most floating point constants are given as arrays
 959  * of octal integers to eliminate decimal to binary conversion
 960  * errors that might be introduced by the compiler.
 961  *
 962  * For computers, such as IBM PC, that follow the IEEE
 963  * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
 964  * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
 965  * These numbers have 53-bit significands.  In this mode, constants
 966  * are provided as arrays of hexadecimal 16 bit integers.
 967  *
 968  * To accommodate other types of computer arithmetic, all
 969  * constants are also provided in a normal decimal radix
 970  * which one can hope are correctly converted to a suitable
 971  * format by the available C language compiler.  To invoke
 972  * this mode, the symbol UNK is defined.
 973  *
 974  * An important difference among these modes is a predefined
 975  * set of machine arithmetic constants for each.  The numbers
 976  * MACHEP (the machine roundoff error), MAXNUM (largest number
 977  * represented), and several other parameters are preset by
 978  * the configuration symbol.  Check the file const.c to
 979  * ensure that these values are correct for your computer.
 980  *
 981  * For ANSI C compatibility, define ANSIC equal to 1.  Currently
 982  * this affects only the atan2 function and others that use it.
 983  */
 984 \f
 985 /* Constant definitions for math error conditions.  */
 986
 987 #define DOMAIN          1       /* argument domain error */
 988 #define SING            2       /* argument singularity */
 989 #define OVERFLOW        3       /* overflow range error */
 990 #define UNDERFLOW       4       /* underflow range error */
 991 #define TLOSS           5       /* total loss of precision */
 992 #define PLOSS           6       /* partial loss of precision */
 993
 994 /*  e type constants used by high precision check routines */
 995
 996 /*include "ehead.h"*/
 997 /* 0.0 */
 998 unsigned EMUSHORT ezero[NE] =
 999 {
1000   0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1001 extern unsigned EMUSHORT ezero[];
1002
1003 /* 5.0E-1 */
1004 unsigned EMUSHORT ehalf[NE] =
1005 {
1006   0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1007 extern unsigned EMUSHORT ehalf[];
1008
1009 /* 1.0E0 */
1010 unsigned EMUSHORT eone[NE] =
1011 {
1012   0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1013 extern unsigned EMUSHORT eone[];
1014
1015 /* 2.0E0 */
1016 unsigned EMUSHORT etwo[NE] =
1017 {
1018   0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1019 extern unsigned EMUSHORT etwo[];
1020
1021 /* 3.2E1 */
1022 unsigned EMUSHORT e32[NE] =
1023 {
1024   0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1025 extern unsigned EMUSHORT e32[];
1026
1027 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1028 unsigned EMUSHORT elog2[NE] =
1029 {
1030   0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1031 extern unsigned EMUSHORT elog2[];
1032
1033 /* 1.41421356237309504880168872420969807856967187537695E0 */
1034 unsigned EMUSHORT esqrt2[NE] =
1035 {
1036   0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1037 extern unsigned EMUSHORT esqrt2[];
1038
1039 /* 2/sqrt (PI) =
1040  * 1.12837916709551257389615890312154517168810125865800E0 */
1041 unsigned EMUSHORT eoneopi[NE] =
1042 {
1043   0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1044 extern unsigned EMUSHORT eoneopi[];
1045
1046 /* 3.14159265358979323846264338327950288419716939937511E0 */
1047 unsigned EMUSHORT epi[NE] =
1048 {
1049   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1050 extern unsigned EMUSHORT epi[];
1051
1052 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1053 unsigned EMUSHORT eeul[NE] =
1054 {
1055   0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1056 extern unsigned EMUSHORT eeul[];
1057
1058 /*
1059 include "ehead.h"
1060 include "mconf.h"
1061 */
1062
1063
1064
1065 /* Control register for rounding precision.
1066  * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1067  */
1068 int rndprc = NBITS;
1069 extern int rndprc;
1070
1071 void eaddm (), esubm (), emdnorm (), asctoeg ();
1072 static void toe24 (), toe53 (), toe64 ();
1073 void eremain (), einit (), eiremain ();
1074 int ecmpm (), edivm (), emulm ();
1075 void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1076 void etodec (), todec (), dectoe ();
1077
1078
1079
1080
1081 void
1082 einit ()
1083 {
1084 }
1085
1086 /*
1087 ; Clear out entire external format number.
1088 ;
1089 ; unsigned EMUSHORT x[];
1090 ; eclear (x);
1091 */
1092
1093 void
1094 eclear (x)
1095      register unsigned EMUSHORT *x;
1096 {
1097   register int i;
1098
1099   for (i = 0; i < NE; i++)
1100     *x++ = 0;
1101 }
1102
1103
1104
1105 /* Move external format number from a to b.
1106  *
1107  * emov (a, b);
1108  */
1109
1110 void
1111 emov (a, b)
1112      register unsigned EMUSHORT *a, *b;
1113 {
1114   register int i;
1115
1116   for (i = 0; i < NE; i++)
1117     *b++ = *a++;
1118 }
1119
1120
1121 /*
1122 ;       Absolute value of external format number
1123 ;
1124 ;       EMUSHORT x[NE];
1125 ;       eabs (x);
1126 */
1127
1128 void
1129 eabs (x)
1130      unsigned EMUSHORT x[];     /* x is the memory address of a short */
1131 {
1132
1133   x[NE - 1] &= 0x7fff;          /* sign is top bit of last word of external format */
1134 }
1135
1136
1137
1138
1139 /*
1140 ;       Negate external format number
1141 ;
1142 ;       unsigned EMUSHORT x[NE];
1143 ;       eneg (x);
1144 */
1145
1146 void
1147 eneg (x)
1148      unsigned EMUSHORT x[];
1149 {
1150
1151   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1152 }
1153
1154
1155
1156 /* Return 1 if external format number is negative,
1157  * else return zero.
1158  */
1159 int
1160 eisneg (x)
1161      unsigned EMUSHORT x[];
1162 {
1163
1164   if (x[NE - 1] & 0x8000)
1165     return (1);
1166   else
1167     return (0);
1168 }
1169
1170
1171 /* Return 1 if external format number has maximum possible exponent,
1172  * else return zero.
1173  */
1174 int
1175 eisinf (x)
1176      unsigned EMUSHORT x[];
1177 {
1178
1179   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1180     return (1);
1181   else
1182     return (0);
1183 }
1184
1185
1186 /*
1187 ; Fill entire number, including exponent and significand, with
1188 ; largest possible number.  These programs implement a saturation
1189 ; value that is an ordinary, legal number.  A special value
1190 ; "infinity" may also be implemented; this would require tests
1191 ; for that value and implementation of special rules for arithmetic
1192 ; operations involving inifinity.
1193 */
1194
1195 void
1196 einfin (x)
1197      register unsigned EMUSHORT *x;
1198 {
1199   register int i;
1200
1201 #ifdef INFINITY
1202   for (i = 0; i < NE - 1; i++)
1203     *x++ = 0;
1204   *x |= 32767;
1205 #else
1206   for (i = 0; i < NE - 1; i++)
1207     *x++ = 0xffff;
1208   *x |= 32766;
1209   if (rndprc < NBITS)
1210     {
1211       if (rndprc == 64)
1212         {
1213           *(x - 5) = 0;
1214         }
1215       if (rndprc == 53)
1216         {
1217           *(x - 4) = 0xf800;
1218         }
1219       else
1220         {
1221           *(x - 4) = 0;
1222           *(x - 3) = 0;
1223           *(x - 2) = 0xff00;
1224         }
1225     }
1226 #endif
1227 }
1228
1229
1230
1231 /* Move in external format number,
1232  * converting it to internal format.
1233  */
1234 void
1235 emovi (a, b)
1236      unsigned EMUSHORT *a, *b;
1237 {
1238   register unsigned EMUSHORT *p, *q;
1239   int i;
1240
1241   q = b;
1242   p = a + (NE - 1);             /* point to last word of external number */
1243   /* get the sign bit */
1244   if (*p & 0x8000)
1245     *q++ = 0xffff;
1246   else
1247     *q++ = 0;
1248   /* get the exponent */
1249   *q = *p--;
1250   *q++ &= 0x7fff;               /* delete the sign bit */
1251 #ifdef INFINITY
1252   if ((*(q - 1) & 0x7fff) == 0x7fff)
1253     {
1254       for (i = 2; i < NI; i++)
1255         *q++ = 0;
1256       return;
1257     }
1258 #endif
1259   /* clear high guard word */
1260   *q++ = 0;
1261   /* move in the significand */
1262   for (i = 0; i < NE - 1; i++)
1263     *q++ = *p--;
1264   /* clear low guard word */
1265   *q = 0;
1266 }
1267
1268
1269 /* Move internal format number out,
1270  * converting it to external format.
1271  */
1272 void
1273 emovo (a, b)
1274      unsigned EMUSHORT *a, *b;
1275 {
1276   register unsigned EMUSHORT *p, *q;
1277   unsigned EMUSHORT i;
1278
1279   p = a;
1280   q = b + (NE - 1);             /* point to output exponent */
1281   /* combine sign and exponent */
1282   i = *p++;
1283   if (i)
1284     *q-- = *p++ | 0x8000;
1285   else
1286     *q-- = *p++;
1287 #ifdef INFINITY
1288   if (*(p - 1) == 0x7fff)
1289     {
1290       einfin (b);
1291       return;
1292     }
1293 #endif
1294   /* skip over guard word */
1295   ++p;
1296   /* move the significand */
1297   for (i = 0; i < NE - 1; i++)
1298     *q-- = *p++;
1299 }
1300
1301
1302
1303
1304 /* Clear out internal format number.
1305  */
1306
1307 void
1308 ecleaz (xi)
1309      register unsigned EMUSHORT *xi;
1310 {
1311   register int i;
1312
1313   for (i = 0; i < NI; i++)
1314     *xi++ = 0;
1315 }
1316
1317
1318 /* same, but don't touch the sign. */
1319
1320 void
1321 ecleazs (xi)
1322      register unsigned EMUSHORT *xi;
1323 {
1324   register int i;
1325
1326   ++xi;
1327   for (i = 0; i < NI - 1; i++)
1328     *xi++ = 0;
1329 }
1330
1331
1332
1333 /* Move internal format number from a to b.
1334  */
1335 void
1336 emovz (a, b)
1337      register unsigned EMUSHORT *a, *b;
1338 {
1339   register int i;
1340
1341   for (i = 0; i < NI - 1; i++)
1342     *b++ = *a++;
1343   /* clear low guard word */
1344   *b = 0;
1345 }
1346
1347
1348 /*
1349 ;       Compare significands of numbers in internal format.
1350 ;       Guard words are included in the comparison.
1351 ;
1352 ;       unsigned EMUSHORT a[NI], b[NI];
1353 ;       cmpm (a, b);
1354 ;
1355 ;       for the significands:
1356 ;       returns +1 if a > b
1357 ;                0 if a == b
1358 ;               -1 if a < b
1359 */
1360 int
1361 ecmpm (a, b)
1362      register unsigned EMUSHORT *a, *b;
1363 {
1364   int i;
1365
1366   a += M;                       /* skip up to significand area */
1367   b += M;
1368   for (i = M; i < NI; i++)
1369     {
1370       if (*a++ != *b++)
1371         goto difrnt;
1372     }
1373   return (0);
1374
1375  difrnt:
1376   if (*(--a) > *(--b))
1377     return (1);
1378   else
1379     return (-1);
1380 }
1381
1382
1383 /*
1384 ;       Shift significand down by 1 bit
1385 */
1386
1387 void
1388 eshdn1 (x)
1389      register unsigned EMUSHORT *x;
1390 {
1391   register unsigned EMUSHORT bits;
1392   int i;
1393
1394   x += M;                       /* point to significand area */
1395
1396   bits = 0;
1397   for (i = M; i < NI; i++)
1398     {
1399       if (*x & 1)
1400         bits |= 1;
1401       *x >>= 1;
1402       if (bits & 2)
1403         *x |= 0x8000;
1404       bits <<= 1;
1405       ++x;
1406     }
1407 }
1408
1409
1410
1411 /*
1412 ;       Shift significand up by 1 bit
1413 */
1414
1415 void
1416 eshup1 (x)
1417      register unsigned EMUSHORT *x;
1418 {
1419   register unsigned EMUSHORT bits;
1420   int i;
1421
1422   x += NI - 1;
1423   bits = 0;
1424
1425   for (i = M; i < NI; i++)
1426     {
1427       if (*x & 0x8000)
1428         bits |= 1;
1429       *x <<= 1;
1430       if (bits & 2)
1431         *x |= 1;
1432       bits <<= 1;
1433       --x;
1434     }
1435 }
1436
1437
1438
1439 /*
1440 ;       Shift significand down by 8 bits
1441 */
1442
1443 void
1444 eshdn8 (x)
1445      register unsigned EMUSHORT *x;
1446 {
1447   register unsigned EMUSHORT newbyt, oldbyt;
1448   int i;
1449
1450   x += M;
1451   oldbyt = 0;
1452   for (i = M; i < NI; i++)
1453     {
1454       newbyt = *x << 8;
1455       *x >>= 8;
1456       *x |= oldbyt;
1457       oldbyt = newbyt;
1458       ++x;
1459     }
1460 }
1461
1462 /*
1463 ;       Shift significand up by 8 bits
1464 */
1465
1466 void
1467 eshup8 (x)
1468      register unsigned EMUSHORT *x;
1469 {
1470   int i;
1471   register unsigned EMUSHORT newbyt, oldbyt;
1472
1473   x += NI - 1;
1474   oldbyt = 0;
1475
1476   for (i = M; i < NI; i++)
1477     {
1478       newbyt = *x >> 8;
1479       *x <<= 8;
1480       *x |= oldbyt;
1481       oldbyt = newbyt;
1482       --x;
1483     }
1484 }
1485
1486 /*
1487 ;       Shift significand up by 16 bits
1488 */
1489
1490 void
1491 eshup6 (x)
1492      register unsigned EMUSHORT *x;
1493 {
1494   int i;
1495   register unsigned EMUSHORT *p;
1496
1497   p = x + M;
1498   x += M + 1;
1499
1500   for (i = M; i < NI - 1; i++)
1501     *p++ = *x++;
1502
1503   *p = 0;
1504 }
1505
1506 /*
1507 ;       Shift significand down by 16 bits
1508 */
1509
1510 void
1511 eshdn6 (x)
1512      register unsigned EMUSHORT *x;
1513 {
1514   int i;
1515   register unsigned EMUSHORT *p;
1516
1517   x += NI - 1;
1518   p = x + 1;
1519
1520   for (i = M; i < NI - 1; i++)
1521     *(--p) = *(--x);
1522
1523   *(--p) = 0;
1524 }
1525 \f
1526 /*
1527 ;       Add significands
1528 ;       x + y replaces y
1529 */
1530
1531 void
1532 eaddm (x, y)
1533      unsigned EMUSHORT *x, *y;
1534 {
1535   register unsigned EMULONG a;
1536   int i;
1537   unsigned int carry;
1538
1539   x += NI - 1;
1540   y += NI - 1;
1541   carry = 0;
1542   for (i = M; i < NI; i++)
1543     {
1544       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1545       if (a & 0x10000)
1546         carry = 1;
1547       else
1548         carry = 0;
1549       *y = (unsigned EMUSHORT) a;
1550       --x;
1551       --y;
1552     }
1553 }
1554
1555 /*
1556 ;       Subtract significands
1557 ;       y - x replaces y
1558 */
1559
1560 void
1561 esubm (x, y)
1562      unsigned EMUSHORT *x, *y;
1563 {
1564   unsigned EMULONG a;
1565   int i;
1566   unsigned int carry;
1567
1568   x += NI - 1;
1569   y += NI - 1;
1570   carry = 0;
1571   for (i = M; i < NI; i++)
1572     {
1573       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1574       if (a & 0x10000)
1575         carry = 1;
1576       else
1577         carry = 0;
1578       *y = (unsigned EMUSHORT) a;
1579       --x;
1580       --y;
1581     }
1582 }
1583
1584
1585 /* Divide significands */
1586
1587 static unsigned EMUSHORT equot[NI];
1588
1589 int
1590 edivm (den, num)
1591      unsigned EMUSHORT den[], num[];
1592 {
1593   int i;
1594   register unsigned EMUSHORT *p, *q;
1595   unsigned EMUSHORT j;
1596
1597   p = &equot[0];
1598   *p++ = num[0];
1599   *p++ = num[1];
1600
1601   for (i = M; i < NI; i++)
1602     {
1603       *p++ = 0;
1604     }
1605
1606   /* Use faster compare and subtraction if denominator
1607    * has only 15 bits of significance.
1608    */
1609   p = &den[M + 2];
1610   if (*p++ == 0)
1611     {
1612       for (i = M + 3; i < NI; i++)
1613         {
1614           if (*p++ != 0)
1615             goto fulldiv;
1616         }
1617       if ((den[M + 1] & 1) != 0)
1618         goto fulldiv;
1619       eshdn1 (num);
1620       eshdn1 (den);
1621
1622       p = &den[M + 1];
1623       q = &num[M + 1];
1624
1625       for (i = 0; i < NBITS + 2; i++)
1626         {
1627           if (*p <= *q)
1628             {
1629               *q -= *p;
1630               j = 1;
1631             }
1632           else
1633             {
1634               j = 0;
1635             }
1636           eshup1 (equot);
1637           equot[NI - 2] |= j;
1638           eshup1 (num);
1639         }
1640       goto divdon;
1641     }
1642
1643   /* The number of quotient bits to calculate is
1644    * NBITS + 1 scaling guard bit + 1 roundoff bit.
1645    */
1646  fulldiv:
1647
1648   p = &equot[NI - 2];
1649   for (i = 0; i < NBITS + 2; i++)
1650     {
1651       if (ecmpm (den, num) <= 0)
1652         {
1653           esubm (den, num);
1654           j = 1;                /* quotient bit = 1 */
1655         }
1656       else
1657         j = 0;
1658       eshup1 (equot);
1659       *p |= j;
1660       eshup1 (num);
1661     }
1662
1663  divdon:
1664
1665   eshdn1 (equot);
1666   eshdn1 (equot);
1667
1668   /* test for nonzero remainder after roundoff bit */
1669   p = &num[M];
1670   j = 0;
1671   for (i = M; i < NI; i++)
1672     {
1673       j |= *p++;
1674     }
1675   if (j)
1676     j = 1;
1677
1678
1679   for (i = 0; i < NI; i++)
1680     num[i] = equot[i];
1681   return ((int) j);
1682 }
1683
1684
1685 /* Multiply significands */
1686 int
1687 emulm (a, b)
1688      unsigned EMUSHORT a[], b[];
1689 {
1690   unsigned EMUSHORT *p, *q;
1691   int i, j, k;
1692
1693   equot[0] = b[0];
1694   equot[1] = b[1];
1695   for (i = M; i < NI; i++)
1696     equot[i] = 0;
1697
1698   p = &a[NI - 2];
1699   k = NBITS;
1700   while (*p == 0)               /* significand is not supposed to be all zero */
1701     {
1702       eshdn6 (a);
1703       k -= 16;
1704     }
1705   if ((*p & 0xff) == 0)
1706     {
1707       eshdn8 (a);
1708       k -= 8;
1709     }
1710
1711   q = &equot[NI - 1];
1712   j = 0;
1713   for (i = 0; i < k; i++)
1714     {
1715       if (*p & 1)
1716         eaddm (b, equot);
1717       /* remember if there were any nonzero bits shifted out */
1718       if (*q & 1)
1719         j |= 1;
1720       eshdn1 (a);
1721       eshdn1 (equot);
1722     }
1723
1724   for (i = 0; i < NI; i++)
1725     b[i] = equot[i];
1726
1727   /* return flag for lost nonzero bits */
1728   return (j);
1729 }
1730
1731
1732
1733 /*
1734  * Normalize and round off.
1735  *
1736  * The internal format number to be rounded is "s".
1737  * Input "lost" indicates whether or not the number is exact.
1738  * This is the so-called sticky bit.
1739  *
1740  * Input "subflg" indicates whether the number was obtained
1741  * by a subtraction operation.  In that case if lost is nonzero
1742  * then the number is slightly smaller than indicated.
1743  *
1744  * Input "exp" is the biased exponent, which may be negative.
1745  * the exponent field of "s" is ignored but is replaced by
1746  * "exp" as adjusted by normalization and rounding.
1747  *
1748  * Input "rcntrl" is the rounding control.
1749  */
1750
1751 static int rlast = -1;
1752 static int rw = 0;
1753 static unsigned EMUSHORT rmsk = 0;
1754 static unsigned EMUSHORT rmbit = 0;
1755 static unsigned EMUSHORT rebit = 0;
1756 static int re = 0;
1757 static unsigned EMUSHORT rbit[NI];
1758
1759 void
1760 emdnorm (s, lost, subflg, exp, rcntrl)
1761      unsigned EMUSHORT s[];
1762      int lost;
1763      int subflg;
1764      EMULONG exp;
1765      int rcntrl;
1766 {
1767   int i, j;
1768   unsigned EMUSHORT r;
1769
1770   /* Normalize */
1771   j = enormlz (s);
1772
1773   /* a blank significand could mean either zero or infinity. */
1774 #ifndef INFINITY
1775   if (j > NBITS)
1776     {
1777       ecleazs (s);
1778       return;
1779     }
1780 #endif
1781   exp -= j;
1782 #ifndef INFINITY
1783   if (exp >= 32767L)
1784     goto overf;
1785 #else
1786   if ((j > NBITS) && (exp < 32767))
1787     {
1788       ecleazs (s);
1789       return;
1790     }
1791 #endif
1792   if (exp < 0L)
1793     {
1794       if (exp > (EMULONG) (-NBITS - 1))
1795         {
1796           j = (int) exp;
1797           i = eshift (s, j);
1798           if (i)
1799             lost = 1;
1800         }
1801       else
1802         {
1803           ecleazs (s);
1804           return;
1805         }
1806     }
1807   /* Round off, unless told not to by rcntrl. */
1808   if (rcntrl == 0)
1809     goto mdfin;
1810   /* Set up rounding parameters if the control register changed. */
1811   if (rndprc != rlast)
1812     {
1813       ecleaz (rbit);
1814       switch (rndprc)
1815         {
1816         default:
1817         case NBITS:
1818           rw = NI - 1;          /* low guard word */
1819           rmsk = 0xffff;
1820           rmbit = 0x8000;
1821           rbit[rw - 1] = 1;
1822           re = NI - 2;
1823           rebit = 1;
1824           break;
1825         case 64:
1826           rw = 7;
1827           rmsk = 0xffff;
1828           rmbit = 0x8000;
1829           rbit[rw - 1] = 1;
1830           re = rw - 1;
1831           rebit = 1;
1832           break;
1833           /* For DEC arithmetic */
1834         case 56:
1835           rw = 6;
1836           rmsk = 0xff;
1837           rmbit = 0x80;
1838           rbit[rw] = 0x100;
1839           re = rw;
1840           rebit = 0x100;
1841           break;
1842         case 53:
1843           rw = 6;
1844           rmsk = 0x7ff;
1845           rmbit = 0x0400;
1846           rbit[rw] = 0x800;
1847           re = rw;
1848           rebit = 0x800;
1849           break;
1850         case 24:
1851           rw = 4;
1852           rmsk = 0xff;
1853           rmbit = 0x80;
1854           rbit[rw] = 0x100;
1855           re = rw;
1856           rebit = 0x100;
1857           break;
1858         }
1859       rlast = rndprc;
1860     }
1861
1862   if (rndprc >= 64)
1863     {
1864       r = s[rw] & rmsk;
1865       if (rndprc == 64)
1866         {
1867           i = rw + 1;
1868           while (i < NI)
1869             {
1870               if (s[i])
1871                 r |= 1;
1872               s[i] = 0;
1873               ++i;
1874             }
1875         }
1876     }
1877   else
1878     {
1879       if (exp <= 0)
1880         eshdn1 (s);
1881       r = s[rw] & rmsk;
1882       /* These tests assume NI = 8 */
1883       i = rw + 1;
1884       while (i < NI)
1885         {
1886           if (s[i])
1887             r |= 1;
1888           s[i] = 0;
1889           ++i;
1890         }
1891       /*
1892          if (rndprc == 24)
1893          {
1894          if (s[5] || s[6])
1895          r |= 1;
1896          s[5] = 0;
1897          s[6] = 0;
1898          }
1899          */
1900       s[rw] &= ~rmsk;
1901     }
1902   if ((r & rmbit) != 0)
1903     {
1904       if (r == rmbit)
1905         {
1906           if (lost == 0)
1907             {                   /* round to even */
1908               if ((s[re] & rebit) == 0)
1909                 goto mddone;
1910             }
1911           else
1912             {
1913               if (subflg != 0)
1914                 goto mddone;
1915             }
1916         }
1917       eaddm (rbit, s);
1918     }
1919  mddone:
1920   if ((rndprc < 64) && (exp <= 0))
1921     {
1922       eshup1 (s);
1923     }
1924   if (s[2] != 0)
1925     {                           /* overflow on roundoff */
1926       eshdn1 (s);
1927       exp += 1;
1928     }
1929  mdfin:
1930   s[NI - 1] = 0;
1931   if (exp >= 32767L)
1932     {
1933 #ifndef INFINITY
1934     overf:
1935 #endif
1936 #ifdef INFINITY
1937       s[1] = 32767;
1938       for (i = 2; i < NI - 1; i++)
1939         s[i] = 0;
1940       if (extra_warnings)
1941         warning ("floating point overflow");
1942 #else
1943       s[1] = 32766;
1944       s[2] = 0;
1945       for (i = M + 1; i < NI - 1; i++)
1946         s[i] = 0xffff;
1947       s[NI - 1] = 0;
1948       if (rndprc < 64)
1949         {
1950           s[rw] &= ~rmsk;
1951           if (rndprc == 24)
1952             {
1953               s[5] = 0;
1954               s[6] = 0;
1955             }
1956         }
1957 #endif
1958       return;
1959     }
1960   if (exp < 0)
1961     s[1] = 0;
1962   else
1963     s[1] = (unsigned EMUSHORT) exp;
1964 }
1965
1966
1967
1968 /*
1969 ;       Subtract external format numbers.
1970 ;
1971 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1972 ;       esub (a, b, c);  c = b - a
1973 */
1974
1975 static int subflg = 0;
1976
1977 void
1978 esub (a, b, c)
1979      unsigned EMUSHORT *a, *b, *c;
1980 {
1981
1982   subflg = 1;
1983   eadd1 (a, b, c);
1984 }
1985
1986
1987 /*
1988 ;       Add.
1989 ;
1990 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1991 ;       eadd (a, b, c);  c = b + a
1992 */
1993 void
1994 eadd (a, b, c)
1995      unsigned EMUSHORT *a, *b, *c;
1996 {
1997
1998   subflg = 0;
1999   eadd1 (a, b, c);
2000 }
2001
2002 void
2003 eadd1 (a, b, c)
2004      unsigned EMUSHORT *a, *b, *c;
2005 {
2006   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2007   int i, lost, j, k;
2008   EMULONG lt, lta, ltb;
2009
2010 #ifdef INFINITY
2011   if (eisinf (a))
2012     {
2013       emov (a, c);
2014       if (subflg)
2015         eneg (c);
2016       return;
2017     }
2018   if (eisinf (b))
2019     {
2020       emov (b, c);
2021       return;
2022     }
2023 #endif
2024   emovi (a, ai);
2025   emovi (b, bi);
2026   if (subflg)
2027     ai[0] = ~ai[0];
2028
2029   /* compare exponents */
2030   lta = ai[E];
2031   ltb = bi[E];
2032   lt = lta - ltb;
2033   if (lt > 0L)
2034     {                           /* put the larger number in bi */
2035       emovz (bi, ci);
2036       emovz (ai, bi);
2037       emovz (ci, ai);
2038       ltb = bi[E];
2039       lt = -lt;
2040     }
2041   lost = 0;
2042   if (lt != 0L)
2043     {
2044       if (lt < (EMULONG) (-NBITS - 1))
2045         goto done;              /* answer same as larger addend */
2046       k = (int) lt;
2047       lost = eshift (ai, k);    /* shift the smaller number down */
2048     }
2049   else
2050     {
2051       /* exponents were the same, so must compare significands */
2052       i = ecmpm (ai, bi);
2053       if (i == 0)
2054         {                       /* the numbers are identical in magnitude */
2055           /* if different signs, result is zero */
2056           if (ai[0] != bi[0])
2057             {
2058               eclear (c);
2059               return;
2060             }
2061           /* if same sign, result is double */
2062           /* double denomalized tiny number */
2063           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2064             {
2065               eshup1 (bi);
2066               goto done;
2067             }
2068           /* add 1 to exponent unless both are zero! */
2069           for (j = 1; j < NI - 1; j++)
2070             {
2071               if (bi[j] != 0)
2072                 {
2073                   /* This could overflow, but let emovo take care of that. */
2074                   ltb += 1;
2075                   break;
2076                 }
2077             }
2078           bi[E] = (unsigned EMUSHORT) ltb;
2079           goto done;
2080         }
2081       if (i > 0)
2082         {                       /* put the larger number in bi */
2083           emovz (bi, ci);
2084           emovz (ai, bi);
2085           emovz (ci, ai);
2086         }
2087     }
2088   if (ai[0] == bi[0])
2089     {
2090       eaddm (ai, bi);
2091       subflg = 0;
2092     }
2093   else
2094     {
2095       esubm (ai, bi);
2096       subflg = 1;
2097     }
2098   emdnorm (bi, lost, subflg, ltb, 64);
2099
2100  done:
2101   emovo (bi, c);
2102 }
2103
2104
2105
2106 /*
2107 ;       Divide.
2108 ;
2109 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2110 ;       ediv (a, b, c); c = b / a
2111 */
2112 void
2113 ediv (a, b, c)
2114      unsigned EMUSHORT *a, *b, *c;
2115 {
2116   unsigned EMUSHORT ai[NI], bi[NI];
2117   int i;
2118   EMULONG lt, lta, ltb;
2119
2120 #ifdef INFINITY
2121   if (eisinf (b))
2122     {
2123       if (eisneg (a) ^ eisneg (b))
2124         *(c + (NE - 1)) = 0x8000;
2125       else
2126         *(c + (NE - 1)) = 0;
2127       einfin (c);
2128       return;
2129     }
2130   if (eisinf (a))
2131     {
2132       eclear (c);
2133       return;
2134     }
2135 #endif
2136   emovi (a, ai);
2137   emovi (b, bi);
2138   lta = ai[E];
2139   ltb = bi[E];
2140   if (bi[E] == 0)
2141     {                           /* See if numerator is zero. */
2142       for (i = 1; i < NI - 1; i++)
2143         {
2144           if (bi[i] != 0)
2145             {
2146               ltb -= enormlz (bi);
2147               goto dnzro1;
2148             }
2149         }
2150       eclear (c);
2151       return;
2152     }
2153  dnzro1:
2154
2155   if (ai[E] == 0)
2156     {                           /* possible divide by zero */
2157       for (i = 1; i < NI - 1; i++)
2158         {
2159           if (ai[i] != 0)
2160             {
2161               lta -= enormlz (ai);
2162               goto dnzro2;
2163             }
2164         }
2165       if (ai[0] == bi[0])
2166         *(c + (NE - 1)) = 0;
2167       else
2168         *(c + (NE - 1)) = 0x8000;
2169       einfin (c);
2170       mtherr ("ediv", SING);
2171       return;
2172     }
2173  dnzro2:
2174
2175   i = edivm (ai, bi);
2176   /* calculate exponent */
2177   lt = ltb - lta + EXONE;
2178   emdnorm (bi, i, 0, lt, 64);
2179   /* set the sign */
2180   if (ai[0] == bi[0])
2181     bi[0] = 0;
2182   else
2183     bi[0] = 0Xffff;
2184   emovo (bi, c);
2185 }
2186
2187
2188
2189 /*
2190 ;       Multiply.
2191 ;
2192 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2193 ;       emul (a, b, c); c = b * a
2194 */
2195 void
2196 emul (a, b, c)
2197      unsigned EMUSHORT *a, *b, *c;
2198 {
2199   unsigned EMUSHORT ai[NI], bi[NI];
2200   int i, j;
2201   EMULONG lt, lta, ltb;
2202
2203 #ifdef INFINITY
2204   if (eisinf (a) || eisinf (b))
2205     {
2206       if (eisneg (a) ^ eisneg (b))
2207         *(c + (NE - 1)) = 0x8000;
2208       else
2209         *(c + (NE - 1)) = 0;
2210       einfin (c);
2211       return;
2212     }
2213 #endif
2214   emovi (a, ai);
2215   emovi (b, bi);
2216   lta = ai[E];
2217   ltb = bi[E];
2218   if (ai[E] == 0)
2219     {
2220       for (i = 1; i < NI - 1; i++)
2221         {
2222           if (ai[i] != 0)
2223             {
2224               lta -= enormlz (ai);
2225               goto mnzer1;
2226             }
2227         }
2228       eclear (c);
2229       return;
2230     }
2231  mnzer1:
2232
2233   if (bi[E] == 0)
2234     {
2235       for (i = 1; i < NI - 1; i++)
2236         {
2237           if (bi[i] != 0)
2238             {
2239               ltb -= enormlz (bi);
2240               goto mnzer2;
2241             }
2242         }
2243       eclear (c);
2244       return;
2245     }
2246  mnzer2:
2247
2248   /* Multiply significands */
2249   j = emulm (ai, bi);
2250   /* calculate exponent */
2251   lt = lta + ltb - (EXONE - 1);
2252   emdnorm (bi, j, 0, lt, 64);
2253   /* calculate sign of product */
2254   if (ai[0] == bi[0])
2255     bi[0] = 0;
2256   else
2257     bi[0] = 0xffff;
2258   emovo (bi, c);
2259 }
2260
2261
2262
2263
2264 /*
2265 ; Convert IEEE double precision to e type
2266 ;       double d;
2267 ;       unsigned EMUSHORT x[N+2];
2268 ;       e53toe (&d, x);
2269 */
2270 void
2271 e53toe (e, y)
2272      unsigned EMUSHORT *e, *y;
2273 {
2274 #ifdef DEC
2275
2276   dectoe (e, y);                /* see etodec.c */
2277
2278 #else
2279
2280   register unsigned EMUSHORT r;
2281   register unsigned EMUSHORT *p;
2282   unsigned EMUSHORT yy[NI];
2283   int denorm, k;
2284
2285   denorm = 0;                   /* flag if denormalized number */
2286   ecleaz (yy);
2287 #ifdef IBMPC
2288   e += 3;
2289 #endif
2290   r = *e;
2291   yy[0] = 0;
2292   if (r & 0x8000)
2293     yy[0] = 0xffff;
2294   yy[M] = (r & 0x0f) | 0x10;
2295   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2296 #ifdef INFINITY
2297   if (r == 0x7ff0)
2298     {
2299       einfin (y);
2300       if (r & 0x8000)
2301         eneg (y);
2302       return;
2303     }
2304 #endif
2305   r >>= 4;
2306   /* If zero exponent, then the significand is denormalized.
2307    * So, take back the understood high significand bit. */
2308   if (r == 0)
2309     {
2310       denorm = 1;
2311       yy[M] &= ~0x10;
2312     }
2313   r += EXONE - 01777;
2314   yy[E] = r;
2315   p = &yy[M + 1];
2316 #ifdef IBMPC
2317   *p++ = *(--e);
2318   *p++ = *(--e);
2319   *p++ = *(--e);
2320 #endif
2321 #ifdef MIEEE
2322   ++e;
2323   *p++ = *e++;
2324   *p++ = *e++;
2325   *p++ = *e++;
2326 #endif
2327   eshift (yy, -5);
2328   if (denorm)
2329     {                           /* if zero exponent, then normalize the significand */
2330       if ((k = enormlz (yy)) > NBITS)
2331         ecleazs (yy);
2332       else
2333         yy[E] -= (unsigned EMUSHORT) (k - 1);
2334     }
2335   emovo (yy, y);
2336 #endif /* not DEC */
2337 }
2338
2339 void
2340 e64toe (e, y)
2341      unsigned EMUSHORT *e, *y;
2342 {
2343   unsigned EMUSHORT yy[NI];
2344   unsigned EMUSHORT *p, *q;
2345   int i;
2346
2347   p = yy;
2348   for (i = 0; i < NE - 5; i++)
2349     *p++ = 0;
2350 #ifdef IBMPC
2351   for (i = 0; i < 5; i++)
2352     *p++ = *e++;
2353 #endif
2354 #ifdef DEC
2355   for (i = 0; i < 5; i++)
2356     *p++ = *e++;
2357 #endif
2358 #ifdef MIEEE
2359   p = &yy[0] + (NE - 1);
2360   *p-- = *e++;
2361   ++e;
2362   for (i = 0; i < 4; i++)
2363     *p-- = *e++;
2364 #endif
2365   p = yy;
2366   q = y;
2367 #ifdef INFINITY
2368   if (*p == 0x7fff)
2369     {
2370       einfin (y);
2371       if (*p & 0x8000)
2372         eneg (y);
2373       return;
2374     }
2375 #endif
2376   for (i = 0; i < NE; i++)
2377     *q++ = *p++;
2378 }
2379
2380
2381 /*
2382 ; Convert IEEE single precision to e type
2383 ;       float d;
2384 ;       unsigned EMUSHORT x[N+2];
2385 ;       dtox (&d, x);
2386 */
2387 void
2388 e24toe (e, y)
2389      unsigned EMUSHORT *e, *y;
2390 {
2391   register unsigned EMUSHORT r;
2392   register unsigned EMUSHORT *p;
2393   unsigned EMUSHORT yy[NI];
2394   int denorm, k;
2395
2396   denorm = 0;                   /* flag if denormalized number */
2397   ecleaz (yy);
2398 #ifdef IBMPC
2399   e += 1;
2400 #endif
2401 #ifdef DEC
2402   e += 1;
2403 #endif
2404   r = *e;
2405   yy[0] = 0;
2406   if (r & 0x8000)
2407     yy[0] = 0xffff;
2408   yy[M] = (r & 0x7f) | 0200;
2409   r &= ~0x807f;                 /* strip sign and 7 significand bits */
2410 #ifdef INFINITY
2411   if (r == 0x7f80)
2412     {
2413       einfin (y);
2414       if (r & 0x8000)
2415         eneg (y);
2416       return;
2417     }
2418 #endif
2419   r >>= 7;
2420   /* If zero exponent, then the significand is denormalized.
2421    * So, take back the understood high significand bit. */
2422   if (r == 0)
2423     {
2424       denorm = 1;
2425       yy[M] &= ~0200;
2426     }
2427   r += EXONE - 0177;
2428   yy[E] = r;
2429   p = &yy[M + 1];
2430 #ifdef IBMPC
2431   *p++ = *(--e);
2432 #endif
2433 #ifdef DEC
2434   *p++ = *(--e);
2435 #endif
2436 #ifdef MIEEE
2437   ++e;
2438   *p++ = *e++;
2439 #endif
2440   eshift (yy, -8);
2441   if (denorm)
2442     {                           /* if zero exponent, then normalize the significand */
2443       if ((k = enormlz (yy)) > NBITS)
2444         ecleazs (yy);
2445       else
2446         yy[E] -= (unsigned EMUSHORT) (k - 1);
2447     }
2448   emovo (yy, y);
2449 }
2450
2451
2452 void
2453 etoe64 (x, e)
2454      unsigned EMUSHORT *x, *e;
2455 {
2456   unsigned EMUSHORT xi[NI];
2457   EMULONG exp;
2458   int rndsav;
2459
2460   emovi (x, xi);
2461   /* adjust exponent for offset */
2462   exp = (EMULONG) xi[E];
2463 #ifdef INFINITY
2464   if (eisinf (x))
2465     goto nonorm;
2466 #endif
2467   /* round off to nearest or even */
2468   rndsav = rndprc;
2469   rndprc = 64;
2470   emdnorm (xi, 0, 0, exp, 64);
2471   rndprc = rndsav;
2472  nonorm:
2473   toe64 (xi, e);
2474 }
2475
2476 /* move out internal format to ieee long double */
2477 static void
2478 toe64 (a, b)
2479      unsigned EMUSHORT *a, *b;
2480 {
2481   register unsigned EMUSHORT *p, *q;
2482   unsigned EMUSHORT i;
2483
2484   p = a;
2485 #ifdef MIEEE
2486   q = b;
2487 #else
2488   q = b + 4;                    /* point to output exponent */
2489 #if LONG_DOUBLE_TYPE_SIZE == 96
2490   /* Clear the last two bytes of 12-byte Intel format */
2491   *(q+1) = 0;
2492 #endif
2493 #endif
2494
2495   /* combine sign and exponent */
2496   i = *p++;
2497 #ifdef MIEEE
2498   if (i)
2499     *q++ = *p++ | 0x8000;
2500   else
2501     *q++ = *p++;
2502   *q++ = 0;
2503 #else
2504   if (i)
2505     *q-- = *p++ | 0x8000;
2506   else
2507     *q-- = *p++;
2508 #endif
2509   /* skip over guard word */
2510   ++p;
2511   /* move the significand */
2512 #ifdef MIEEE
2513   for (i = 0; i < 4; i++)
2514     *q++ = *p++;
2515 #else
2516   for (i = 0; i < 4; i++)
2517     *q-- = *p++;
2518 #endif
2519 }
2520
2521
2522 /*
2523 ; e type to IEEE double precision
2524 ;       double d;
2525 ;       unsigned EMUSHORT x[NE];
2526 ;       etoe53 (x, &d);
2527 */
2528
2529 #ifdef DEC
2530
2531 void
2532 etoe53 (x, e)
2533      unsigned EMUSHORT *x, *e;
2534 {
2535   etodec (x, e);                /* see etodec.c */
2536 }
2537
2538 static void
2539 toe53 (x, y)
2540      unsigned EMUSHORT *x, *y;
2541 {
2542   todec (x, y);
2543 }
2544
2545 #else
2546
2547 void
2548 etoe53 (x, e)
2549      unsigned EMUSHORT *x, *e;
2550 {
2551   unsigned EMUSHORT xi[NI];
2552   EMULONG exp;
2553   int rndsav;
2554
2555   emovi (x, xi);
2556   /* adjust exponent for offsets */
2557   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2558 #ifdef INFINITY
2559   if (eisinf (x))
2560     goto nonorm;
2561 #endif
2562   /* round off to nearest or even */
2563   rndsav = rndprc;
2564   rndprc = 53;
2565   emdnorm (xi, 0, 0, exp, 64);
2566   rndprc = rndsav;
2567  nonorm:
2568   toe53 (xi, e);
2569 }
2570
2571
2572 static void
2573 toe53 (x, y)
2574      unsigned EMUSHORT *x, *y;
2575 {
2576   unsigned EMUSHORT i;
2577   unsigned EMUSHORT *p;
2578
2579
2580   p = &x[0];
2581 #ifdef IBMPC
2582   y += 3;
2583 #endif
2584   *y = 0;                       /* output high order */
2585   if (*p++)
2586     *y = 0x8000;                /* output sign bit */
2587
2588   i = *p++;
2589   if (i >= (unsigned int) 2047)
2590     {                           /* Saturate at largest number less than infinity. */
2591 #ifdef INFINITY
2592       *y |= 0x7ff0;
2593 #ifdef IBMPC
2594       *(--y) = 0;
2595       *(--y) = 0;
2596       *(--y) = 0;
2597 #endif
2598 #ifdef MIEEE
2599       ++y;
2600       *y++ = 0;
2601       *y++ = 0;
2602       *y++ = 0;
2603 #endif
2604 #else
2605       *y |= (unsigned EMUSHORT) 0x7fef;
2606 #ifdef IBMPC
2607       *(--y) = 0xffff;
2608       *(--y) = 0xffff;
2609       *(--y) = 0xffff;
2610 #endif
2611 #ifdef MIEEE
2612       ++y;
2613       *y++ = 0xffff;
2614       *y++ = 0xffff;
2615       *y++ = 0xffff;
2616 #endif
2617 #endif
2618       return;
2619     }
2620   if (i == 0)
2621     {
2622       eshift (x, 4);
2623     }
2624   else
2625     {
2626       i <<= 4;
2627       eshift (x, 5);
2628     }
2629   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
2630   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
2631 #ifdef IBMPC
2632   *(--y) = *p++;
2633   *(--y) = *p++;
2634   *(--y) = *p;
2635 #endif
2636 #ifdef MIEEE
2637   ++y;
2638   *y++ = *p++;
2639   *y++ = *p++;
2640   *y++ = *p++;
2641 #endif
2642 }
2643
2644 #endif /* not DEC */
2645
2646
2647
2648 /*
2649 ; e type to IEEE single precision
2650 ;       float d;
2651 ;       unsigned EMUSHORT x[N+2];
2652 ;       xtod (x, &d);
2653 */
2654 void
2655 etoe24 (x, e)
2656      unsigned EMUSHORT *x, *e;
2657 {
2658   EMULONG exp;
2659   unsigned EMUSHORT xi[NI];
2660   int rndsav;
2661
2662   emovi (x, xi);
2663   /* adjust exponent for offsets */
2664   exp = (EMULONG) xi[E] - (EXONE - 0177);
2665 #ifdef INFINITY
2666   if (eisinf (x))
2667     goto nonorm;
2668 #endif
2669   /* round off to nearest or even */
2670   rndsav = rndprc;
2671   rndprc = 24;
2672   emdnorm (xi, 0, 0, exp, 64);
2673   rndprc = rndsav;
2674  nonorm:
2675   toe24 (xi, e);
2676 }
2677
2678 static void
2679 toe24 (x, y)
2680      unsigned EMUSHORT *x, *y;
2681 {
2682   unsigned EMUSHORT i;
2683   unsigned EMUSHORT *p;
2684
2685   p = &x[0];
2686 #ifdef IBMPC
2687   y += 1;
2688 #endif
2689 #ifdef DEC
2690   y += 1;
2691 #endif
2692   *y = 0;                       /* output high order */
2693   if (*p++)
2694     *y = 0x8000;                /* output sign bit */
2695
2696   i = *p++;
2697 /* Handle overflow cases. */
2698   if (i >= 255)
2699     {
2700 #ifdef INFINITY
2701       *y |= (unsigned EMUSHORT) 0x7f80;
2702 #ifdef IBMPC
2703       *(--y) = 0;
2704 #endif
2705 #ifdef DEC
2706       *(--y) = 0;
2707 #endif
2708 #ifdef MIEEE
2709       ++y;
2710       *y = 0;
2711 #endif
2712 #else  /* no INFINITY */
2713       *y |= (unsigned EMUSHORT) 0x7f7f;
2714 #ifdef IBMPC
2715       *(--y) = 0xffff;
2716 #endif
2717 #ifdef DEC
2718       *(--y) = 0xffff;
2719 #endif
2720 #ifdef MIEEE
2721       ++y;
2722       *y = 0xffff;
2723 #endif
2724 #ifdef ERANGE
2725       errno = ERANGE;
2726 #endif
2727 #endif  /* no INFINITY */
2728       return;
2729     }
2730   if (i == 0)
2731     {
2732       eshift (x, 7);
2733     }
2734   else
2735     {
2736       i <<= 7;
2737       eshift (x, 8);
2738     }
2739   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
2740   *y |= i;                      /* high order output already has sign bit set */
2741 #ifdef IBMPC
2742   *(--y) = *p;
2743 #endif
2744 #ifdef DEC
2745   *(--y) = *p;
2746 #endif
2747 #ifdef MIEEE
2748   ++y;
2749   *y = *p;
2750 #endif
2751 }
2752
2753
2754 /* Compare two e type numbers.
2755  *
2756  * unsigned EMUSHORT a[NE], b[NE];
2757  * ecmp (a, b);
2758  *
2759  *  returns +1 if a > b
2760  *           0 if a == b
2761  *          -1 if a < b
2762  */
2763 int
2764 ecmp (a, b)
2765      unsigned EMUSHORT *a, *b;
2766 {
2767   unsigned EMUSHORT ai[NI], bi[NI];
2768   register unsigned EMUSHORT *p, *q;
2769   register int i;
2770   int msign;
2771
2772   emovi (a, ai);
2773   p = ai;
2774   emovi (b, bi);
2775   q = bi;
2776
2777   if (*p != *q)
2778     {                           /* the signs are different */
2779       /* -0 equals + 0 */
2780       for (i = 1; i < NI - 1; i++)
2781         {
2782           if (ai[i] != 0)
2783             goto nzro;
2784           if (bi[i] != 0)
2785             goto nzro;
2786         }
2787       return (0);
2788     nzro:
2789       if (*p == 0)
2790         return (1);
2791       else
2792         return (-1);
2793     }
2794   /* both are the same sign */
2795   if (*p == 0)
2796     msign = 1;
2797   else
2798     msign = -1;
2799   i = NI - 1;
2800   do
2801     {
2802       if (*p++ != *q++)
2803         {
2804           goto diff;
2805         }
2806     }
2807   while (--i > 0);
2808
2809   return (0);                   /* equality */
2810
2811
2812
2813  diff:
2814
2815   if (*(--p) > *(--q))
2816     return (msign);             /* p is bigger */
2817   else
2818     return (-msign);            /* p is littler */
2819 }
2820
2821
2822
2823
2824 /* Find nearest integer to x = floor (x + 0.5)
2825  *
2826  * unsigned EMUSHORT x[NE], y[NE]
2827  * eround (x, y);
2828  */
2829 void
2830 eround (x, y)
2831      unsigned EMUSHORT *x, *y;
2832 {
2833   eadd (ehalf, x, y);
2834   efloor (y, y);
2835 }
2836
2837
2838
2839
2840 /*
2841 ; convert long integer to e type
2842 ;
2843 ;       long l;
2844 ;       unsigned EMUSHORT x[NE];
2845 ;       ltoe (&l, x);
2846 ; note &l is the memory address of l
2847 */
2848 void
2849 ltoe (lp, y)
2850      long *lp;                  /* lp is the memory address of a long integer */
2851      unsigned EMUSHORT *y;              /* y is the address of a short */
2852 {
2853   unsigned EMUSHORT yi[NI];
2854   unsigned long ll;
2855   int k;
2856
2857   ecleaz (yi);
2858   if (*lp < 0)
2859     {
2860       /* make it positive */
2861       ll = (unsigned long) (-(*lp));
2862       yi[0] = 0xffff;           /* put correct sign in the e type number */
2863     }
2864   else
2865     {
2866       ll = (unsigned long) (*lp);
2867     }
2868   /* move the long integer to yi significand area */
2869   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2870   yi[M + 1] = (unsigned EMUSHORT) ll;
2871
2872   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2873   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2874     ecleaz (yi);                /* it was zero */
2875   else
2876     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
2877   emovo (yi, y);                /* output the answer */
2878 }
2879
2880 /*
2881 ; convert unsigned long integer to e type
2882 ;
2883 ;       unsigned long l;
2884 ;       unsigned EMUSHORT x[NE];
2885 ;       ltox (&l, x);
2886 ; note &l is the memory address of l
2887 */
2888 void
2889 ultoe (lp, y)
2890      unsigned long *lp;         /* lp is the memory address of a long integer */
2891      unsigned EMUSHORT *y;              /* y is the address of a short */
2892 {
2893   unsigned EMUSHORT yi[NI];
2894   unsigned long ll;
2895   int k;
2896
2897   ecleaz (yi);
2898   ll = *lp;
2899
2900   /* move the long integer to ayi significand area */
2901   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2902   yi[M + 1] = (unsigned EMUSHORT) ll;
2903
2904   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2905   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2906     ecleaz (yi);                /* it was zero */
2907   else
2908     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
2909   emovo (yi, y);                /* output the answer */
2910 }
2911
2912
2913 /*
2914 ;       Find long integer and fractional parts
2915
2916 ;       long i;
2917 ;       unsigned EMUSHORT x[NE], frac[NE];
2918 ;       xifrac (x, &i, frac);
2919
2920   The integer output has the sign of the input.  The fraction is
2921 the positive fractional part of abs (x).
2922 */
2923 void
2924 eifrac (x, i, frac)
2925      unsigned EMUSHORT *x;
2926      long *i;
2927      unsigned EMUSHORT *frac;
2928 {
2929   unsigned EMUSHORT xi[NI];
2930   int k;
2931
2932   emovi (x, xi);
2933   k = (int) xi[E] - (EXONE - 1);
2934   if (k <= 0)
2935     {
2936       /* if exponent <= 0, integer = 0 and real output is fraction */
2937       *i = 0L;
2938       emovo (xi, frac);
2939       return;
2940     }
2941   if (k > (HOST_BITS_PER_LONG - 1))
2942     {
2943       /*
2944          ;      long integer overflow: output large integer
2945          ;      and correct fraction
2946          */
2947       if (xi[0])
2948         *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
2949       else
2950         *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
2951       eshift (xi, k);
2952       if (extra_warnings)
2953         warning ("overflow on truncation to integer");
2954       goto lab11;
2955     }
2956
2957   if (k > 16)
2958     {
2959       /*
2960          ; shift more than 16 bits: shift up k-16, output the integer,
2961          ; then complete the shift to get the fraction.
2962          */
2963       k -= 16;
2964       eshift (xi, k);
2965
2966       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
2967       eshup6 (xi);
2968       goto lab10;
2969     }
2970
2971   /* shift not more than 16 bits */
2972   eshift (xi, k);
2973   *i = (long) xi[M] & 0xffff;
2974
2975  lab10:
2976
2977   if (xi[0])
2978     *i = -(*i);
2979  lab11:
2980
2981   xi[0] = 0;
2982   xi[E] = EXONE - 1;
2983   xi[M] = 0;
2984   if ((k = enormlz (xi)) > NBITS)
2985     ecleaz (xi);
2986   else
2987     xi[E] -= (unsigned EMUSHORT) k;
2988
2989   emovo (xi, frac);
2990 }
2991
2992
2993 /*
2994 ;       Find unsigned long integer and fractional parts
2995
2996 ;       unsigned long i;
2997 ;       unsigned EMUSHORT x[NE], frac[NE];
2998 ;       xifrac (x, &i, frac);
2999
3000   A negative e type input yields integer output = 0
3001   but correct fraction.
3002 */
3003 void
3004 euifrac (x, i, frac)
3005      unsigned EMUSHORT *x;
3006      long *i;
3007      unsigned EMUSHORT *frac;
3008 {
3009   unsigned EMUSHORT xi[NI];
3010   int k;
3011
3012   emovi (x, xi);
3013   k = (int) xi[E] - (EXONE - 1);
3014   if (k <= 0)
3015     {
3016       /* if exponent <= 0, integer = 0 and argument is fraction */
3017       *i = 0L;
3018       emovo (xi, frac);
3019       return;
3020     }
3021   if (k > 32)
3022     {
3023       /*
3024          ;      long integer overflow: output large integer
3025          ;      and correct fraction
3026          */
3027       *i = ~(0L);
3028       eshift (xi, k);
3029       if (extra_warnings)
3030         warning ("overflow on truncation to unsigned integer");
3031       goto lab10;
3032     }
3033
3034   if (k > 16)
3035     {
3036       /*
3037          ; shift more than 16 bits: shift up k-16, output the integer,
3038          ; then complete the shift to get the fraction.
3039          */
3040       k -= 16;
3041       eshift (xi, k);
3042
3043       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
3044       eshup6 (xi);
3045       goto lab10;
3046     }
3047
3048   /* shift not more than 16 bits */
3049   eshift (xi, k);
3050   *i = (long) xi[M] & 0xffff;
3051
3052  lab10:
3053
3054   if (xi[0])
3055     *i = 0L;
3056
3057   xi[0] = 0;
3058   xi[E] = EXONE - 1;
3059   xi[M] = 0;
3060   if ((k = enormlz (xi)) > NBITS)
3061     ecleaz (xi);
3062   else
3063     xi[E] -= (unsigned EMUSHORT) k;
3064
3065   emovo (xi, frac);
3066 }
3067
3068
3069
3070 /*
3071 ;       Shift significand
3072 ;
3073 ;       Shifts significand area up or down by the number of bits
3074 ;       given by the variable sc.
3075 */
3076 int
3077 eshift (x, sc)
3078      unsigned EMUSHORT *x;
3079      int sc;
3080 {
3081   unsigned EMUSHORT lost;
3082   unsigned EMUSHORT *p;
3083
3084   if (sc == 0)
3085     return (0);
3086
3087   lost = 0;
3088   p = x + NI - 1;
3089
3090   if (sc < 0)
3091     {
3092       sc = -sc;
3093       while (sc >= 16)
3094         {
3095           lost |= *p;           /* remember lost bits */
3096           eshdn6 (x);
3097           sc -= 16;
3098         }
3099
3100       while (sc >= 8)
3101         {
3102           lost |= *p & 0xff;
3103           eshdn8 (x);
3104           sc -= 8;
3105         }
3106
3107       while (sc > 0)
3108         {
3109           lost |= *p & 1;
3110           eshdn1 (x);
3111           sc -= 1;
3112         }
3113     }
3114   else
3115     {
3116       while (sc >= 16)
3117         {
3118           eshup6 (x);
3119           sc -= 16;
3120         }
3121
3122       while (sc >= 8)
3123         {
3124           eshup8 (x);
3125           sc -= 8;
3126         }
3127
3128       while (sc > 0)
3129         {
3130           eshup1 (x);
3131           sc -= 1;
3132         }
3133     }
3134   if (lost)
3135     lost = 1;
3136   return ((int) lost);
3137 }
3138
3139
3140
3141 /*
3142 ;       normalize
3143 ;
3144 ; Shift normalizes the significand area pointed to by argument
3145 ; shift count (up = positive) is returned.
3146 */
3147 int
3148 enormlz (x)
3149      unsigned EMUSHORT x[];
3150 {
3151   register unsigned EMUSHORT *p;
3152   int sc;
3153
3154   sc = 0;
3155   p = &x[M];
3156   if (*p != 0)
3157     goto normdn;
3158   ++p;
3159   if (*p & 0x8000)
3160     return (0);                 /* already normalized */
3161   while (*p == 0)
3162     {
3163       eshup6 (x);
3164       sc += 16;
3165       /* With guard word, there are NBITS+16 bits available.
3166        * return true if all are zero.
3167        */
3168       if (sc > NBITS)
3169         return (sc);
3170     }
3171   /* see if high byte is zero */
3172   while ((*p & 0xff00) == 0)
3173     {
3174       eshup8 (x);
3175       sc += 8;
3176     }
3177   /* now shift 1 bit at a time */
3178   while ((*p & 0x8000) == 0)
3179     {
3180       eshup1 (x);
3181       sc += 1;
3182       if (sc > NBITS)
3183         {
3184           mtherr ("enormlz", UNDERFLOW);
3185           return (sc);
3186         }
3187     }
3188   return (sc);
3189
3190   /* Normalize by shifting down out of the high guard word
3191      of the significand */
3192  normdn:
3193
3194   if (*p & 0xff00)
3195     {
3196       eshdn8 (x);
3197       sc -= 8;
3198     }
3199   while (*p != 0)
3200     {
3201       eshdn1 (x);
3202       sc -= 1;
3203
3204       if (sc < -NBITS)
3205         {
3206           mtherr ("enormlz", OVERFLOW);
3207           return (sc);
3208         }
3209     }
3210   return (sc);
3211 }
3212
3213
3214
3215
3216 /* Convert e type number to decimal format ASCII string.
3217  * The constants are for 64 bit precision.
3218  */
3219
3220 #define NTEN 12
3221 #define MAXP 4096
3222
3223 static unsigned EMUSHORT etens[NTEN + 1][NE] =
3224 {
3225   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
3226   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
3227   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3228   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3229   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3230   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3231   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3232   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3233   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3234   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3235   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3236   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3237   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
3238 };
3239
3240 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
3241 {
3242   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
3243   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
3244   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3245   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3246   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3247   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3248   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3249   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3250   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3251   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3252   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3253   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3254   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
3255 };
3256
3257 void
3258 e24toasc (x, string, ndigs)
3259      unsigned EMUSHORT x[];
3260      char *string;
3261      int ndigs;
3262 {
3263   unsigned EMUSHORT w[NI];
3264
3265 #ifdef INFINITY
3266 #ifdef IBMPC
3267   if ((x[1] & 0x7f80) == 0x7f80)
3268 #else
3269   if ((x[0] & 0x7f80) == 0x7f80)
3270 #endif
3271     {
3272 #ifdef IBMPC
3273       if (x[1] & 0x8000)
3274 #else
3275       if (x[0] & 0x8000)
3276 #endif
3277         sprintf (string, " -Infinity ");
3278       else
3279         sprintf (string, " Infinity ");
3280       return;
3281     }
3282 #endif
3283   e24toe (x, w);
3284   etoasc (w, string, ndigs);
3285 }
3286
3287
3288 void
3289 e53toasc (x, string, ndigs)
3290      unsigned EMUSHORT x[];
3291      char *string;
3292      int ndigs;
3293 {
3294   unsigned EMUSHORT w[NI];
3295
3296 #ifdef INFINITY
3297 #ifdef IBMPC
3298   if ((x[3] & 0x7ff0) == 0x7ff0)
3299 #else
3300   if ((x[0] & 0x7ff0) == 0x7ff0)
3301 #endif
3302     {
3303 #ifdef IBMPC
3304       if (x[3] & 0x8000)
3305 #else
3306       if (x[0] & 0x8000)
3307 #endif
3308         sprintf (string, " -Infinity ");
3309       else
3310         sprintf (string, " Infinity ");
3311       return;
3312     }
3313 #endif
3314   e53toe (x, w);
3315   etoasc (w, string, ndigs);
3316 }
3317
3318
3319 void
3320 e64toasc (x, string, ndigs)
3321      unsigned EMUSHORT x[];
3322      char *string;
3323      int ndigs;
3324 {
3325   unsigned EMUSHORT w[NI];
3326
3327 #ifdef INFINITY
3328 #ifdef IBMPC
3329   if ((x[4] & 0x7fff) == 0x7fff)
3330 #else
3331   if ((x[0] & 0x7fff) == 0x7fff)
3332 #endif
3333     {
3334 #ifdef IBMPC
3335       if (x[4] & 0x8000)
3336 #else
3337       if (x[0] & 0x8000)
3338 #endif
3339         sprintf (string, " -Infinity ");
3340       else
3341         sprintf (string, " Infinity ");
3342       return;
3343     }
3344 #endif
3345   e64toe (x, w);
3346   etoasc (w, string, ndigs);
3347 }
3348
3349
3350 static char wstring[80];        /* working storage for ASCII output */
3351
3352 void
3353 etoasc (x, string, ndigs)
3354      unsigned EMUSHORT x[];
3355      char *string;
3356      int ndigs;
3357 {
3358   EMUSHORT digit;
3359   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
3360   unsigned EMUSHORT *p, *r, *ten;
3361   unsigned EMUSHORT sign;
3362   int i, j, k, expon, rndsav;
3363   char *s, *ss;
3364   unsigned EMUSHORT m;
3365
3366   ss = string;
3367   s = wstring;
3368   while ((*s++ = *ss++) != '\0')
3369     ;
3370   rndsav = rndprc;
3371   rndprc = NBITS;               /* set to full precision */
3372   emov (x, y);                  /* retain external format */
3373   if (y[NE - 1] & 0x8000)
3374     {
3375       sign = 0xffff;
3376       y[NE - 1] &= 0x7fff;
3377     }
3378   else
3379     {
3380       sign = 0;
3381     }
3382   expon = 0;
3383   ten = &etens[NTEN][0];
3384   emov (eone, t);
3385   /* Test for zero exponent */
3386   if (y[NE - 1] == 0)
3387     {
3388       for (k = 0; k < NE - 1; k++)
3389         {
3390           if (y[k] != 0)
3391             goto tnzro;         /* denormalized number */
3392         }
3393       goto isone;               /* legal all zeros */
3394     }
3395  tnzro:
3396
3397   /* Test for infinity.  Don't bother with illegal infinities.
3398    */
3399   if (y[NE - 1] == 0x7fff)
3400     {
3401       if (sign)
3402         sprintf (wstring, " -Infinity ");
3403       else
3404         sprintf (wstring, " Infinity ");
3405       goto bxit;
3406     }
3407
3408   /* Test for exponent nonzero but significand denormalized.
3409    * This is an error condition.
3410    */
3411   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3412     {
3413       mtherr ("etoasc", DOMAIN);
3414       sprintf (wstring, "NaN");
3415       goto bxit;
3416     }
3417
3418   /* Compare to 1.0 */
3419   i = ecmp (eone, y);
3420   if (i == 0)
3421     goto isone;
3422
3423   if (i < 0)
3424     {                           /* Number is greater than 1 */
3425       /* Convert significand to an integer and strip trailing decimal zeros. */
3426       emov (y, u);
3427       u[NE - 1] = EXONE + NBITS - 1;
3428
3429       p = &etens[NTEN - 4][0];
3430       m = 16;
3431       do
3432         {
3433           ediv (p, u, t);
3434           efloor (t, w);
3435           for (j = 0; j < NE - 1; j++)
3436             {
3437               if (t[j] != w[j])
3438                 goto noint;
3439             }
3440           emov (t, u);
3441           expon += (int) m;
3442         noint:
3443           p += NE;
3444           m >>= 1;
3445         }
3446       while (m != 0);
3447
3448       /* Rescale from integer significand */
3449       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3450       emov (u, y);
3451       /* Find power of 10 */
3452       emov (eone, t);
3453       m = MAXP;
3454       p = &etens[0][0];
3455       while (ecmp (ten, u) <= 0)
3456         {
3457           if (ecmp (p, u) <= 0)
3458             {
3459               ediv (p, u, u);
3460               emul (p, t, t);
3461               expon += (int) m;
3462             }
3463           m >>= 1;
3464           if (m == 0)
3465             break;
3466           p += NE;
3467         }
3468     }
3469   else
3470     {                           /* Number is less than 1.0 */
3471       /* Pad significand with trailing decimal zeros. */
3472       if (y[NE - 1] == 0)
3473         {
3474           while ((y[NE - 2] & 0x8000) == 0)
3475             {
3476               emul (ten, y, y);
3477               expon -= 1;
3478             }
3479         }
3480       else
3481         {
3482           emovi (y, w);
3483           for (i = 0; i < NDEC + 1; i++)
3484             {
3485               if ((w[NI - 1] & 0x7) != 0)
3486                 break;
3487               /* multiply by 10 */
3488               emovz (w, u);
3489               eshdn1 (u);
3490               eshdn1 (u);
3491               eaddm (w, u);
3492               u[1] += 3;
3493               while (u[2] != 0)
3494                 {
3495                   eshdn1 (u);
3496                   u[1] += 1;
3497                 }
3498               if (u[NI - 1] != 0)
3499                 break;
3500               if (eone[NE - 1] <= u[1])
3501                 break;
3502               emovz (u, w);
3503               expon -= 1;
3504             }
3505           emovo (w, y);
3506         }
3507       k = -MAXP;
3508       p = &emtens[0][0];
3509       r = &etens[0][0];
3510       emov (y, w);
3511       emov (eone, t);
3512       while (ecmp (eone, w) > 0)
3513         {
3514           if (ecmp (p, w) >= 0)
3515             {
3516               emul (r, w, w);
3517               emul (r, t, t);
3518               expon += k;
3519             }
3520           k /= 2;
3521           if (k == 0)
3522             break;
3523           p += NE;
3524           r += NE;
3525         }
3526       ediv (t, eone, t);
3527     }
3528  isone:
3529   /* Find the first (leading) digit. */
3530   emovi (t, w);
3531   emovz (w, t);
3532   emovi (y, w);
3533   emovz (w, y);
3534   eiremain (t, y);
3535   digit = equot[NI - 1];
3536   while ((digit == 0) && (ecmp (y, ezero) != 0))
3537     {
3538       eshup1 (y);
3539       emovz (y, u);
3540       eshup1 (u);
3541       eshup1 (u);
3542       eaddm (u, y);
3543       eiremain (t, y);
3544       digit = equot[NI - 1];
3545       expon -= 1;
3546     }
3547   s = wstring;
3548   if (sign)
3549     *s++ = '-';
3550   else
3551     *s++ = ' ';
3552   /* Examine number of digits requested by caller. */
3553   if (ndigs < 0)
3554     ndigs = 0;
3555   if (ndigs > NDEC)
3556     ndigs = NDEC;
3557   if (digit == 10)
3558     {
3559       *s++ = '1';
3560       *s++ = '.';
3561       if (ndigs > 0)
3562         {
3563           *s++ = '0';
3564           ndigs -= 1;
3565         }
3566       expon += 1;
3567     }
3568   else
3569     {
3570       *s++ = (char )digit + '0';
3571       *s++ = '.';
3572     }
3573   /* Generate digits after the decimal point. */
3574   for (k = 0; k <= ndigs; k++)
3575     {
3576       /* multiply current number by 10, without normalizing */
3577       eshup1 (y);
3578       emovz (y, u);
3579       eshup1 (u);
3580       eshup1 (u);
3581       eaddm (u, y);
3582       eiremain (t, y);
3583       *s++ = (char) equot[NI - 1] + '0';
3584     }
3585   digit = equot[NI - 1];
3586   --s;
3587   ss = s;
3588   /* round off the ASCII string */
3589   if (digit > 4)
3590     {
3591       /* Test for critical rounding case in ASCII output. */
3592       if (digit == 5)
3593         {
3594           emovo (y, t);
3595           if (ecmp (t, ezero) != 0)
3596             goto roun;          /* round to nearest */
3597           if ((*(s - 1) & 1) == 0)
3598             goto doexp;         /* round to even */
3599         }
3600       /* Round up and propagate carry-outs */
3601     roun:
3602       --s;
3603       k = *s & 0x7f;
3604       /* Carry out to most significant digit? */
3605       if (k == '.')
3606         {
3607           --s;
3608           k = *s;
3609           k += 1;
3610           *s = (char) k;
3611           /* Most significant digit carries to 10? */
3612           if (k > '9')
3613             {
3614               expon += 1;
3615               *s = '1';
3616             }
3617           goto doexp;
3618         }
3619       /* Round up and carry out from less significant digits */
3620       k += 1;
3621       *s = (char) k;
3622       if (k > '9')
3623         {
3624           *s = '0';
3625           goto roun;
3626         }
3627     }
3628  doexp:
3629   /*
3630      if (expon >= 0)
3631      sprintf (ss, "e+%d", expon);
3632      else
3633      sprintf (ss, "e%d", expon);
3634      */
3635   sprintf (ss, "e%d", expon);
3636  bxit:
3637   rndprc = rndsav;
3638   /* copy out the working string */
3639   s = string;
3640   ss = wstring;
3641   while (*ss == ' ')            /* strip possible leading space */
3642     ++ss;
3643   while ((*s++ = *ss++) != '\0')
3644     ;
3645 }
3646
3647
3648
3649
3650 /*
3651 ;                                                               ASCTOQ
3652 ;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
3653 ;                                       SLM, 3 JAN 78
3654 ;
3655 ;       Convert ASCII string to quadruple precision floating point
3656 ;
3657 ;               Numeric input is free field decimal number
3658 ;               with max of 15 digits with or without
3659 ;               decimal point entered as ASCII from teletype.
3660 ;       Entering E after the number followed by a second
3661 ;       number causes the second number to be interpreted
3662 ;       as a power of 10 to be multiplied by the first number
3663 ;       (i.e., "scientific" notation).
3664 ;
3665 ;       Usage:
3666 ;               asctoq (string, q);
3667 */
3668
3669 /* ASCII to single */
3670 void
3671 asctoe24 (s, y)
3672      char *s;
3673      unsigned EMUSHORT *y;
3674 {
3675   asctoeg (s, y, 24);
3676 }
3677
3678
3679 /* ASCII to double */
3680 void
3681 asctoe53 (s, y)
3682      char *s;
3683      unsigned EMUSHORT *y;
3684 {
3685 #ifdef DEC
3686   asctoeg (s, y, 56);
3687 #else
3688   asctoeg (s, y, 53);
3689 #endif
3690 }
3691
3692
3693 /* ASCII to long double */
3694 void
3695 asctoe64 (s, y)
3696      char *s;
3697      unsigned EMUSHORT *y;
3698 {
3699   asctoeg (s, y, 64);
3700 }
3701
3702 /* ASCII to super double */
3703 void
3704 asctoe (s, y)
3705      char *s;
3706      unsigned EMUSHORT *y;
3707 {
3708   asctoeg (s, y, NBITS);
3709 }
3710
3711 /* Space to make a copy of the input string: */
3712 static char lstr[82];
3713
3714 void
3715 asctoeg (ss, y, oprec)
3716      char *ss;
3717      unsigned EMUSHORT *y;
3718      int oprec;
3719 {
3720   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
3721   int esign, decflg, sgnflg, nexp, exp, prec, lost;
3722   int k, trail, c, rndsav;
3723   EMULONG lexp;
3724   unsigned EMUSHORT nsign, *p;
3725   char *sp, *s;
3726
3727   /* Copy the input string. */
3728   s = ss;
3729   while (*s == ' ')             /* skip leading spaces */
3730     ++s;
3731   sp = lstr;
3732   for (k = 0; k < 79; k++)
3733     {
3734       if ((*sp++ = *s++) == '\0')
3735         break;
3736     }
3737   *sp = '\0';
3738   s = lstr;
3739
3740   rndsav = rndprc;
3741   rndprc = NBITS;               /* Set to full precision */
3742   lost = 0;
3743   nsign = 0;
3744   decflg = 0;
3745   sgnflg = 0;
3746   nexp = 0;
3747   exp = 0;
3748   prec = 0;
3749   ecleaz (yy);
3750   trail = 0;
3751
3752  nxtcom:
3753   k = *s - '0';
3754   if ((k >= 0) && (k <= 9))
3755     {
3756       /* Ignore leading zeros */
3757       if ((prec == 0) && (decflg == 0) && (k == 0))
3758         goto donchr;
3759       /* Identify and strip trailing zeros after the decimal point. */
3760       if ((trail == 0) && (decflg != 0))
3761         {
3762           sp = s;
3763           while ((*sp >= '0') && (*sp <= '9'))
3764             ++sp;
3765           /* Check for syntax error */
3766           c = *sp & 0x7f;
3767           if ((c != 'e') && (c != 'E') && (c != '\0')
3768               && (c != '\n') && (c != '\r') && (c != ' ')
3769               && (c != ','))
3770             goto error;
3771           --sp;
3772           while (*sp == '0')
3773             *sp-- = 'z';
3774           trail = 1;
3775           if (*s == 'z')
3776             goto donchr;
3777         }
3778       /* If enough digits were given to more than fill up the yy register,
3779        * continuing until overflow into the high guard word yy[2]
3780        * guarantees that there will be a roundoff bit at the top
3781        * of the low guard word after normalization.
3782        */
3783       if (yy[2] == 0)
3784         {
3785           if (decflg)
3786             nexp += 1;          /* count digits after decimal point */
3787           eshup1 (yy);          /* multiply current number by 10 */
3788           emovz (yy, xt);
3789           eshup1 (xt);
3790           eshup1 (xt);
3791           eaddm (xt, yy);
3792           ecleaz (xt);
3793           xt[NI - 2] = (unsigned EMUSHORT) k;
3794           eaddm (xt, yy);
3795         }
3796       else
3797         {
3798           lost |= k;
3799         }
3800       prec += 1;
3801       goto donchr;
3802     }
3803
3804   switch (*s)
3805     {
3806     case 'z':
3807       break;
3808     case 'E':
3809     case 'e':
3810       goto expnt;
3811     case '.':                   /* decimal point */
3812       if (decflg)
3813         goto error;
3814       ++decflg;
3815       break;
3816     case '-':
3817       nsign = 0xffff;
3818       if (sgnflg)
3819         goto error;
3820       ++sgnflg;
3821       break;
3822     case '+':
3823       if (sgnflg)
3824         goto error;
3825       ++sgnflg;
3826       break;
3827     case ',':
3828     case ' ':
3829     case '\0':
3830     case '\n':
3831     case '\r':
3832       goto daldone;
3833     case 'i':
3834     case 'I':
3835       goto infinite;
3836     default:
3837     error:
3838       mtherr ("asctoe", DOMAIN);
3839       eclear (y);
3840       goto aexit;
3841     }
3842  donchr:
3843   ++s;
3844   goto nxtcom;
3845
3846   /* Exponent interpretation */
3847  expnt:
3848
3849   esign = 1;
3850   exp = 0;
3851   ++s;
3852   /* check for + or - */
3853   if (*s == '-')
3854     {
3855       esign = -1;
3856       ++s;
3857     }
3858   if (*s == '+')
3859     ++s;
3860   while ((*s >= '0') && (*s <= '9'))
3861     {
3862       exp *= 10;
3863       exp += *s++ - '0';
3864       if (exp > 4956)
3865         {
3866           if (esign < 0)
3867             goto zero;
3868           else
3869             goto infinite;
3870         }
3871     }
3872   if (esign < 0)
3873     exp = -exp;
3874   if (exp > 4932)
3875     {
3876  infinite:
3877       ecleaz (yy);
3878       yy[E] = 0x7fff;           /* infinity */
3879       goto aexit;
3880     }
3881   if (exp < -4956)
3882     {
3883  zero:
3884       ecleaz (yy);
3885       goto aexit;
3886     }
3887
3888  daldone:
3889   nexp = exp - nexp;
3890   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3891   while ((nexp > 0) && (yy[2] == 0))
3892     {
3893       emovz (yy, xt);
3894       eshup1 (xt);
3895       eshup1 (xt);
3896       eaddm (yy, xt);
3897       eshup1 (xt);
3898       if (xt[2] != 0)
3899         break;
3900       nexp -= 1;
3901       emovz (xt, yy);
3902     }
3903   if ((k = enormlz (yy)) > NBITS)
3904     {
3905       ecleaz (yy);
3906       goto aexit;
3907     }
3908   lexp = (EXONE - 1 + NBITS) - k;
3909   emdnorm (yy, lost, 0, lexp, 64);
3910   /* convert to external format */
3911
3912
3913   /* Multiply by 10**nexp.  If precision is 64 bits,
3914    * the maximum relative error incurred in forming 10**n
3915    * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3916    * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3917    * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3918    */
3919   lexp = yy[E];
3920   if (nexp == 0)
3921     {
3922       k = 0;
3923       goto expdon;
3924     }
3925   esign = 1;
3926   if (nexp < 0)
3927     {
3928       nexp = -nexp;
3929       esign = -1;
3930       if (nexp > 4096)
3931         {                       /* Punt.  Can't handle this without 2 divides. */
3932           emovi (etens[0], tt);
3933           lexp -= tt[E];
3934           k = edivm (tt, yy);
3935           lexp += EXONE;
3936           nexp -= 4096;
3937         }
3938     }
3939   p = &etens[NTEN][0];
3940   emov (eone, xt);
3941   exp = 1;
3942   do
3943     {
3944       if (exp & nexp)
3945         emul (p, xt, xt);
3946       p -= NE;
3947       exp = exp + exp;
3948     }
3949   while (exp <= MAXP);
3950
3951   emovi (xt, tt);
3952   if (esign < 0)
3953     {
3954       lexp -= tt[E];
3955       k = edivm (tt, yy);
3956       lexp += EXONE;
3957     }
3958   else
3959     {
3960       lexp += tt[E];
3961       k = emulm (tt, yy);
3962       lexp -= EXONE - 1;
3963     }
3964
3965  expdon:
3966
3967   /* Round and convert directly to the destination type */
3968   if (oprec == 53)
3969     lexp -= EXONE - 0x3ff;
3970   else if (oprec == 24)
3971     lexp -= EXONE - 0177;
3972 #ifdef DEC
3973   else if (oprec == 56)
3974     lexp -= EXONE - 0201;
3975 #endif
3976   rndprc = oprec;
3977   emdnorm (yy, k, 0, lexp, 64);
3978
3979  aexit:
3980
3981   rndprc = rndsav;
3982   yy[0] = nsign;
3983   switch (oprec)
3984     {
3985 #ifdef DEC
3986     case 56:
3987       todec (yy, y);            /* see etodec.c */
3988       break;
3989 #endif
3990     case 53:
3991       toe53 (yy, y);
3992       break;
3993     case 24:
3994       toe24 (yy, y);
3995       break;
3996     case 64:
3997       toe64 (yy, y);
3998       break;
3999     case NBITS:
4000       emovo (yy, y);
4001       break;
4002     }
4003 }
4004
4005
4006
4007 /* y = largest integer not greater than x
4008  * (truncated toward minus infinity)
4009  *
4010  * unsigned EMUSHORT x[NE], y[NE]
4011  *
4012  * efloor (x, y);
4013  */
4014 static unsigned EMUSHORT bmask[] =
4015 {
4016   0xffff,
4017   0xfffe,
4018   0xfffc,
4019   0xfff8,
4020   0xfff0,
4021   0xffe0,
4022   0xffc0,
4023   0xff80,
4024   0xff00,
4025   0xfe00,
4026   0xfc00,
4027   0xf800,
4028   0xf000,
4029   0xe000,
4030   0xc000,
4031   0x8000,
4032   0x0000,
4033 };
4034
4035 void
4036 efloor (x, y)
4037      unsigned EMUSHORT x[], y[];
4038 {
4039   register unsigned EMUSHORT *p;
4040   int e, expon, i;
4041   unsigned EMUSHORT f[NE];
4042
4043   emov (x, f);                  /* leave in external format */
4044   expon = (int) f[NE - 1];
4045   e = (expon & 0x7fff) - (EXONE - 1);
4046   if (e <= 0)
4047     {
4048       eclear (y);
4049       goto isitneg;
4050     }
4051   /* number of bits to clear out */
4052   e = NBITS - e;
4053   emov (f, y);
4054   if (e <= 0)
4055     return;
4056
4057   p = &y[0];
4058   while (e >= 16)
4059     {
4060       *p++ = 0;
4061       e -= 16;
4062     }
4063   /* clear the remaining bits */
4064   *p &= bmask[e];
4065   /* truncate negatives toward minus infinity */
4066  isitneg:
4067
4068   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4069     {
4070       for (i = 0; i < NE - 1; i++)
4071         {
4072           if (f[i] != y[i])
4073             {
4074               esub (eone, y, y);
4075               break;
4076             }
4077         }
4078     }
4079 }
4080
4081
4082 /* unsigned EMUSHORT x[], s[];
4083  * int *exp;
4084  *
4085  * efrexp (x, exp, s);
4086  *
4087  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
4088  * For example, 1.1 = 0.55 * 2**1
4089  * Handles denormalized numbers properly using long integer exp.
4090  */
4091 void
4092 efrexp (x, exp, s)
4093      unsigned EMUSHORT x[];
4094      int *exp;
4095      unsigned EMUSHORT s[];
4096 {
4097   unsigned EMUSHORT xi[NI];
4098   EMULONG li;
4099
4100   emovi (x, xi);
4101   li = (EMULONG) ((EMUSHORT) xi[1]);
4102
4103   if (li == 0)
4104     {
4105       li -= enormlz (xi);
4106     }
4107   xi[1] = 0x3ffe;
4108   emovo (xi, s);
4109   *exp = (int) (li - 0x3ffe);
4110 }
4111
4112
4113
4114 /* unsigned EMUSHORT x[], y[];
4115  * long pwr2;
4116  *
4117  * eldexp (x, pwr2, y);
4118  *
4119  * Returns y = x * 2**pwr2.
4120  */
4121 void
4122 eldexp (x, pwr2, y)
4123      unsigned EMUSHORT x[];
4124      int pwr2;
4125      unsigned EMUSHORT y[];
4126 {
4127   unsigned EMUSHORT xi[NI];
4128   EMULONG li;
4129   int i;
4130
4131   emovi (x, xi);
4132   li = xi[1];
4133   li += pwr2;
4134   i = 0;
4135   emdnorm (xi, i, i, li, 64);
4136   emovo (xi, y);
4137 }
4138
4139
4140 /* c = remainder after dividing b by a
4141  * Least significant integer quotient bits left in equot[].
4142  */
4143 void
4144 eremain (a, b, c)
4145      unsigned EMUSHORT a[], b[], c[];
4146 {
4147   unsigned EMUSHORT den[NI], num[NI];
4148
4149   if (ecmp (a, ezero) == 0)
4150     {
4151       mtherr ("eremain", SING);
4152       eclear (c);
4153       return;
4154     }
4155   emovi (a, den);
4156   emovi (b, num);
4157   eiremain (den, num);
4158   /* Sign of remainder = sign of quotient */
4159   if (a[0] == b[0])
4160     num[0] = 0;
4161   else
4162     num[0] = 0xffff;
4163   emovo (num, c);
4164 }
4165
4166 void
4167 eiremain (den, num)
4168      unsigned EMUSHORT den[], num[];
4169 {
4170   EMULONG ld, ln;
4171   unsigned EMUSHORT j;
4172
4173   ld = den[E];
4174   ld -= enormlz (den);
4175   ln = num[E];
4176   ln -= enormlz (num);
4177   ecleaz (equot);
4178   while (ln >= ld)
4179     {
4180       if (ecmpm (den, num) <= 0)
4181         {
4182           esubm (den, num);
4183           j = 1;
4184         }
4185       else
4186         {
4187           j = 0;
4188         }
4189       eshup1 (equot);
4190       equot[NI - 1] |= j;
4191       eshup1 (num);
4192       ln -= 1;
4193     }
4194   emdnorm (num, 0, 0, ln, 0);
4195 }
4196
4197 /*                                                      mtherr.c
4198  *
4199  *      Library common error handling routine
4200  *
4201  *
4202  *
4203  * SYNOPSIS:
4204  *
4205  * char *fctnam;
4206  * int code;
4207  * void mtherr ();
4208  *
4209  * mtherr (fctnam, code);
4210  *
4211  *
4212  *
4213  * DESCRIPTION:
4214  *
4215  * This routine may be called to report one of the following
4216  * error conditions (in the include file mconf.h).
4217  *
4218  *   Mnemonic        Value          Significance
4219  *
4220  *    DOMAIN            1       argument domain error
4221  *    SING              2       function singularity
4222  *    OVERFLOW          3       overflow range error
4223  *    UNDERFLOW         4       underflow range error
4224  *    TLOSS             5       total loss of precision
4225  *    PLOSS             6       partial loss of precision
4226  *    EDOM             33       Unix domain error code
4227  *    ERANGE           34       Unix range error code
4228  *
4229  * The default version of the file prints the function name,
4230  * passed to it by the pointer fctnam, followed by the
4231  * error condition.  The display is directed to the standard
4232  * output device.  The routine then returns to the calling
4233  * program.  Users may wish to modify the program to abort by
4234  * calling exit under severe error conditions such as domain
4235  * errors.
4236  *
4237  * Since all error conditions pass control to this function,
4238  * the display may be easily changed, eliminated, or directed
4239  * to an error logging device.
4240  *
4241  * SEE ALSO:
4242  *
4243  * mconf.h
4244  *
4245  */
4246 \f
4247 /*
4248 Cephes Math Library Release 2.0:  April, 1987
4249 Copyright 1984, 1987 by Stephen L. Moshier
4250 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4251 */
4252
4253 /* include "mconf.h" */
4254
4255 /* Notice: the order of appearance of the following
4256  * messages is bound to the error codes defined
4257  * in mconf.h.
4258  */
4259 static char *ermsg[7] =
4260 {
4261   "unknown",                    /* error code 0 */
4262   "domain",                     /* error code 1 */
4263   "singularity",                /* et seq.      */
4264   "overflow",
4265   "underflow",
4266   "total loss of precision",
4267   "partial loss of precision"
4268 };
4269
4270 int merror = 0;
4271 extern int merror;
4272
4273 void
4274 mtherr (name, code)
4275      char *name;
4276      int code;
4277 {
4278   char errstr[80];
4279
4280   /* Display string passed by calling program,
4281    * which is supposed to be the name of the
4282    * function in which the error occurred.
4283    */
4284
4285   /* Display error message defined
4286    * by the code argument.
4287    */
4288   if ((code <= 0) || (code >= 6))
4289     code = 0;
4290   sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
4291   if (extra_warnings)
4292     warning (errstr);
4293   /* Set global error message word */
4294   merror = code + 1;
4295
4296   /* Return to calling
4297    * program
4298    */
4299 }
4300
4301 /* Here is etodec.c .
4302  *
4303  */
4304
4305 /*
4306 ;       convert DEC double precision to e type
4307 ;       double d;
4308 ;       EMUSHORT e[NE];
4309 ;       dectoe (&d, e);
4310 */
4311 void
4312 dectoe (d, e)
4313      unsigned EMUSHORT *d;
4314      unsigned EMUSHORT *e;
4315 {
4316   unsigned EMUSHORT y[NI];
4317   register unsigned EMUSHORT r, *p;
4318
4319   ecleaz (y);                   /* start with a zero */
4320   p = y;                        /* point to our number */
4321   r = *d;                       /* get DEC exponent word */
4322   if (*d & (unsigned int) 0x8000)
4323     *p = 0xffff;                /* fill in our sign */
4324   ++p;                          /* bump pointer to our exponent word */
4325   r &= 0x7fff;                  /* strip the sign bit */
4326   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
4327     goto done;
4328
4329
4330   r >>= 7;                      /* shift exponent word down 7 bits */
4331   r += EXONE - 0201;            /* subtract DEC exponent offset */
4332   /* add our e type exponent offset */
4333   *p++ = r;                     /* to form our exponent */
4334
4335   r = *d++;                     /* now do the high order mantissa */
4336   r &= 0177;                    /* strip off the DEC exponent and sign bits */
4337   r |= 0200;                    /* the DEC understood high order mantissa bit */
4338   *p++ = r;                     /* put result in our high guard word */
4339
4340   *p++ = *d++;                  /* fill in the rest of our mantissa */
4341   *p++ = *d++;
4342   *p = *d;
4343
4344   eshdn8 (y);                   /* shift our mantissa down 8 bits */
4345  done:
4346   emovo (y, e);
4347 }
4348
4349
4350
4351 /*
4352 ;       convert e type to DEC double precision
4353 ;       double d;
4354 ;       EMUSHORT e[NE];
4355 ;       etodec (e, &d);
4356 */
4357 #if 0
4358 static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4359
4360 void
4361 etodec (x, d)
4362      unsigned EMUSHORT *x, *d;
4363 {
4364   unsigned EMUSHORT xi[NI];
4365   register unsigned EMUSHORT r;
4366   int i, j;
4367
4368   emovi (x, xi);
4369   *d = 0;
4370   if (xi[0] != 0)
4371     *d = 0100000;
4372   r = xi[E];
4373   if (r < (EXONE - 128))
4374     goto zout;
4375   i = xi[M + 4];
4376   if ((i & 0200) != 0)
4377     {
4378       if ((i & 0377) == 0200)
4379         {
4380           if ((i & 0400) != 0)
4381             {
4382               /* check all less significant bits */
4383               for (j = M + 5; j < NI; j++)
4384                 {
4385                   if (xi[j] != 0)
4386                     goto yesrnd;
4387                 }
4388             }
4389           goto nornd;
4390         }
4391     yesrnd:
4392       eaddm (decbit, xi);
4393       r -= enormlz (xi);
4394     }
4395
4396  nornd:
4397
4398   r -= EXONE;
4399   r += 0201;
4400   if (r < 0)
4401     {
4402     zout:
4403       *d++ = 0;
4404       *d++ = 0;
4405       *d++ = 0;
4406       *d++ = 0;
4407       return;
4408     }
4409   if (r >= 0377)
4410     {
4411       *d++ = 077777;
4412       *d++ = -1;
4413       *d++ = -1;
4414       *d++ = -1;
4415       return;
4416     }
4417   r &= 0377;
4418   r <<= 7;
4419   eshup8 (xi);
4420   xi[M] &= 0177;
4421   r |= xi[M];
4422   *d++ |= r;
4423   *d++ = xi[M + 1];
4424   *d++ = xi[M + 2];
4425   *d++ = xi[M + 3];
4426 }
4427
4428 #else
4429
4430 void
4431 etodec (x, d)
4432      unsigned EMUSHORT *x, *d;
4433 {
4434   unsigned EMUSHORT xi[NI];
4435   EMULONG exp;
4436   int rndsav;
4437
4438   emovi (x, xi);
4439   exp = (EMULONG) xi[E] - (EXONE - 0201);       /* adjust exponent for offsets */
4440 /* round off to nearest or even */
4441   rndsav = rndprc;
4442   rndprc = 56;
4443   emdnorm (xi, 0, 0, exp, 64);
4444   rndprc = rndsav;
4445   todec (xi, d);
4446 }
4447
4448 void
4449 todec (x, y)
4450      unsigned EMUSHORT *x, *y;
4451 {
4452   unsigned EMUSHORT i;
4453   unsigned EMUSHORT *p;
4454
4455   p = x;
4456   *y = 0;
4457   if (*p++)
4458     *y = 0100000;
4459   i = *p++;
4460   if (i == 0)
4461     {
4462       *y++ = 0;
4463       *y++ = 0;
4464       *y++ = 0;
4465       *y++ = 0;
4466       return;
4467     }
4468   if (i > 0377)
4469     {
4470       *y++ |= 077777;
4471       *y++ = 0xffff;
4472       *y++ = 0xffff;
4473       *y++ = 0xffff;
4474 #ifdef ERANGE
4475       errno = ERANGE;
4476 #endif
4477       return;
4478     }
4479   i &= 0377;
4480   i <<= 7;
4481   eshup8 (x);
4482   x[M] &= 0177;
4483   i |= x[M];
4484   *y++ |= i;
4485   *y++ = x[M + 1];
4486   *y++ = x[M + 2];
4487   *y++ = x[M + 3];
4488 }
4489
4490 #endif /* not 0 */
4491
4492 #endif /* EMU_NON_COMPILE not defined */