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 #define EDOM            33
 995 #define ERANGE          34
 996
 997 /*  e type constants used by high precision check routines */
 998
 999 /*include "ehead.h"*/
1000 /* 0.0 */
1001 unsigned EMUSHORT ezero[NE] =
1002 {
1003   0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1004 extern unsigned EMUSHORT ezero[];
1005
1006 /* 5.0E-1 */
1007 unsigned EMUSHORT ehalf[NE] =
1008 {
1009   0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1010 extern unsigned EMUSHORT ehalf[];
1011
1012 /* 1.0E0 */
1013 unsigned EMUSHORT eone[NE] =
1014 {
1015   0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1016 extern unsigned EMUSHORT eone[];
1017
1018 /* 2.0E0 */
1019 unsigned EMUSHORT etwo[NE] =
1020 {
1021   0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1022 extern unsigned EMUSHORT etwo[];
1023
1024 /* 3.2E1 */
1025 unsigned EMUSHORT e32[NE] =
1026 {
1027   0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1028 extern unsigned EMUSHORT e32[];
1029
1030 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1031 unsigned EMUSHORT elog2[NE] =
1032 {
1033   0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1034 extern unsigned EMUSHORT elog2[];
1035
1036 /* 1.41421356237309504880168872420969807856967187537695E0 */
1037 unsigned EMUSHORT esqrt2[NE] =
1038 {
1039   0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1040 extern unsigned EMUSHORT esqrt2[];
1041
1042 /* 2/sqrt (PI) =
1043  * 1.12837916709551257389615890312154517168810125865800E0 */
1044 unsigned EMUSHORT eoneopi[NE] =
1045 {
1046   0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1047 extern unsigned EMUSHORT eoneopi[];
1048
1049 /* 3.14159265358979323846264338327950288419716939937511E0 */
1050 unsigned EMUSHORT epi[NE] =
1051 {
1052   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1053 extern unsigned EMUSHORT epi[];
1054
1055 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1056 unsigned EMUSHORT eeul[NE] =
1057 {
1058   0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1059 extern unsigned EMUSHORT eeul[];
1060
1061 /*
1062 include "ehead.h"
1063 include "mconf.h"
1064 */
1065
1066
1067
1068 /* Control register for rounding precision.
1069  * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1070  */
1071 int rndprc = NBITS;
1072 extern int rndprc;
1073
1074 void eaddm (), esubm (), emdnorm (), asctoeg ();
1075 static void toe24 (), toe53 (), toe64 ();
1076 void eremain (), einit (), eiremain ();
1077 int ecmpm (), edivm (), emulm ();
1078 void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1079 void etodec (), todec (), dectoe ();
1080
1081
1082
1083
1084 void
1085 einit ()
1086 {
1087 }
1088
1089 /*
1090 ; Clear out entire external format number.
1091 ;
1092 ; unsigned EMUSHORT x[];
1093 ; eclear (x);
1094 */
1095
1096 void
1097 eclear (x)
1098      register unsigned EMUSHORT *x;
1099 {
1100   register int i;
1101
1102   for (i = 0; i < NE; i++)
1103     *x++ = 0;
1104 }
1105
1106
1107
1108 /* Move external format number from a to b.
1109  *
1110  * emov (a, b);
1111  */
1112
1113 void
1114 emov (a, b)
1115      register unsigned EMUSHORT *a, *b;
1116 {
1117   register int i;
1118
1119   for (i = 0; i < NE; i++)
1120     *b++ = *a++;
1121 }
1122
1123
1124 /*
1125 ;       Absolute value of external format number
1126 ;
1127 ;       EMUSHORT x[NE];
1128 ;       eabs (x);
1129 */
1130
1131 void
1132 eabs (x)
1133      unsigned EMUSHORT x[];     /* x is the memory address of a short */
1134 {
1135
1136   x[NE - 1] &= 0x7fff;          /* sign is top bit of last word of external format */
1137 }
1138
1139
1140
1141
1142 /*
1143 ;       Negate external format number
1144 ;
1145 ;       unsigned EMUSHORT x[NE];
1146 ;       eneg (x);
1147 */
1148
1149 void
1150 eneg (x)
1151      unsigned EMUSHORT x[];
1152 {
1153
1154   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1155 }
1156
1157
1158
1159 /* Return 1 if external format number is negative,
1160  * else return zero.
1161  */
1162 int
1163 eisneg (x)
1164      unsigned EMUSHORT x[];
1165 {
1166
1167   if (x[NE - 1] & 0x8000)
1168     return (1);
1169   else
1170     return (0);
1171 }
1172
1173
1174 /* Return 1 if external format number has maximum possible exponent,
1175  * else return zero.
1176  */
1177 int
1178 eisinf (x)
1179      unsigned EMUSHORT x[];
1180 {
1181
1182   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1183     return (1);
1184   else
1185     return (0);
1186 }
1187
1188
1189 /*
1190 ; Fill entire number, including exponent and significand, with
1191 ; largest possible number.  These programs implement a saturation
1192 ; value that is an ordinary, legal number.  A special value
1193 ; "infinity" may also be implemented; this would require tests
1194 ; for that value and implementation of special rules for arithmetic
1195 ; operations involving inifinity.
1196 */
1197
1198 void
1199 einfin (x)
1200      register unsigned EMUSHORT *x;
1201 {
1202   register int i;
1203
1204 #ifdef INFINITY
1205   for (i = 0; i < NE - 1; i++)
1206     *x++ = 0;
1207   *x |= 32767;
1208 #else
1209   for (i = 0; i < NE - 1; i++)
1210     *x++ = 0xffff;
1211   *x |= 32766;
1212   if (rndprc < NBITS)
1213     {
1214       if (rndprc == 64)
1215         {
1216           *(x - 5) = 0;
1217         }
1218       if (rndprc == 53)
1219         {
1220           *(x - 4) = 0xf800;
1221         }
1222       else
1223         {
1224           *(x - 4) = 0;
1225           *(x - 3) = 0;
1226           *(x - 2) = 0xff00;
1227         }
1228     }
1229 #endif
1230 }
1231
1232
1233
1234 /* Move in external format number,
1235  * converting it to internal format.
1236  */
1237 void
1238 emovi (a, b)
1239      unsigned EMUSHORT *a, *b;
1240 {
1241   register unsigned EMUSHORT *p, *q;
1242   int i;
1243
1244   q = b;
1245   p = a + (NE - 1);             /* point to last word of external number */
1246   /* get the sign bit */
1247   if (*p & 0x8000)
1248     *q++ = 0xffff;
1249   else
1250     *q++ = 0;
1251   /* get the exponent */
1252   *q = *p--;
1253   *q++ &= 0x7fff;               /* delete the sign bit */
1254 #ifdef INFINITY
1255   if ((*(q - 1) & 0x7fff) == 0x7fff)
1256     {
1257       for (i = 2; i < NI; i++)
1258         *q++ = 0;
1259       return;
1260     }
1261 #endif
1262   /* clear high guard word */
1263   *q++ = 0;
1264   /* move in the significand */
1265   for (i = 0; i < NE - 1; i++)
1266     *q++ = *p--;
1267   /* clear low guard word */
1268   *q = 0;
1269 }
1270
1271
1272 /* Move internal format number out,
1273  * converting it to external format.
1274  */
1275 void
1276 emovo (a, b)
1277      unsigned EMUSHORT *a, *b;
1278 {
1279   register unsigned EMUSHORT *p, *q;
1280   unsigned EMUSHORT i;
1281
1282   p = a;
1283   q = b + (NE - 1);             /* point to output exponent */
1284   /* combine sign and exponent */
1285   i = *p++;
1286   if (i)
1287     *q-- = *p++ | 0x8000;
1288   else
1289     *q-- = *p++;
1290 #ifdef INFINITY
1291   if (*(p - 1) == 0x7fff)
1292     {
1293       einfin (b);
1294       return;
1295     }
1296 #endif
1297   /* skip over guard word */
1298   ++p;
1299   /* move the significand */
1300   for (i = 0; i < NE - 1; i++)
1301     *q-- = *p++;
1302 }
1303
1304
1305
1306
1307 /* Clear out internal format number.
1308  */
1309
1310 void
1311 ecleaz (xi)
1312      register unsigned EMUSHORT *xi;
1313 {
1314   register int i;
1315
1316   for (i = 0; i < NI; i++)
1317     *xi++ = 0;
1318 }
1319
1320
1321 /* same, but don't touch the sign. */
1322
1323 void
1324 ecleazs (xi)
1325      register unsigned EMUSHORT *xi;
1326 {
1327   register int i;
1328
1329   ++xi;
1330   for (i = 0; i < NI - 1; i++)
1331     *xi++ = 0;
1332 }
1333
1334
1335
1336 /* Move internal format number from a to b.
1337  */
1338 void
1339 emovz (a, b)
1340      register unsigned EMUSHORT *a, *b;
1341 {
1342   register int i;
1343
1344   for (i = 0; i < NI - 1; i++)
1345     *b++ = *a++;
1346   /* clear low guard word */
1347   *b = 0;
1348 }
1349
1350
1351 /*
1352 ;       Compare significands of numbers in internal format.
1353 ;       Guard words are included in the comparison.
1354 ;
1355 ;       unsigned EMUSHORT a[NI], b[NI];
1356 ;       cmpm (a, b);
1357 ;
1358 ;       for the significands:
1359 ;       returns +1 if a > b
1360 ;                0 if a == b
1361 ;               -1 if a < b
1362 */
1363 int
1364 ecmpm (a, b)
1365      register unsigned EMUSHORT *a, *b;
1366 {
1367   int i;
1368
1369   a += M;                       /* skip up to significand area */
1370   b += M;
1371   for (i = M; i < NI; i++)
1372     {
1373       if (*a++ != *b++)
1374         goto difrnt;
1375     }
1376   return (0);
1377
1378  difrnt:
1379   if (*(--a) > *(--b))
1380     return (1);
1381   else
1382     return (-1);
1383 }
1384
1385
1386 /*
1387 ;       Shift significand down by 1 bit
1388 */
1389
1390 void
1391 eshdn1 (x)
1392      register unsigned EMUSHORT *x;
1393 {
1394   register unsigned EMUSHORT bits;
1395   int i;
1396
1397   x += M;                       /* point to significand area */
1398
1399   bits = 0;
1400   for (i = M; i < NI; i++)
1401     {
1402       if (*x & 1)
1403         bits |= 1;
1404       *x >>= 1;
1405       if (bits & 2)
1406         *x |= 0x8000;
1407       bits <<= 1;
1408       ++x;
1409     }
1410 }
1411
1412
1413
1414 /*
1415 ;       Shift significand up by 1 bit
1416 */
1417
1418 void
1419 eshup1 (x)
1420      register unsigned EMUSHORT *x;
1421 {
1422   register unsigned EMUSHORT bits;
1423   int i;
1424
1425   x += NI - 1;
1426   bits = 0;
1427
1428   for (i = M; i < NI; i++)
1429     {
1430       if (*x & 0x8000)
1431         bits |= 1;
1432       *x <<= 1;
1433       if (bits & 2)
1434         *x |= 1;
1435       bits <<= 1;
1436       --x;
1437     }
1438 }
1439
1440
1441
1442 /*
1443 ;       Shift significand down by 8 bits
1444 */
1445
1446 void
1447 eshdn8 (x)
1448      register unsigned EMUSHORT *x;
1449 {
1450   register unsigned EMUSHORT newbyt, oldbyt;
1451   int i;
1452
1453   x += M;
1454   oldbyt = 0;
1455   for (i = M; i < NI; i++)
1456     {
1457       newbyt = *x << 8;
1458       *x >>= 8;
1459       *x |= oldbyt;
1460       oldbyt = newbyt;
1461       ++x;
1462     }
1463 }
1464
1465 /*
1466 ;       Shift significand up by 8 bits
1467 */
1468
1469 void
1470 eshup8 (x)
1471      register unsigned EMUSHORT *x;
1472 {
1473   int i;
1474   register unsigned EMUSHORT newbyt, oldbyt;
1475
1476   x += NI - 1;
1477   oldbyt = 0;
1478
1479   for (i = M; i < NI; i++)
1480     {
1481       newbyt = *x >> 8;
1482       *x <<= 8;
1483       *x |= oldbyt;
1484       oldbyt = newbyt;
1485       --x;
1486     }
1487 }
1488
1489 /*
1490 ;       Shift significand up by 16 bits
1491 */
1492
1493 void
1494 eshup6 (x)
1495      register unsigned EMUSHORT *x;
1496 {
1497   int i;
1498   register unsigned EMUSHORT *p;
1499
1500   p = x + M;
1501   x += M + 1;
1502
1503   for (i = M; i < NI - 1; i++)
1504     *p++ = *x++;
1505
1506   *p = 0;
1507 }
1508
1509 /*
1510 ;       Shift significand down by 16 bits
1511 */
1512
1513 void
1514 eshdn6 (x)
1515      register unsigned EMUSHORT *x;
1516 {
1517   int i;
1518   register unsigned EMUSHORT *p;
1519
1520   x += NI - 1;
1521   p = x + 1;
1522
1523   for (i = M; i < NI - 1; i++)
1524     *(--p) = *(--x);
1525
1526   *(--p) = 0;
1527 }
1528 \f
1529 /*
1530 ;       Add significands
1531 ;       x + y replaces y
1532 */
1533
1534 void
1535 eaddm (x, y)
1536      unsigned EMUSHORT *x, *y;
1537 {
1538   register unsigned EMULONG a;
1539   int i;
1540   unsigned int carry;
1541
1542   x += NI - 1;
1543   y += NI - 1;
1544   carry = 0;
1545   for (i = M; i < NI; i++)
1546     {
1547       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1548       if (a & 0x10000)
1549         carry = 1;
1550       else
1551         carry = 0;
1552       *y = (unsigned EMUSHORT) a;
1553       --x;
1554       --y;
1555     }
1556 }
1557
1558 /*
1559 ;       Subtract significands
1560 ;       y - x replaces y
1561 */
1562
1563 void
1564 esubm (x, y)
1565      unsigned EMUSHORT *x, *y;
1566 {
1567   unsigned EMULONG a;
1568   int i;
1569   unsigned int carry;
1570
1571   x += NI - 1;
1572   y += NI - 1;
1573   carry = 0;
1574   for (i = M; i < NI; i++)
1575     {
1576       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1577       if (a & 0x10000)
1578         carry = 1;
1579       else
1580         carry = 0;
1581       *y = (unsigned EMUSHORT) a;
1582       --x;
1583       --y;
1584     }
1585 }
1586
1587
1588 /* Divide significands */
1589
1590 static unsigned EMUSHORT equot[NI];
1591
1592 int
1593 edivm (den, num)
1594      unsigned EMUSHORT den[], num[];
1595 {
1596   int i;
1597   register unsigned EMUSHORT *p, *q;
1598   unsigned EMUSHORT j;
1599
1600   p = &equot[0];
1601   *p++ = num[0];
1602   *p++ = num[1];
1603
1604   for (i = M; i < NI; i++)
1605     {
1606       *p++ = 0;
1607     }
1608
1609   /* Use faster compare and subtraction if denominator
1610    * has only 15 bits of significance.
1611    */
1612   p = &den[M + 2];
1613   if (*p++ == 0)
1614     {
1615       for (i = M + 3; i < NI; i++)
1616         {
1617           if (*p++ != 0)
1618             goto fulldiv;
1619         }
1620       if ((den[M + 1] & 1) != 0)
1621         goto fulldiv;
1622       eshdn1 (num);
1623       eshdn1 (den);
1624
1625       p = &den[M + 1];
1626       q = &num[M + 1];
1627
1628       for (i = 0; i < NBITS + 2; i++)
1629         {
1630           if (*p <= *q)
1631             {
1632               *q -= *p;
1633               j = 1;
1634             }
1635           else
1636             {
1637               j = 0;
1638             }
1639           eshup1 (equot);
1640           equot[NI - 2] |= j;
1641           eshup1 (num);
1642         }
1643       goto divdon;
1644     }
1645
1646   /* The number of quotient bits to calculate is
1647    * NBITS + 1 scaling guard bit + 1 roundoff bit.
1648    */
1649  fulldiv:
1650
1651   p = &equot[NI - 2];
1652   for (i = 0; i < NBITS + 2; i++)
1653     {
1654       if (ecmpm (den, num) <= 0)
1655         {
1656           esubm (den, num);
1657           j = 1;                /* quotient bit = 1 */
1658         }
1659       else
1660         j = 0;
1661       eshup1 (equot);
1662       *p |= j;
1663       eshup1 (num);
1664     }
1665
1666  divdon:
1667
1668   eshdn1 (equot);
1669   eshdn1 (equot);
1670
1671   /* test for nonzero remainder after roundoff bit */
1672   p = &num[M];
1673   j = 0;
1674   for (i = M; i < NI; i++)
1675     {
1676       j |= *p++;
1677     }
1678   if (j)
1679     j = 1;
1680
1681
1682   for (i = 0; i < NI; i++)
1683     num[i] = equot[i];
1684   return ((int) j);
1685 }
1686
1687
1688 /* Multiply significands */
1689 int
1690 emulm (a, b)
1691      unsigned EMUSHORT a[], b[];
1692 {
1693   unsigned EMUSHORT *p, *q;
1694   int i, j, k;
1695
1696   equot[0] = b[0];
1697   equot[1] = b[1];
1698   for (i = M; i < NI; i++)
1699     equot[i] = 0;
1700
1701   p = &a[NI - 2];
1702   k = NBITS;
1703   while (*p == 0)               /* significand is not supposed to be all zero */
1704     {
1705       eshdn6 (a);
1706       k -= 16;
1707     }
1708   if ((*p & 0xff) == 0)
1709     {
1710       eshdn8 (a);
1711       k -= 8;
1712     }
1713
1714   q = &equot[NI - 1];
1715   j = 0;
1716   for (i = 0; i < k; i++)
1717     {
1718       if (*p & 1)
1719         eaddm (b, equot);
1720       /* remember if there were any nonzero bits shifted out */
1721       if (*q & 1)
1722         j |= 1;
1723       eshdn1 (a);
1724       eshdn1 (equot);
1725     }
1726
1727   for (i = 0; i < NI; i++)
1728     b[i] = equot[i];
1729
1730   /* return flag for lost nonzero bits */
1731   return (j);
1732 }
1733
1734
1735
1736 /*
1737  * Normalize and round off.
1738  *
1739  * The internal format number to be rounded is "s".
1740  * Input "lost" indicates whether or not the number is exact.
1741  * This is the so-called sticky bit.
1742  *
1743  * Input "subflg" indicates whether the number was obtained
1744  * by a subtraction operation.  In that case if lost is nonzero
1745  * then the number is slightly smaller than indicated.
1746  *
1747  * Input "exp" is the biased exponent, which may be negative.
1748  * the exponent field of "s" is ignored but is replaced by
1749  * "exp" as adjusted by normalization and rounding.
1750  *
1751  * Input "rcntrl" is the rounding control.
1752  */
1753
1754 static int rlast = -1;
1755 static int rw = 0;
1756 static unsigned EMUSHORT rmsk = 0;
1757 static unsigned EMUSHORT rmbit = 0;
1758 static unsigned EMUSHORT rebit = 0;
1759 static int re = 0;
1760 static unsigned EMUSHORT rbit[NI];
1761
1762 void
1763 emdnorm (s, lost, subflg, exp, rcntrl)
1764      unsigned EMUSHORT s[];
1765      int lost;
1766      int subflg;
1767      EMULONG exp;
1768      int rcntrl;
1769 {
1770   int i, j;
1771   unsigned EMUSHORT r;
1772
1773   /* Normalize */
1774   j = enormlz (s);
1775
1776   /* a blank significand could mean either zero or infinity. */
1777 #ifndef INFINITY
1778   if (j > NBITS)
1779     {
1780       ecleazs (s);
1781       return;
1782     }
1783 #endif
1784   exp -= j;
1785 #ifndef INFINITY
1786   if (exp >= 32767L)
1787     goto overf;
1788 #else
1789   if ((j > NBITS) && (exp < 32767))
1790     {
1791       ecleazs (s);
1792       return;
1793     }
1794 #endif
1795   if (exp < 0L)
1796     {
1797       if (exp > (EMULONG) (-NBITS - 1))
1798         {
1799           j = (int) exp;
1800           i = eshift (s, j);
1801           if (i)
1802             lost = 1;
1803         }
1804       else
1805         {
1806           ecleazs (s);
1807           return;
1808         }
1809     }
1810   /* Round off, unless told not to by rcntrl. */
1811   if (rcntrl == 0)
1812     goto mdfin;
1813   /* Set up rounding parameters if the control register changed. */
1814   if (rndprc != rlast)
1815     {
1816       ecleaz (rbit);
1817       switch (rndprc)
1818         {
1819         default:
1820         case NBITS:
1821           rw = NI - 1;          /* low guard word */
1822           rmsk = 0xffff;
1823           rmbit = 0x8000;
1824           rbit[rw - 1] = 1;
1825           re = NI - 2;
1826           rebit = 1;
1827           break;
1828         case 64:
1829           rw = 7;
1830           rmsk = 0xffff;
1831           rmbit = 0x8000;
1832           rbit[rw - 1] = 1;
1833           re = rw - 1;
1834           rebit = 1;
1835           break;
1836           /* For DEC arithmetic */
1837         case 56:
1838           rw = 6;
1839           rmsk = 0xff;
1840           rmbit = 0x80;
1841           rbit[rw] = 0x100;
1842           re = rw;
1843           rebit = 0x100;
1844           break;
1845         case 53:
1846           rw = 6;
1847           rmsk = 0x7ff;
1848           rmbit = 0x0400;
1849           rbit[rw] = 0x800;
1850           re = rw;
1851           rebit = 0x800;
1852           break;
1853         case 24:
1854           rw = 4;
1855           rmsk = 0xff;
1856           rmbit = 0x80;
1857           rbit[rw] = 0x100;
1858           re = rw;
1859           rebit = 0x100;
1860           break;
1861         }
1862       rlast = rndprc;
1863     }
1864
1865   if (rndprc >= 64)
1866     {
1867       r = s[rw] & rmsk;
1868       if (rndprc == 64)
1869         {
1870           i = rw + 1;
1871           while (i < NI)
1872             {
1873               if (s[i])
1874                 r |= 1;
1875               s[i] = 0;
1876               ++i;
1877             }
1878         }
1879     }
1880   else
1881     {
1882       if (exp <= 0)
1883         eshdn1 (s);
1884       r = s[rw] & rmsk;
1885       /* These tests assume NI = 8 */
1886       i = rw + 1;
1887       while (i < NI)
1888         {
1889           if (s[i])
1890             r |= 1;
1891           s[i] = 0;
1892           ++i;
1893         }
1894       /*
1895          if (rndprc == 24)
1896          {
1897          if (s[5] || s[6])
1898          r |= 1;
1899          s[5] = 0;
1900          s[6] = 0;
1901          }
1902          */
1903       s[rw] &= ~rmsk;
1904     }
1905   if ((r & rmbit) != 0)
1906     {
1907       if (r == rmbit)
1908         {
1909           if (lost == 0)
1910             {                   /* round to even */
1911               if ((s[re] & rebit) == 0)
1912                 goto mddone;
1913             }
1914           else
1915             {
1916               if (subflg != 0)
1917                 goto mddone;
1918             }
1919         }
1920       eaddm (rbit, s);
1921     }
1922  mddone:
1923   if ((rndprc < 64) && (exp <= 0))
1924     {
1925       eshup1 (s);
1926     }
1927   if (s[2] != 0)
1928     {                           /* overflow on roundoff */
1929       eshdn1 (s);
1930       exp += 1;
1931     }
1932  mdfin:
1933   s[NI - 1] = 0;
1934   if (exp >= 32767L)
1935     {
1936 #ifndef INFINITY
1937     overf:
1938 #endif
1939 #ifdef INFINITY
1940       s[1] = 32767;
1941       for (i = 2; i < NI - 1; i++)
1942         s[i] = 0;
1943       if (extra_warnings)
1944         warning ("floating point overflow");
1945 #else
1946       s[1] = 32766;
1947       s[2] = 0;
1948       for (i = M + 1; i < NI - 1; i++)
1949         s[i] = 0xffff;
1950       s[NI - 1] = 0;
1951       if (rndprc < 64)
1952         {
1953           s[rw] &= ~rmsk;
1954           if (rndprc == 24)
1955             {
1956               s[5] = 0;
1957               s[6] = 0;
1958             }
1959         }
1960 #endif
1961       return;
1962     }
1963   if (exp < 0)
1964     s[1] = 0;
1965   else
1966     s[1] = (unsigned EMUSHORT) exp;
1967 }
1968
1969
1970
1971 /*
1972 ;       Subtract external format numbers.
1973 ;
1974 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1975 ;       esub (a, b, c);  c = b - a
1976 */
1977
1978 static int subflg = 0;
1979
1980 void
1981 esub (a, b, c)
1982      unsigned EMUSHORT *a, *b, *c;
1983 {
1984
1985   subflg = 1;
1986   eadd1 (a, b, c);
1987 }
1988
1989
1990 /*
1991 ;       Add.
1992 ;
1993 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1994 ;       eadd (a, b, c);  c = b + a
1995 */
1996 void
1997 eadd (a, b, c)
1998      unsigned EMUSHORT *a, *b, *c;
1999 {
2000
2001   subflg = 0;
2002   eadd1 (a, b, c);
2003 }
2004
2005 void
2006 eadd1 (a, b, c)
2007      unsigned EMUSHORT *a, *b, *c;
2008 {
2009   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2010   int i, lost, j, k;
2011   EMULONG lt, lta, ltb;
2012
2013 #ifdef INFINITY
2014   if (eisinf (a))
2015     {
2016       emov (a, c);
2017       if (subflg)
2018         eneg (c);
2019       return;
2020     }
2021   if (eisinf (b))
2022     {
2023       emov (b, c);
2024       return;
2025     }
2026 #endif
2027   emovi (a, ai);
2028   emovi (b, bi);
2029   if (subflg)
2030     ai[0] = ~ai[0];
2031
2032   /* compare exponents */
2033   lta = ai[E];
2034   ltb = bi[E];
2035   lt = lta - ltb;
2036   if (lt > 0L)
2037     {                           /* put the larger number in bi */
2038       emovz (bi, ci);
2039       emovz (ai, bi);
2040       emovz (ci, ai);
2041       ltb = bi[E];
2042       lt = -lt;
2043     }
2044   lost = 0;
2045   if (lt != 0L)
2046     {
2047       if (lt < (EMULONG) (-NBITS - 1))
2048         goto done;              /* answer same as larger addend */
2049       k = (int) lt;
2050       lost = eshift (ai, k);    /* shift the smaller number down */
2051     }
2052   else
2053     {
2054       /* exponents were the same, so must compare significands */
2055       i = ecmpm (ai, bi);
2056       if (i == 0)
2057         {                       /* the numbers are identical in magnitude */
2058           /* if different signs, result is zero */
2059           if (ai[0] != bi[0])
2060             {
2061               eclear (c);
2062               return;
2063             }
2064           /* if same sign, result is double */
2065           /* double denomalized tiny number */
2066           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2067             {
2068               eshup1 (bi);
2069               goto done;
2070             }
2071           /* add 1 to exponent unless both are zero! */
2072           for (j = 1; j < NI - 1; j++)
2073             {
2074               if (bi[j] != 0)
2075                 {
2076                   /* This could overflow, but let emovo take care of that. */
2077                   ltb += 1;
2078                   break;
2079                 }
2080             }
2081           bi[E] = (unsigned EMUSHORT) ltb;
2082           goto done;
2083         }
2084       if (i > 0)
2085         {                       /* put the larger number in bi */
2086           emovz (bi, ci);
2087           emovz (ai, bi);
2088           emovz (ci, ai);
2089         }
2090     }
2091   if (ai[0] == bi[0])
2092     {
2093       eaddm (ai, bi);
2094       subflg = 0;
2095     }
2096   else
2097     {
2098       esubm (ai, bi);
2099       subflg = 1;
2100     }
2101   emdnorm (bi, lost, subflg, ltb, 64);
2102
2103  done:
2104   emovo (bi, c);
2105 }
2106
2107
2108
2109 /*
2110 ;       Divide.
2111 ;
2112 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2113 ;       ediv (a, b, c); c = b / a
2114 */
2115 void
2116 ediv (a, b, c)
2117      unsigned EMUSHORT *a, *b, *c;
2118 {
2119   unsigned EMUSHORT ai[NI], bi[NI];
2120   int i;
2121   EMULONG lt, lta, ltb;
2122
2123 #ifdef INFINITY
2124   if (eisinf (b))
2125     {
2126       if (eisneg (a) ^ eisneg (b))
2127         *(c + (NE - 1)) = 0x8000;
2128       else
2129         *(c + (NE - 1)) = 0;
2130       einfin (c);
2131       return;
2132     }
2133   if (eisinf (a))
2134     {
2135       eclear (c);
2136       return;
2137     }
2138 #endif
2139   emovi (a, ai);
2140   emovi (b, bi);
2141   lta = ai[E];
2142   ltb = bi[E];
2143   if (bi[E] == 0)
2144     {                           /* See if numerator is zero. */
2145       for (i = 1; i < NI - 1; i++)
2146         {
2147           if (bi[i] != 0)
2148             {
2149               ltb -= enormlz (bi);
2150               goto dnzro1;
2151             }
2152         }
2153       eclear (c);
2154       return;
2155     }
2156  dnzro1:
2157
2158   if (ai[E] == 0)
2159     {                           /* possible divide by zero */
2160       for (i = 1; i < NI - 1; i++)
2161         {
2162           if (ai[i] != 0)
2163             {
2164               lta -= enormlz (ai);
2165               goto dnzro2;
2166             }
2167         }
2168       if (ai[0] == bi[0])
2169         *(c + (NE - 1)) = 0;
2170       else
2171         *(c + (NE - 1)) = 0x8000;
2172       einfin (c);
2173       mtherr ("ediv", SING);
2174       return;
2175     }
2176  dnzro2:
2177
2178   i = edivm (ai, bi);
2179   /* calculate exponent */
2180   lt = ltb - lta + EXONE;
2181   emdnorm (bi, i, 0, lt, 64);
2182   /* set the sign */
2183   if (ai[0] == bi[0])
2184     bi[0] = 0;
2185   else
2186     bi[0] = 0Xffff;
2187   emovo (bi, c);
2188 }
2189
2190
2191
2192 /*
2193 ;       Multiply.
2194 ;
2195 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2196 ;       emul (a, b, c); c = b * a
2197 */
2198 void
2199 emul (a, b, c)
2200      unsigned EMUSHORT *a, *b, *c;
2201 {
2202   unsigned EMUSHORT ai[NI], bi[NI];
2203   int i, j;
2204   EMULONG lt, lta, ltb;
2205
2206 #ifdef INFINITY
2207   if (eisinf (a) || eisinf (b))
2208     {
2209       if (eisneg (a) ^ eisneg (b))
2210         *(c + (NE - 1)) = 0x8000;
2211       else
2212         *(c + (NE - 1)) = 0;
2213       einfin (c);
2214       return;
2215     }
2216 #endif
2217   emovi (a, ai);
2218   emovi (b, bi);
2219   lta = ai[E];
2220   ltb = bi[E];
2221   if (ai[E] == 0)
2222     {
2223       for (i = 1; i < NI - 1; i++)
2224         {
2225           if (ai[i] != 0)
2226             {
2227               lta -= enormlz (ai);
2228               goto mnzer1;
2229             }
2230         }
2231       eclear (c);
2232       return;
2233     }
2234  mnzer1:
2235
2236   if (bi[E] == 0)
2237     {
2238       for (i = 1; i < NI - 1; i++)
2239         {
2240           if (bi[i] != 0)
2241             {
2242               ltb -= enormlz (bi);
2243               goto mnzer2;
2244             }
2245         }
2246       eclear (c);
2247       return;
2248     }
2249  mnzer2:
2250
2251   /* Multiply significands */
2252   j = emulm (ai, bi);
2253   /* calculate exponent */
2254   lt = lta + ltb - (EXONE - 1);
2255   emdnorm (bi, j, 0, lt, 64);
2256   /* calculate sign of product */
2257   if (ai[0] == bi[0])
2258     bi[0] = 0;
2259   else
2260     bi[0] = 0xffff;
2261   emovo (bi, c);
2262 }
2263
2264
2265
2266
2267 /*
2268 ; Convert IEEE double precision to e type
2269 ;       double d;
2270 ;       unsigned EMUSHORT x[N+2];
2271 ;       e53toe (&d, x);
2272 */
2273 void
2274 e53toe (e, y)
2275      unsigned EMUSHORT *e, *y;
2276 {
2277 #ifdef DEC
2278
2279   dectoe (e, y);                /* see etodec.c */
2280
2281 #else
2282
2283   register unsigned EMUSHORT r;
2284   register unsigned EMUSHORT *p;
2285   unsigned EMUSHORT yy[NI];
2286   int denorm, k;
2287
2288   denorm = 0;                   /* flag if denormalized number */
2289   ecleaz (yy);
2290 #ifdef IBMPC
2291   e += 3;
2292 #endif
2293   r = *e;
2294   yy[0] = 0;
2295   if (r & 0x8000)
2296     yy[0] = 0xffff;
2297   yy[M] = (r & 0x0f) | 0x10;
2298   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2299 #ifdef INFINITY
2300   if (r == 0x7ff0)
2301     {
2302       einfin (y);
2303       if (r & 0x8000)
2304         eneg (y);
2305       return;
2306     }
2307 #endif
2308   r >>= 4;
2309   /* If zero exponent, then the significand is denormalized.
2310    * So, take back the understood high significand bit. */
2311   if (r == 0)
2312     {
2313       denorm = 1;
2314       yy[M] &= ~0x10;
2315     }
2316   r += EXONE - 01777;
2317   yy[E] = r;
2318   p = &yy[M + 1];
2319 #ifdef IBMPC
2320   *p++ = *(--e);
2321   *p++ = *(--e);
2322   *p++ = *(--e);
2323 #endif
2324 #ifdef MIEEE
2325   ++e;
2326   *p++ = *e++;
2327   *p++ = *e++;
2328   *p++ = *e++;
2329 #endif
2330   eshift (yy, -5);
2331   if (denorm)
2332     {                           /* if zero exponent, then normalize the significand */
2333       if ((k = enormlz (yy)) > NBITS)
2334         ecleazs (yy);
2335       else
2336         yy[E] -= (unsigned EMUSHORT) (k - 1);
2337     }
2338   emovo (yy, y);
2339 #endif /* not DEC */
2340 }
2341
2342 void
2343 e64toe (e, y)
2344      unsigned EMUSHORT *e, *y;
2345 {
2346   unsigned EMUSHORT yy[NI];
2347   unsigned EMUSHORT *p, *q;
2348   int i;
2349
2350   p = yy;
2351   for (i = 0; i < NE - 5; i++)
2352     *p++ = 0;
2353 #ifdef IBMPC
2354   for (i = 0; i < 5; i++)
2355     *p++ = *e++;
2356 #endif
2357 #ifdef DEC
2358   for (i = 0; i < 5; i++)
2359     *p++ = *e++;
2360 #endif
2361 #ifdef MIEEE
2362   p = &yy[0] + (NE - 1);
2363   *p-- = *e++;
2364   ++e;
2365   for (i = 0; i < 4; i++)
2366     *p-- = *e++;
2367 #endif
2368   p = yy;
2369   q = y;
2370 #ifdef INFINITY
2371   if (*p == 0x7fff)
2372     {
2373       einfin (y);
2374       if (*p & 0x8000)
2375         eneg (y);
2376       return;
2377     }
2378 #endif
2379   for (i = 0; i < NE; i++)
2380     *q++ = *p++;
2381 }
2382
2383
2384 /*
2385 ; Convert IEEE single precision to e type
2386 ;       float d;
2387 ;       unsigned EMUSHORT x[N+2];
2388 ;       dtox (&d, x);
2389 */
2390 void
2391 e24toe (e, y)
2392      unsigned EMUSHORT *e, *y;
2393 {
2394   register unsigned EMUSHORT r;
2395   register unsigned EMUSHORT *p;
2396   unsigned EMUSHORT yy[NI];
2397   int denorm, k;
2398
2399   denorm = 0;                   /* flag if denormalized number */
2400   ecleaz (yy);
2401 #ifdef IBMPC
2402   e += 1;
2403 #endif
2404 #ifdef DEC
2405   e += 1;
2406 #endif
2407   r = *e;
2408   yy[0] = 0;
2409   if (r & 0x8000)
2410     yy[0] = 0xffff;
2411   yy[M] = (r & 0x7f) | 0200;
2412   r &= ~0x807f;                 /* strip sign and 7 significand bits */
2413 #ifdef INFINITY
2414   if (r == 0x7f80)
2415     {
2416       einfin (y);
2417       if (r & 0x8000)
2418         eneg (y);
2419       return;
2420     }
2421 #endif
2422   r >>= 7;
2423   /* If zero exponent, then the significand is denormalized.
2424    * So, take back the understood high significand bit. */
2425   if (r == 0)
2426     {
2427       denorm = 1;
2428       yy[M] &= ~0200;
2429     }
2430   r += EXONE - 0177;
2431   yy[E] = r;
2432   p = &yy[M + 1];
2433 #ifdef IBMPC
2434   *p++ = *(--e);
2435 #endif
2436 #ifdef DEC
2437   *p++ = *(--e);
2438 #endif
2439 #ifdef MIEEE
2440   ++e;
2441   *p++ = *e++;
2442 #endif
2443   eshift (yy, -8);
2444   if (denorm)
2445     {                           /* if zero exponent, then normalize the significand */
2446       if ((k = enormlz (yy)) > NBITS)
2447         ecleazs (yy);
2448       else
2449         yy[E] -= (unsigned EMUSHORT) (k - 1);
2450     }
2451   emovo (yy, y);
2452 }
2453
2454
2455 void
2456 etoe64 (x, e)
2457      unsigned EMUSHORT *x, *e;
2458 {
2459   unsigned EMUSHORT xi[NI];
2460   EMULONG exp;
2461   int rndsav;
2462
2463   emovi (x, xi);
2464   /* adjust exponent for offset */
2465   exp = (EMULONG) xi[E];
2466 #ifdef INFINITY
2467   if (eisinf (x))
2468     goto nonorm;
2469 #endif
2470   /* round off to nearest or even */
2471   rndsav = rndprc;
2472   rndprc = 64;
2473   emdnorm (xi, 0, 0, exp, 64);
2474   rndprc = rndsav;
2475  nonorm:
2476   toe64 (xi, e);
2477 }
2478
2479 /* move out internal format to ieee long double */
2480 static void
2481 toe64 (a, b)
2482      unsigned EMUSHORT *a, *b;
2483 {
2484   register unsigned EMUSHORT *p, *q;
2485   unsigned EMUSHORT i;
2486
2487   p = a;
2488 #ifdef MIEEE
2489   q = b;
2490 #else
2491   q = b + 4;                    /* point to output exponent */
2492 #if LONG_DOUBLE_TYPE_SIZE == 96
2493   /* Clear the last two bytes of 12-byte Intel format */
2494   *(q+1) = 0;
2495 #endif
2496 #endif
2497
2498   /* combine sign and exponent */
2499   i = *p++;
2500 #ifdef MIEEE
2501   if (i)
2502     *q++ = *p++ | 0x8000;
2503   else
2504     *q++ = *p++;
2505   *q++ = 0;
2506 #else
2507   if (i)
2508     *q-- = *p++ | 0x8000;
2509   else
2510     *q-- = *p++;
2511 #endif
2512   /* skip over guard word */
2513   ++p;
2514   /* move the significand */
2515 #ifdef MIEEE
2516   for (i = 0; i < 4; i++)
2517     *q++ = *p++;
2518 #else
2519   for (i = 0; i < 4; i++)
2520     *q-- = *p++;
2521 #endif
2522 }
2523
2524
2525 /*
2526 ; e type to IEEE double precision
2527 ;       double d;
2528 ;       unsigned EMUSHORT x[NE];
2529 ;       etoe53 (x, &d);
2530 */
2531
2532 #ifdef DEC
2533
2534 void
2535 etoe53 (x, e)
2536      unsigned EMUSHORT *x, *e;
2537 {
2538   etodec (x, e);                /* see etodec.c */
2539 }
2540
2541 static void
2542 toe53 (x, y)
2543      unsigned EMUSHORT *x, *y;
2544 {
2545   todec (x, y);
2546 }
2547
2548 #else
2549
2550 void
2551 etoe53 (x, e)
2552      unsigned EMUSHORT *x, *e;
2553 {
2554   unsigned EMUSHORT xi[NI];
2555   EMULONG exp;
2556   int rndsav;
2557
2558   emovi (x, xi);
2559   /* adjust exponent for offsets */
2560   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2561 #ifdef INFINITY
2562   if (eisinf (x))
2563     goto nonorm;
2564 #endif
2565   /* round off to nearest or even */
2566   rndsav = rndprc;
2567   rndprc = 53;
2568   emdnorm (xi, 0, 0, exp, 64);
2569   rndprc = rndsav;
2570  nonorm:
2571   toe53 (xi, e);
2572 }
2573
2574
2575 static void
2576 toe53 (x, y)
2577      unsigned EMUSHORT *x, *y;
2578 {
2579   unsigned EMUSHORT i;
2580   unsigned EMUSHORT *p;
2581
2582
2583   p = &x[0];
2584 #ifdef IBMPC
2585   y += 3;
2586 #endif
2587   *y = 0;                       /* output high order */
2588   if (*p++)
2589     *y = 0x8000;                /* output sign bit */
2590
2591   i = *p++;
2592   if (i >= (unsigned int) 2047)
2593     {                           /* Saturate at largest number less than infinity. */
2594 #ifdef INFINITY
2595       *y |= 0x7ff0;
2596 #ifdef IBMPC
2597       *(--y) = 0;
2598       *(--y) = 0;
2599       *(--y) = 0;
2600 #endif
2601 #ifdef MIEEE
2602       ++y;
2603       *y++ = 0;
2604       *y++ = 0;
2605       *y++ = 0;
2606 #endif
2607 #else
2608       *y |= (unsigned EMUSHORT) 0x7fef;
2609 #ifdef IBMPC
2610       *(--y) = 0xffff;
2611       *(--y) = 0xffff;
2612       *(--y) = 0xffff;
2613 #endif
2614 #ifdef MIEEE
2615       ++y;
2616       *y++ = 0xffff;
2617       *y++ = 0xffff;
2618       *y++ = 0xffff;
2619 #endif
2620 #endif
2621       return;
2622     }
2623   if (i == 0)
2624     {
2625       eshift (x, 4);
2626     }
2627   else
2628     {
2629       i <<= 4;
2630       eshift (x, 5);
2631     }
2632   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
2633   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
2634 #ifdef IBMPC
2635   *(--y) = *p++;
2636   *(--y) = *p++;
2637   *(--y) = *p;
2638 #endif
2639 #ifdef MIEEE
2640   ++y;
2641   *y++ = *p++;
2642   *y++ = *p++;
2643   *y++ = *p++;
2644 #endif
2645 }
2646
2647 #endif /* not DEC */
2648
2649
2650
2651 /*
2652 ; e type to IEEE single precision
2653 ;       float d;
2654 ;       unsigned EMUSHORT x[N+2];
2655 ;       xtod (x, &d);
2656 */
2657 void
2658 etoe24 (x, e)
2659      unsigned EMUSHORT *x, *e;
2660 {
2661   EMULONG exp;
2662   unsigned EMUSHORT xi[NI];
2663   int rndsav;
2664
2665   emovi (x, xi);
2666   /* adjust exponent for offsets */
2667   exp = (EMULONG) xi[E] - (EXONE - 0177);
2668 #ifdef INFINITY
2669   if (eisinf (x))
2670     goto nonorm;
2671 #endif
2672   /* round off to nearest or even */
2673   rndsav = rndprc;
2674   rndprc = 24;
2675   emdnorm (xi, 0, 0, exp, 64);
2676   rndprc = rndsav;
2677  nonorm:
2678   toe24 (xi, e);
2679 }
2680
2681 static void
2682 toe24 (x, y)
2683      unsigned EMUSHORT *x, *y;
2684 {
2685   unsigned EMUSHORT i;
2686   unsigned EMUSHORT *p;
2687
2688   p = &x[0];
2689 #ifdef IBMPC
2690   y += 1;
2691 #endif
2692 #ifdef DEC
2693   y += 1;
2694 #endif
2695   *y = 0;                       /* output high order */
2696   if (*p++)
2697     *y = 0x8000;                /* output sign bit */
2698
2699   i = *p++;
2700 /* Handle overflow cases. */
2701   if (i >= 255)
2702     {
2703 #ifdef INFINITY
2704       *y |= (unsigned EMUSHORT) 0x7f80;
2705 #ifdef IBMPC
2706       *(--y) = 0;
2707 #endif
2708 #ifdef DEC
2709       *(--y) = 0;
2710 #endif
2711 #ifdef MIEEE
2712       ++y;
2713       *y = 0;
2714 #endif
2715 #else  /* no INFINITY */
2716       *y |= (unsigned EMUSHORT) 0x7f7f;
2717 #ifdef IBMPC
2718       *(--y) = 0xffff;
2719 #endif
2720 #ifdef DEC
2721       *(--y) = 0xffff;
2722 #endif
2723 #ifdef MIEEE
2724       ++y;
2725       *y = 0xffff;
2726 #endif
2727 #ifdef ERANGE
2728       errno = ERANGE;
2729 #endif
2730 #endif  /* no INFINITY */
2731       return;
2732     }
2733   if (i == 0)
2734     {
2735       eshift (x, 7);
2736     }
2737   else
2738     {
2739       i <<= 7;
2740       eshift (x, 8);
2741     }
2742   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
2743   *y |= i;                      /* high order output already has sign bit set */
2744 #ifdef IBMPC
2745   *(--y) = *p;
2746 #endif
2747 #ifdef DEC
2748   *(--y) = *p;
2749 #endif
2750 #ifdef MIEEE
2751   ++y;
2752   *y = *p;
2753 #endif
2754 }
2755
2756
2757 /* Compare two e type numbers.
2758  *
2759  * unsigned EMUSHORT a[NE], b[NE];
2760  * ecmp (a, b);
2761  *
2762  *  returns +1 if a > b
2763  *           0 if a == b
2764  *          -1 if a < b
2765  */
2766 int
2767 ecmp (a, b)
2768      unsigned EMUSHORT *a, *b;
2769 {
2770   unsigned EMUSHORT ai[NI], bi[NI];
2771   register unsigned EMUSHORT *p, *q;
2772   register int i;
2773   int msign;
2774
2775   emovi (a, ai);
2776   p = ai;
2777   emovi (b, bi);
2778   q = bi;
2779
2780   if (*p != *q)
2781     {                           /* the signs are different */
2782       /* -0 equals + 0 */
2783       for (i = 1; i < NI - 1; i++)
2784         {
2785           if (ai[i] != 0)
2786             goto nzro;
2787           if (bi[i] != 0)
2788             goto nzro;
2789         }
2790       return (0);
2791     nzro:
2792       if (*p == 0)
2793         return (1);
2794       else
2795         return (-1);
2796     }
2797   /* both are the same sign */
2798   if (*p == 0)
2799     msign = 1;
2800   else
2801     msign = -1;
2802   i = NI - 1;
2803   do
2804     {
2805       if (*p++ != *q++)
2806         {
2807           goto diff;
2808         }
2809     }
2810   while (--i > 0);
2811
2812   return (0);                   /* equality */
2813
2814
2815
2816  diff:
2817
2818   if (*(--p) > *(--q))
2819     return (msign);             /* p is bigger */
2820   else
2821     return (-msign);            /* p is littler */
2822 }
2823
2824
2825
2826
2827 /* Find nearest integer to x = floor (x + 0.5)
2828  *
2829  * unsigned EMUSHORT x[NE], y[NE]
2830  * eround (x, y);
2831  */
2832 void
2833 eround (x, y)
2834      unsigned EMUSHORT *x, *y;
2835 {
2836   eadd (ehalf, x, y);
2837   efloor (y, y);
2838 }
2839
2840
2841
2842
2843 /*
2844 ; convert long integer to e type
2845 ;
2846 ;       long l;
2847 ;       unsigned EMUSHORT x[NE];
2848 ;       ltoe (&l, x);
2849 ; note &l is the memory address of l
2850 */
2851 void
2852 ltoe (lp, y)
2853      long *lp;                  /* lp is the memory address of a long integer */
2854      unsigned EMUSHORT *y;              /* y is the address of a short */
2855 {
2856   unsigned EMUSHORT yi[NI];
2857   unsigned long ll;
2858   int k;
2859
2860   ecleaz (yi);
2861   if (*lp < 0)
2862     {
2863       /* make it positive */
2864       ll = (unsigned long) (-(*lp));
2865       yi[0] = 0xffff;           /* put correct sign in the e type number */
2866     }
2867   else
2868     {
2869       ll = (unsigned long) (*lp);
2870     }
2871   /* move the long integer to yi significand area */
2872   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2873   yi[M + 1] = (unsigned EMUSHORT) ll;
2874
2875   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2876   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2877     ecleaz (yi);                /* it was zero */
2878   else
2879     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
2880   emovo (yi, y);                /* output the answer */
2881 }
2882
2883 /*
2884 ; convert unsigned long integer to e type
2885 ;
2886 ;       unsigned long l;
2887 ;       unsigned EMUSHORT x[NE];
2888 ;       ltox (&l, x);
2889 ; note &l is the memory address of l
2890 */
2891 void
2892 ultoe (lp, y)
2893      unsigned long *lp;         /* lp is the memory address of a long integer */
2894      unsigned EMUSHORT *y;              /* y is the address of a short */
2895 {
2896   unsigned EMUSHORT yi[NI];
2897   unsigned long ll;
2898   int k;
2899
2900   ecleaz (yi);
2901   ll = *lp;
2902
2903   /* move the long integer to ayi significand area */
2904   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2905   yi[M + 1] = (unsigned EMUSHORT) ll;
2906
2907   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2908   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2909     ecleaz (yi);                /* it was zero */
2910   else
2911     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
2912   emovo (yi, y);                /* output the answer */
2913 }
2914
2915
2916 /*
2917 ;       Find long integer and fractional parts
2918
2919 ;       long i;
2920 ;       unsigned EMUSHORT x[NE], frac[NE];
2921 ;       xifrac (x, &i, frac);
2922
2923   The integer output has the sign of the input.  The fraction is
2924 the positive fractional part of abs (x).
2925 */
2926 void
2927 eifrac (x, i, frac)
2928      unsigned EMUSHORT *x;
2929      long *i;
2930      unsigned EMUSHORT *frac;
2931 {
2932   unsigned EMUSHORT xi[NI];
2933   int k;
2934
2935   emovi (x, xi);
2936   k = (int) xi[E] - (EXONE - 1);
2937   if (k <= 0)
2938     {
2939       /* if exponent <= 0, integer = 0 and real output is fraction */
2940       *i = 0L;
2941       emovo (xi, frac);
2942       return;
2943     }
2944   if (k > (HOST_BITS_PER_LONG - 1))
2945     {
2946       /*
2947          ;      long integer overflow: output large integer
2948          ;      and correct fraction
2949          */
2950       if (xi[0])
2951         *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
2952       else
2953         *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
2954       eshift (xi, k);
2955       if (extra_warnings)
2956         warning ("overflow on truncation to integer");
2957       goto lab11;
2958     }
2959
2960   if (k > 16)
2961     {
2962       /*
2963          ; shift more than 16 bits: shift up k-16, output the integer,
2964          ; then complete the shift to get the fraction.
2965          */
2966       k -= 16;
2967       eshift (xi, k);
2968
2969       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
2970       eshup6 (xi);
2971       goto lab10;
2972     }
2973
2974   /* shift not more than 16 bits */
2975   eshift (xi, k);
2976   *i = (long) xi[M] & 0xffff;
2977
2978  lab10:
2979
2980   if (xi[0])
2981     *i = -(*i);
2982  lab11:
2983
2984   xi[0] = 0;
2985   xi[E] = EXONE - 1;
2986   xi[M] = 0;
2987   if ((k = enormlz (xi)) > NBITS)
2988     ecleaz (xi);
2989   else
2990     xi[E] -= (unsigned EMUSHORT) k;
2991
2992   emovo (xi, frac);
2993 }
2994
2995
2996 /*
2997 ;       Find unsigned long integer and fractional parts
2998
2999 ;       unsigned long i;
3000 ;       unsigned EMUSHORT x[NE], frac[NE];
3001 ;       xifrac (x, &i, frac);
3002
3003   A negative e type input yields integer output = 0
3004   but correct fraction.
3005 */
3006 void
3007 euifrac (x, i, frac)
3008      unsigned EMUSHORT *x;
3009      long *i;
3010      unsigned EMUSHORT *frac;
3011 {
3012   unsigned EMUSHORT xi[NI];
3013   int k;
3014
3015   emovi (x, xi);
3016   k = (int) xi[E] - (EXONE - 1);
3017   if (k <= 0)
3018     {
3019       /* if exponent <= 0, integer = 0 and argument is fraction */
3020       *i = 0L;
3021       emovo (xi, frac);
3022       return;
3023     }
3024   if (k > 32)
3025     {
3026       /*
3027          ;      long integer overflow: output large integer
3028          ;      and correct fraction
3029          */
3030       *i = ~(0L);
3031       eshift (xi, k);
3032       if (extra_warnings)
3033         warning ("overflow on truncation to unsigned integer");
3034       goto lab10;
3035     }
3036
3037   if (k > 16)
3038     {
3039       /*
3040          ; shift more than 16 bits: shift up k-16, output the integer,
3041          ; then complete the shift to get the fraction.
3042          */
3043       k -= 16;
3044       eshift (xi, k);
3045
3046       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
3047       eshup6 (xi);
3048       goto lab10;
3049     }
3050
3051   /* shift not more than 16 bits */
3052   eshift (xi, k);
3053   *i = (long) xi[M] & 0xffff;
3054
3055  lab10:
3056
3057   if (xi[0])
3058     *i = 0L;
3059
3060   xi[0] = 0;
3061   xi[E] = EXONE - 1;
3062   xi[M] = 0;
3063   if ((k = enormlz (xi)) > NBITS)
3064     ecleaz (xi);
3065   else
3066     xi[E] -= (unsigned EMUSHORT) k;
3067
3068   emovo (xi, frac);
3069 }
3070
3071
3072
3073 /*
3074 ;       Shift significand
3075 ;
3076 ;       Shifts significand area up or down by the number of bits
3077 ;       given by the variable sc.
3078 */
3079 int
3080 eshift (x, sc)
3081      unsigned EMUSHORT *x;
3082      int sc;
3083 {
3084   unsigned EMUSHORT lost;
3085   unsigned EMUSHORT *p;
3086
3087   if (sc == 0)
3088     return (0);
3089
3090   lost = 0;
3091   p = x + NI - 1;
3092
3093   if (sc < 0)
3094     {
3095       sc = -sc;
3096       while (sc >= 16)
3097         {
3098           lost |= *p;           /* remember lost bits */
3099           eshdn6 (x);
3100           sc -= 16;
3101         }
3102
3103       while (sc >= 8)
3104         {
3105           lost |= *p & 0xff;
3106           eshdn8 (x);
3107           sc -= 8;
3108         }
3109
3110       while (sc > 0)
3111         {
3112           lost |= *p & 1;
3113           eshdn1 (x);
3114           sc -= 1;
3115         }
3116     }
3117   else
3118     {
3119       while (sc >= 16)
3120         {
3121           eshup6 (x);
3122           sc -= 16;
3123         }
3124
3125       while (sc >= 8)
3126         {
3127           eshup8 (x);
3128           sc -= 8;
3129         }
3130
3131       while (sc > 0)
3132         {
3133           eshup1 (x);
3134           sc -= 1;
3135         }
3136     }
3137   if (lost)
3138     lost = 1;
3139   return ((int) lost);
3140 }
3141
3142
3143
3144 /*
3145 ;       normalize
3146 ;
3147 ; Shift normalizes the significand area pointed to by argument
3148 ; shift count (up = positive) is returned.
3149 */
3150 int
3151 enormlz (x)
3152      unsigned EMUSHORT x[];
3153 {
3154   register unsigned EMUSHORT *p;
3155   int sc;
3156
3157   sc = 0;
3158   p = &x[M];
3159   if (*p != 0)
3160     goto normdn;
3161   ++p;
3162   if (*p & 0x8000)
3163     return (0);                 /* already normalized */
3164   while (*p == 0)
3165     {
3166       eshup6 (x);
3167       sc += 16;
3168       /* With guard word, there are NBITS+16 bits available.
3169        * return true if all are zero.
3170        */
3171       if (sc > NBITS)
3172         return (sc);
3173     }
3174   /* see if high byte is zero */
3175   while ((*p & 0xff00) == 0)
3176     {
3177       eshup8 (x);
3178       sc += 8;
3179     }
3180   /* now shift 1 bit at a time */
3181   while ((*p & 0x8000) == 0)
3182     {
3183       eshup1 (x);
3184       sc += 1;
3185       if (sc > NBITS)
3186         {
3187           mtherr ("enormlz", UNDERFLOW);
3188           return (sc);
3189         }
3190     }
3191   return (sc);
3192
3193   /* Normalize by shifting down out of the high guard word
3194      of the significand */
3195  normdn:
3196
3197   if (*p & 0xff00)
3198     {
3199       eshdn8 (x);
3200       sc -= 8;
3201     }
3202   while (*p != 0)
3203     {
3204       eshdn1 (x);
3205       sc -= 1;
3206
3207       if (sc < -NBITS)
3208         {
3209           mtherr ("enormlz", OVERFLOW);
3210           return (sc);
3211         }
3212     }
3213   return (sc);
3214 }
3215
3216
3217
3218
3219 /* Convert e type number to decimal format ASCII string.
3220  * The constants are for 64 bit precision.
3221  */
3222
3223 #define NTEN 12
3224 #define MAXP 4096
3225
3226 static unsigned EMUSHORT etens[NTEN + 1][NE] =
3227 {
3228   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
3229   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
3230   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3231   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3232   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3233   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3234   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3235   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3236   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3237   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3238   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3239   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3240   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
3241 };
3242
3243 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
3244 {
3245   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
3246   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
3247   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3248   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3249   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3250   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3251   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3252   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3253   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3254   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3255   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3256   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3257   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
3258 };
3259
3260 void
3261 e24toasc (x, string, ndigs)
3262      unsigned EMUSHORT x[];
3263      char *string;
3264      int ndigs;
3265 {
3266   unsigned EMUSHORT w[NI];
3267
3268 #ifdef INFINITY
3269 #ifdef IBMPC
3270   if ((x[1] & 0x7f80) == 0x7f80)
3271 #else
3272   if ((x[0] & 0x7f80) == 0x7f80)
3273 #endif
3274     {
3275 #ifdef IBMPC
3276       if (x[1] & 0x8000)
3277 #else
3278       if (x[0] & 0x8000)
3279 #endif
3280         sprintf (string, " -Infinity ");
3281       else
3282         sprintf (string, " Infinity ");
3283       return;
3284     }
3285 #endif
3286   e24toe (x, w);
3287   etoasc (w, string, ndigs);
3288 }
3289
3290
3291 void
3292 e53toasc (x, string, ndigs)
3293      unsigned EMUSHORT x[];
3294      char *string;
3295      int ndigs;
3296 {
3297   unsigned EMUSHORT w[NI];
3298
3299 #ifdef INFINITY
3300 #ifdef IBMPC
3301   if ((x[3] & 0x7ff0) == 0x7ff0)
3302 #else
3303   if ((x[0] & 0x7ff0) == 0x7ff0)
3304 #endif
3305     {
3306 #ifdef IBMPC
3307       if (x[3] & 0x8000)
3308 #else
3309       if (x[0] & 0x8000)
3310 #endif
3311         sprintf (string, " -Infinity ");
3312       else
3313         sprintf (string, " Infinity ");
3314       return;
3315     }
3316 #endif
3317   e53toe (x, w);
3318   etoasc (w, string, ndigs);
3319 }
3320
3321
3322 void
3323 e64toasc (x, string, ndigs)
3324      unsigned EMUSHORT x[];
3325      char *string;
3326      int ndigs;
3327 {
3328   unsigned EMUSHORT w[NI];
3329
3330 #ifdef INFINITY
3331 #ifdef IBMPC
3332   if ((x[4] & 0x7fff) == 0x7fff)
3333 #else
3334   if ((x[0] & 0x7fff) == 0x7fff)
3335 #endif
3336     {
3337 #ifdef IBMPC
3338       if (x[4] & 0x8000)
3339 #else
3340       if (x[0] & 0x8000)
3341 #endif
3342         sprintf (string, " -Infinity ");
3343       else
3344         sprintf (string, " Infinity ");
3345       return;
3346     }
3347 #endif
3348   e64toe (x, w);
3349   etoasc (w, string, ndigs);
3350 }
3351
3352
3353 static char wstring[80];        /* working storage for ASCII output */
3354
3355 void
3356 etoasc (x, string, ndigs)
3357      unsigned EMUSHORT x[];
3358      char *string;
3359      int ndigs;
3360 {
3361   EMUSHORT digit;
3362   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
3363   unsigned EMUSHORT *p, *r, *ten;
3364   unsigned EMUSHORT sign;
3365   int i, j, k, expon, rndsav;
3366   char *s, *ss;
3367   unsigned EMUSHORT m;
3368
3369   ss = string;
3370   s = wstring;
3371   while ((*s++ = *ss++) != '\0')
3372     ;
3373   rndsav = rndprc;
3374   rndprc = NBITS;               /* set to full precision */
3375   emov (x, y);                  /* retain external format */
3376   if (y[NE - 1] & 0x8000)
3377     {
3378       sign = 0xffff;
3379       y[NE - 1] &= 0x7fff;
3380     }
3381   else
3382     {
3383       sign = 0;
3384     }
3385   expon = 0;
3386   ten = &etens[NTEN][0];
3387   emov (eone, t);
3388   /* Test for zero exponent */
3389   if (y[NE - 1] == 0)
3390     {
3391       for (k = 0; k < NE - 1; k++)
3392         {
3393           if (y[k] != 0)
3394             goto tnzro;         /* denormalized number */
3395         }
3396       goto isone;               /* legal all zeros */
3397     }
3398  tnzro:
3399
3400   /* Test for infinity.  Don't bother with illegal infinities.
3401    */
3402   if (y[NE - 1] == 0x7fff)
3403     {
3404       if (sign)
3405         sprintf (wstring, " -Infinity ");
3406       else
3407         sprintf (wstring, " Infinity ");
3408       goto bxit;
3409     }
3410
3411   /* Test for exponent nonzero but significand denormalized.
3412    * This is an error condition.
3413    */
3414   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3415     {
3416       mtherr ("etoasc", DOMAIN);
3417       sprintf (wstring, "NaN");
3418       goto bxit;
3419     }
3420
3421   /* Compare to 1.0 */
3422   i = ecmp (eone, y);
3423   if (i == 0)
3424     goto isone;
3425
3426   if (i < 0)
3427     {                           /* Number is greater than 1 */
3428       /* Convert significand to an integer and strip trailing decimal zeros. */
3429       emov (y, u);
3430       u[NE - 1] = EXONE + NBITS - 1;
3431
3432       p = &etens[NTEN - 4][0];
3433       m = 16;
3434       do
3435         {
3436           ediv (p, u, t);
3437           efloor (t, w);
3438           for (j = 0; j < NE - 1; j++)
3439             {
3440               if (t[j] != w[j])
3441                 goto noint;
3442             }
3443           emov (t, u);
3444           expon += (int) m;
3445         noint:
3446           p += NE;
3447           m >>= 1;
3448         }
3449       while (m != 0);
3450
3451       /* Rescale from integer significand */
3452       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3453       emov (u, y);
3454       /* Find power of 10 */
3455       emov (eone, t);
3456       m = MAXP;
3457       p = &etens[0][0];
3458       while (ecmp (ten, u) <= 0)
3459         {
3460           if (ecmp (p, u) <= 0)
3461             {
3462               ediv (p, u, u);
3463               emul (p, t, t);
3464               expon += (int) m;
3465             }
3466           m >>= 1;
3467           if (m == 0)
3468             break;
3469           p += NE;
3470         }
3471     }
3472   else
3473     {                           /* Number is less than 1.0 */
3474       /* Pad significand with trailing decimal zeros. */
3475       if (y[NE - 1] == 0)
3476         {
3477           while ((y[NE - 2] & 0x8000) == 0)
3478             {
3479               emul (ten, y, y);
3480               expon -= 1;
3481             }
3482         }
3483       else
3484         {
3485           emovi (y, w);
3486           for (i = 0; i < NDEC + 1; i++)
3487             {
3488               if ((w[NI - 1] & 0x7) != 0)
3489                 break;
3490               /* multiply by 10 */
3491               emovz (w, u);
3492               eshdn1 (u);
3493               eshdn1 (u);
3494               eaddm (w, u);
3495               u[1] += 3;
3496               while (u[2] != 0)
3497                 {
3498                   eshdn1 (u);
3499                   u[1] += 1;
3500                 }
3501               if (u[NI - 1] != 0)
3502                 break;
3503               if (eone[NE - 1] <= u[1])
3504                 break;
3505               emovz (u, w);
3506               expon -= 1;
3507             }
3508           emovo (w, y);
3509         }
3510       k = -MAXP;
3511       p = &emtens[0][0];
3512       r = &etens[0][0];
3513       emov (y, w);
3514       emov (eone, t);
3515       while (ecmp (eone, w) > 0)
3516         {
3517           if (ecmp (p, w) >= 0)
3518             {
3519               emul (r, w, w);
3520               emul (r, t, t);
3521               expon += k;
3522             }
3523           k /= 2;
3524           if (k == 0)
3525             break;
3526           p += NE;
3527           r += NE;
3528         }
3529       ediv (t, eone, t);
3530     }
3531  isone:
3532   /* Find the first (leading) digit. */
3533   emovi (t, w);
3534   emovz (w, t);
3535   emovi (y, w);
3536   emovz (w, y);
3537   eiremain (t, y);
3538   digit = equot[NI - 1];
3539   while ((digit == 0) && (ecmp (y, ezero) != 0))
3540     {
3541       eshup1 (y);
3542       emovz (y, u);
3543       eshup1 (u);
3544       eshup1 (u);
3545       eaddm (u, y);
3546       eiremain (t, y);
3547       digit = equot[NI - 1];
3548       expon -= 1;
3549     }
3550   s = wstring;
3551   if (sign)
3552     *s++ = '-';
3553   else
3554     *s++ = ' ';
3555   /* Examine number of digits requested by caller. */
3556   if (ndigs < 0)
3557     ndigs = 0;
3558   if (ndigs > NDEC)
3559     ndigs = NDEC;
3560   if (digit == 10)
3561     {
3562       *s++ = '1';
3563       *s++ = '.';
3564       if (ndigs > 0)
3565         {
3566           *s++ = '0';
3567           ndigs -= 1;
3568         }
3569       expon += 1;
3570     }
3571   else
3572     {
3573       *s++ = (char )digit + '0';
3574       *s++ = '.';
3575     }
3576   /* Generate digits after the decimal point. */
3577   for (k = 0; k <= ndigs; k++)
3578     {
3579       /* multiply current number by 10, without normalizing */
3580       eshup1 (y);
3581       emovz (y, u);
3582       eshup1 (u);
3583       eshup1 (u);
3584       eaddm (u, y);
3585       eiremain (t, y);
3586       *s++ = (char) equot[NI - 1] + '0';
3587     }
3588   digit = equot[NI - 1];
3589   --s;
3590   ss = s;
3591   /* round off the ASCII string */
3592   if (digit > 4)
3593     {
3594       /* Test for critical rounding case in ASCII output. */
3595       if (digit == 5)
3596         {
3597           emovo (y, t);
3598           if (ecmp (t, ezero) != 0)
3599             goto roun;          /* round to nearest */
3600           if ((*(s - 1) & 1) == 0)
3601             goto doexp;         /* round to even */
3602         }
3603       /* Round up and propagate carry-outs */
3604     roun:
3605       --s;
3606       k = *s & 0x7f;
3607       /* Carry out to most significant digit? */
3608       if (k == '.')
3609         {
3610           --s;
3611           k = *s;
3612           k += 1;
3613           *s = (char) k;
3614           /* Most significant digit carries to 10? */
3615           if (k > '9')
3616             {
3617               expon += 1;
3618               *s = '1';
3619             }
3620           goto doexp;
3621         }
3622       /* Round up and carry out from less significant digits */
3623       k += 1;
3624       *s = (char) k;
3625       if (k > '9')
3626         {
3627           *s = '0';
3628           goto roun;
3629         }
3630     }
3631  doexp:
3632   /*
3633      if (expon >= 0)
3634      sprintf (ss, "e+%d", expon);
3635      else
3636      sprintf (ss, "e%d", expon);
3637      */
3638   sprintf (ss, "e%d", expon);
3639  bxit:
3640   rndprc = rndsav;
3641   /* copy out the working string */
3642   s = string;
3643   ss = wstring;
3644   while (*ss == ' ')            /* strip possible leading space */
3645     ++ss;
3646   while ((*s++ = *ss++) != '\0')
3647     ;
3648 }
3649
3650
3651
3652
3653 /*
3654 ;                                                               ASCTOQ
3655 ;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
3656 ;                                       SLM, 3 JAN 78
3657 ;
3658 ;       Convert ASCII string to quadruple precision floating point
3659 ;
3660 ;               Numeric input is free field decimal number
3661 ;               with max of 15 digits with or without
3662 ;               decimal point entered as ASCII from teletype.
3663 ;       Entering E after the number followed by a second
3664 ;       number causes the second number to be interpreted
3665 ;       as a power of 10 to be multiplied by the first number
3666 ;       (i.e., "scientific" notation).
3667 ;
3668 ;       Usage:
3669 ;               asctoq (string, q);
3670 */
3671
3672 /* ASCII to single */
3673 void
3674 asctoe24 (s, y)
3675      char *s;
3676      unsigned EMUSHORT *y;
3677 {
3678   asctoeg (s, y, 24);
3679 }
3680
3681
3682 /* ASCII to double */
3683 void
3684 asctoe53 (s, y)
3685      char *s;
3686      unsigned EMUSHORT *y;
3687 {
3688 #ifdef DEC
3689   asctoeg (s, y, 56);
3690 #else
3691   asctoeg (s, y, 53);
3692 #endif
3693 }
3694
3695
3696 /* ASCII to long double */
3697 void
3698 asctoe64 (s, y)
3699      char *s;
3700      unsigned EMUSHORT *y;
3701 {
3702   asctoeg (s, y, 64);
3703 }
3704
3705 /* ASCII to super double */
3706 void
3707 asctoe (s, y)
3708      char *s;
3709      unsigned EMUSHORT *y;
3710 {
3711   asctoeg (s, y, NBITS);
3712 }
3713
3714 /* Space to make a copy of the input string: */
3715 static char lstr[82];
3716
3717 void
3718 asctoeg (ss, y, oprec)
3719      char *ss;
3720      unsigned EMUSHORT *y;
3721      int oprec;
3722 {
3723   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
3724   int esign, decflg, sgnflg, nexp, exp, prec, lost;
3725   int k, trail, c, rndsav;
3726   EMULONG lexp;
3727   unsigned EMUSHORT nsign, *p;
3728   char *sp, *s;
3729
3730   /* Copy the input string. */
3731   s = ss;
3732   while (*s == ' ')             /* skip leading spaces */
3733     ++s;
3734   sp = lstr;
3735   for (k = 0; k < 79; k++)
3736     {
3737       if ((*sp++ = *s++) == '\0')
3738         break;
3739     }
3740   *sp = '\0';
3741   s = lstr;
3742
3743   rndsav = rndprc;
3744   rndprc = NBITS;               /* Set to full precision */
3745   lost = 0;
3746   nsign = 0;
3747   decflg = 0;
3748   sgnflg = 0;
3749   nexp = 0;
3750   exp = 0;
3751   prec = 0;
3752   ecleaz (yy);
3753   trail = 0;
3754
3755  nxtcom:
3756   k = *s - '0';
3757   if ((k >= 0) && (k <= 9))
3758     {
3759       /* Ignore leading zeros */
3760       if ((prec == 0) && (decflg == 0) && (k == 0))
3761         goto donchr;
3762       /* Identify and strip trailing zeros after the decimal point. */
3763       if ((trail == 0) && (decflg != 0))
3764         {
3765           sp = s;
3766           while ((*sp >= '0') && (*sp <= '9'))
3767             ++sp;
3768           /* Check for syntax error */
3769           c = *sp & 0x7f;
3770           if ((c != 'e') && (c != 'E') && (c != '\0')
3771               && (c != '\n') && (c != '\r') && (c != ' ')
3772               && (c != ','))
3773             goto error;
3774           --sp;
3775           while (*sp == '0')
3776             *sp-- = 'z';
3777           trail = 1;
3778           if (*s == 'z')
3779             goto donchr;
3780         }
3781       /* If enough digits were given to more than fill up the yy register,
3782        * continuing until overflow into the high guard word yy[2]
3783        * guarantees that there will be a roundoff bit at the top
3784        * of the low guard word after normalization.
3785        */
3786       if (yy[2] == 0)
3787         {
3788           if (decflg)
3789             nexp += 1;          /* count digits after decimal point */
3790           eshup1 (yy);          /* multiply current number by 10 */
3791           emovz (yy, xt);
3792           eshup1 (xt);
3793           eshup1 (xt);
3794           eaddm (xt, yy);
3795           ecleaz (xt);
3796           xt[NI - 2] = (unsigned EMUSHORT) k;
3797           eaddm (xt, yy);
3798         }
3799       else
3800         {
3801           lost |= k;
3802         }
3803       prec += 1;
3804       goto donchr;
3805     }
3806
3807   switch (*s)
3808     {
3809     case 'z':
3810       break;
3811     case 'E':
3812     case 'e':
3813       goto expnt;
3814     case '.':                   /* decimal point */
3815       if (decflg)
3816         goto error;
3817       ++decflg;
3818       break;
3819     case '-':
3820       nsign = 0xffff;
3821       if (sgnflg)
3822         goto error;
3823       ++sgnflg;
3824       break;
3825     case '+':
3826       if (sgnflg)
3827         goto error;
3828       ++sgnflg;
3829       break;
3830     case ',':
3831     case ' ':
3832     case '\0':
3833     case '\n':
3834     case '\r':
3835       goto daldone;
3836     case 'i':
3837     case 'I':
3838       goto infinite;
3839     default:
3840     error:
3841       mtherr ("asctoe", DOMAIN);
3842       eclear (y);
3843       goto aexit;
3844     }
3845  donchr:
3846   ++s;
3847   goto nxtcom;
3848
3849   /* Exponent interpretation */
3850  expnt:
3851
3852   esign = 1;
3853   exp = 0;
3854   ++s;
3855   /* check for + or - */
3856   if (*s == '-')
3857     {
3858       esign = -1;
3859       ++s;
3860     }
3861   if (*s == '+')
3862     ++s;
3863   while ((*s >= '0') && (*s <= '9'))
3864     {
3865       exp *= 10;
3866       exp += *s++ - '0';
3867       if (exp > 4956)
3868         {
3869           if (esign < 0)
3870             goto zero;
3871           else
3872             goto infinite;
3873         }
3874     }
3875   if (esign < 0)
3876     exp = -exp;
3877   if (exp > 4932)
3878     {
3879  infinite:
3880       ecleaz (yy);
3881       yy[E] = 0x7fff;           /* infinity */
3882       goto aexit;
3883     }
3884   if (exp < -4956)
3885     {
3886  zero:
3887       ecleaz (yy);
3888       goto aexit;
3889     }
3890
3891  daldone:
3892   nexp = exp - nexp;
3893   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3894   while ((nexp > 0) && (yy[2] == 0))
3895     {
3896       emovz (yy, xt);
3897       eshup1 (xt);
3898       eshup1 (xt);
3899       eaddm (yy, xt);
3900       eshup1 (xt);
3901       if (xt[2] != 0)
3902         break;
3903       nexp -= 1;
3904       emovz (xt, yy);
3905     }
3906   if ((k = enormlz (yy)) > NBITS)
3907     {
3908       ecleaz (yy);
3909       goto aexit;
3910     }
3911   lexp = (EXONE - 1 + NBITS) - k;
3912   emdnorm (yy, lost, 0, lexp, 64);
3913   /* convert to external format */
3914
3915
3916   /* Multiply by 10**nexp.  If precision is 64 bits,
3917    * the maximum relative error incurred in forming 10**n
3918    * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3919    * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3920    * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3921    */
3922   lexp = yy[E];
3923   if (nexp == 0)
3924     {
3925       k = 0;
3926       goto expdon;
3927     }
3928   esign = 1;
3929   if (nexp < 0)
3930     {
3931       nexp = -nexp;
3932       esign = -1;
3933       if (nexp > 4096)
3934         {                       /* Punt.  Can't handle this without 2 divides. */
3935           emovi (etens[0], tt);
3936           lexp -= tt[E];
3937           k = edivm (tt, yy);
3938           lexp += EXONE;
3939           nexp -= 4096;
3940         }
3941     }
3942   p = &etens[NTEN][0];
3943   emov (eone, xt);
3944   exp = 1;
3945   do
3946     {
3947       if (exp & nexp)
3948         emul (p, xt, xt);
3949       p -= NE;
3950       exp = exp + exp;
3951     }
3952   while (exp <= MAXP);
3953
3954   emovi (xt, tt);
3955   if (esign < 0)
3956     {
3957       lexp -= tt[E];
3958       k = edivm (tt, yy);
3959       lexp += EXONE;
3960     }
3961   else
3962     {
3963       lexp += tt[E];
3964       k = emulm (tt, yy);
3965       lexp -= EXONE - 1;
3966     }
3967
3968  expdon:
3969
3970   /* Round and convert directly to the destination type */
3971   if (oprec == 53)
3972     lexp -= EXONE - 0x3ff;
3973   else if (oprec == 24)
3974     lexp -= EXONE - 0177;
3975 #ifdef DEC
3976   else if (oprec == 56)
3977     lexp -= EXONE - 0201;
3978 #endif
3979   rndprc = oprec;
3980   emdnorm (yy, k, 0, lexp, 64);
3981
3982  aexit:
3983
3984   rndprc = rndsav;
3985   yy[0] = nsign;
3986   switch (oprec)
3987     {
3988 #ifdef DEC
3989     case 56:
3990       todec (yy, y);            /* see etodec.c */
3991       break;
3992 #endif
3993     case 53:
3994       toe53 (yy, y);
3995       break;
3996     case 24:
3997       toe24 (yy, y);
3998       break;
3999     case 64:
4000       toe64 (yy, y);
4001       break;
4002     case NBITS:
4003       emovo (yy, y);
4004       break;
4005     }
4006 }
4007
4008
4009
4010 /* y = largest integer not greater than x
4011  * (truncated toward minus infinity)
4012  *
4013  * unsigned EMUSHORT x[NE], y[NE]
4014  *
4015  * efloor (x, y);
4016  */
4017 static unsigned EMUSHORT bmask[] =
4018 {
4019   0xffff,
4020   0xfffe,
4021   0xfffc,
4022   0xfff8,
4023   0xfff0,
4024   0xffe0,
4025   0xffc0,
4026   0xff80,
4027   0xff00,
4028   0xfe00,
4029   0xfc00,
4030   0xf800,
4031   0xf000,
4032   0xe000,
4033   0xc000,
4034   0x8000,
4035   0x0000,
4036 };
4037
4038 void
4039 efloor (x, y)
4040      unsigned EMUSHORT x[], y[];
4041 {
4042   register unsigned EMUSHORT *p;
4043   int e, expon, i;
4044   unsigned EMUSHORT f[NE];
4045
4046   emov (x, f);                  /* leave in external format */
4047   expon = (int) f[NE - 1];
4048   e = (expon & 0x7fff) - (EXONE - 1);
4049   if (e <= 0)
4050     {
4051       eclear (y);
4052       goto isitneg;
4053     }
4054   /* number of bits to clear out */
4055   e = NBITS - e;
4056   emov (f, y);
4057   if (e <= 0)
4058     return;
4059
4060   p = &y[0];
4061   while (e >= 16)
4062     {
4063       *p++ = 0;
4064       e -= 16;
4065     }
4066   /* clear the remaining bits */
4067   *p &= bmask[e];
4068   /* truncate negatives toward minus infinity */
4069  isitneg:
4070
4071   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4072     {
4073       for (i = 0; i < NE - 1; i++)
4074         {
4075           if (f[i] != y[i])
4076             {
4077               esub (eone, y, y);
4078               break;
4079             }
4080         }
4081     }
4082 }
4083
4084
4085 /* unsigned EMUSHORT x[], s[];
4086  * int *exp;
4087  *
4088  * efrexp (x, exp, s);
4089  *
4090  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
4091  * For example, 1.1 = 0.55 * 2**1
4092  * Handles denormalized numbers properly using long integer exp.
4093  */
4094 void
4095 efrexp (x, exp, s)
4096      unsigned EMUSHORT x[];
4097      int *exp;
4098      unsigned EMUSHORT s[];
4099 {
4100   unsigned EMUSHORT xi[NI];
4101   EMULONG li;
4102
4103   emovi (x, xi);
4104   li = (EMULONG) ((EMUSHORT) xi[1]);
4105
4106   if (li == 0)
4107     {
4108       li -= enormlz (xi);
4109     }
4110   xi[1] = 0x3ffe;
4111   emovo (xi, s);
4112   *exp = (int) (li - 0x3ffe);
4113 }
4114
4115
4116
4117 /* unsigned EMUSHORT x[], y[];
4118  * long pwr2;
4119  *
4120  * eldexp (x, pwr2, y);
4121  *
4122  * Returns y = x * 2**pwr2.
4123  */
4124 void
4125 eldexp (x, pwr2, y)
4126      unsigned EMUSHORT x[];
4127      int pwr2;
4128      unsigned EMUSHORT y[];
4129 {
4130   unsigned EMUSHORT xi[NI];
4131   EMULONG li;
4132   int i;
4133
4134   emovi (x, xi);
4135   li = xi[1];
4136   li += pwr2;
4137   i = 0;
4138   emdnorm (xi, i, i, li, 64);
4139   emovo (xi, y);
4140 }
4141
4142
4143 /* c = remainder after dividing b by a
4144  * Least significant integer quotient bits left in equot[].
4145  */
4146 void
4147 eremain (a, b, c)
4148      unsigned EMUSHORT a[], b[], c[];
4149 {
4150   unsigned EMUSHORT den[NI], num[NI];
4151
4152   if (ecmp (a, ezero) == 0)
4153     {
4154       mtherr ("eremain", SING);
4155       eclear (c);
4156       return;
4157     }
4158   emovi (a, den);
4159   emovi (b, num);
4160   eiremain (den, num);
4161   /* Sign of remainder = sign of quotient */
4162   if (a[0] == b[0])
4163     num[0] = 0;
4164   else
4165     num[0] = 0xffff;
4166   emovo (num, c);
4167 }
4168
4169 void
4170 eiremain (den, num)
4171      unsigned EMUSHORT den[], num[];
4172 {
4173   EMULONG ld, ln;
4174   unsigned EMUSHORT j;
4175
4176   ld = den[E];
4177   ld -= enormlz (den);
4178   ln = num[E];
4179   ln -= enormlz (num);
4180   ecleaz (equot);
4181   while (ln >= ld)
4182     {
4183       if (ecmpm (den, num) <= 0)
4184         {
4185           esubm (den, num);
4186           j = 1;
4187         }
4188       else
4189         {
4190           j = 0;
4191         }
4192       eshup1 (equot);
4193       equot[NI - 1] |= j;
4194       eshup1 (num);
4195       ln -= 1;
4196     }
4197   emdnorm (num, 0, 0, ln, 0);
4198 }
4199
4200 /*                                                      mtherr.c
4201  *
4202  *      Library common error handling routine
4203  *
4204  *
4205  *
4206  * SYNOPSIS:
4207  *
4208  * char *fctnam;
4209  * int code;
4210  * void mtherr ();
4211  *
4212  * mtherr (fctnam, code);
4213  *
4214  *
4215  *
4216  * DESCRIPTION:
4217  *
4218  * This routine may be called to report one of the following
4219  * error conditions (in the include file mconf.h).
4220  *
4221  *   Mnemonic        Value          Significance
4222  *
4223  *    DOMAIN            1       argument domain error
4224  *    SING              2       function singularity
4225  *    OVERFLOW          3       overflow range error
4226  *    UNDERFLOW         4       underflow range error
4227  *    TLOSS             5       total loss of precision
4228  *    PLOSS             6       partial loss of precision
4229  *    EDOM             33       Unix domain error code
4230  *    ERANGE           34       Unix range error code
4231  *
4232  * The default version of the file prints the function name,
4233  * passed to it by the pointer fctnam, followed by the
4234  * error condition.  The display is directed to the standard
4235  * output device.  The routine then returns to the calling
4236  * program.  Users may wish to modify the program to abort by
4237  * calling exit under severe error conditions such as domain
4238  * errors.
4239  *
4240  * Since all error conditions pass control to this function,
4241  * the display may be easily changed, eliminated, or directed
4242  * to an error logging device.
4243  *
4244  * SEE ALSO:
4245  *
4246  * mconf.h
4247  *
4248  */
4249 \f
4250 /*
4251 Cephes Math Library Release 2.0:  April, 1987
4252 Copyright 1984, 1987 by Stephen L. Moshier
4253 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4254 */
4255
4256 /* include "mconf.h" */
4257
4258 /* Notice: the order of appearance of the following
4259  * messages is bound to the error codes defined
4260  * in mconf.h.
4261  */
4262 static char *ermsg[7] =
4263 {
4264   "unknown",                    /* error code 0 */
4265   "domain",                     /* error code 1 */
4266   "singularity",                /* et seq.      */
4267   "overflow",
4268   "underflow",
4269   "total loss of precision",
4270   "partial loss of precision"
4271 };
4272
4273 int merror = 0;
4274 extern int merror;
4275
4276 void
4277 mtherr (name, code)
4278      char *name;
4279      int code;
4280 {
4281   char errstr[80];
4282
4283   /* Display string passed by calling program,
4284    * which is supposed to be the name of the
4285    * function in which the error occurred.
4286    */
4287
4288   /* Display error message defined
4289    * by the code argument.
4290    */
4291   if ((code <= 0) || (code >= 6))
4292     code = 0;
4293   sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
4294   if (extra_warnings)
4295     warning (errstr);
4296   /* Set global error message word */
4297   merror = code + 1;
4298
4299   /* Return to calling
4300    * program
4301    */
4302 }
4303
4304 /* Here is etodec.c .
4305  *
4306  */
4307
4308 /*
4309 ;       convert DEC double precision to e type
4310 ;       double d;
4311 ;       EMUSHORT e[NE];
4312 ;       dectoe (&d, e);
4313 */
4314 void
4315 dectoe (d, e)
4316      unsigned EMUSHORT *d;
4317      unsigned EMUSHORT *e;
4318 {
4319   unsigned EMUSHORT y[NI];
4320   register unsigned EMUSHORT r, *p;
4321
4322   ecleaz (y);                   /* start with a zero */
4323   p = y;                        /* point to our number */
4324   r = *d;                       /* get DEC exponent word */
4325   if (*d & (unsigned int) 0x8000)
4326     *p = 0xffff;                /* fill in our sign */
4327   ++p;                          /* bump pointer to our exponent word */
4328   r &= 0x7fff;                  /* strip the sign bit */
4329   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
4330     goto done;
4331
4332
4333   r >>= 7;                      /* shift exponent word down 7 bits */
4334   r += EXONE - 0201;            /* subtract DEC exponent offset */
4335   /* add our e type exponent offset */
4336   *p++ = r;                     /* to form our exponent */
4337
4338   r = *d++;                     /* now do the high order mantissa */
4339   r &= 0177;                    /* strip off the DEC exponent and sign bits */
4340   r |= 0200;                    /* the DEC understood high order mantissa bit */
4341   *p++ = r;                     /* put result in our high guard word */
4342
4343   *p++ = *d++;                  /* fill in the rest of our mantissa */
4344   *p++ = *d++;
4345   *p = *d;
4346
4347   eshdn8 (y);                   /* shift our mantissa down 8 bits */
4348  done:
4349   emovo (y, e);
4350 }
4351
4352
4353
4354 /*
4355 ;       convert e type to DEC double precision
4356 ;       double d;
4357 ;       EMUSHORT e[NE];
4358 ;       etodec (e, &d);
4359 */
4360 #if 0
4361 static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4362
4363 void
4364 etodec (x, d)
4365      unsigned EMUSHORT *x, *d;
4366 {
4367   unsigned EMUSHORT xi[NI];
4368   register unsigned EMUSHORT r;
4369   int i, j;
4370
4371   emovi (x, xi);
4372   *d = 0;
4373   if (xi[0] != 0)
4374     *d = 0100000;
4375   r = xi[E];
4376   if (r < (EXONE - 128))
4377     goto zout;
4378   i = xi[M + 4];
4379   if ((i & 0200) != 0)
4380     {
4381       if ((i & 0377) == 0200)
4382         {
4383           if ((i & 0400) != 0)
4384             {
4385               /* check all less significant bits */
4386               for (j = M + 5; j < NI; j++)
4387                 {
4388                   if (xi[j] != 0)
4389                     goto yesrnd;
4390                 }
4391             }
4392           goto nornd;
4393         }
4394     yesrnd:
4395       eaddm (decbit, xi);
4396       r -= enormlz (xi);
4397     }
4398
4399  nornd:
4400
4401   r -= EXONE;
4402   r += 0201;
4403   if (r < 0)
4404     {
4405     zout:
4406       *d++ = 0;
4407       *d++ = 0;
4408       *d++ = 0;
4409       *d++ = 0;
4410       return;
4411     }
4412   if (r >= 0377)
4413     {
4414       *d++ = 077777;
4415       *d++ = -1;
4416       *d++ = -1;
4417       *d++ = -1;
4418       return;
4419     }
4420   r &= 0377;
4421   r <<= 7;
4422   eshup8 (xi);
4423   xi[M] &= 0177;
4424   r |= xi[M];
4425   *d++ |= r;
4426   *d++ = xi[M + 1];
4427   *d++ = xi[M + 2];
4428   *d++ = xi[M + 3];
4429 }
4430
4431 #else
4432
4433 void
4434 etodec (x, d)
4435      unsigned EMUSHORT *x, *d;
4436 {
4437   unsigned EMUSHORT xi[NI];
4438   EMULONG exp;
4439   int rndsav;
4440
4441   emovi (x, xi);
4442   exp = (EMULONG) xi[E] - (EXONE - 0201);       /* adjust exponent for offsets */
4443 /* round off to nearest or even */
4444   rndsav = rndprc;
4445   rndprc = 56;
4446   emdnorm (xi, 0, 0, exp, 64);
4447   rndprc = rndsav;
4448   todec (xi, d);
4449 }
4450
4451 void
4452 todec (x, y)
4453      unsigned EMUSHORT *x, *y;
4454 {
4455   unsigned EMUSHORT i;
4456   unsigned EMUSHORT *p;
4457
4458   p = x;
4459   *y = 0;
4460   if (*p++)
4461     *y = 0100000;
4462   i = *p++;
4463   if (i == 0)
4464     {
4465       *y++ = 0;
4466       *y++ = 0;
4467       *y++ = 0;
4468       *y++ = 0;
4469       return;
4470     }
4471   if (i > 0377)
4472     {
4473       *y++ |= 077777;
4474       *y++ = 0xffff;
4475       *y++ = 0xffff;
4476       *y++ = 0xffff;
4477 #ifdef ERANGE
4478       errno = ERANGE;
4479 #endif
4480       return;
4481     }
4482   i &= 0377;
4483   i <<= 7;
4484   eshup8 (x);
4485   x[M] &= 0177;
4486   i |= x[M];
4487   *y++ |= i;
4488   *y++ = x[M + 1];
4489   *y++ = x[M + 2];
4490   *y++ = x[M + 3];
4491 }
4492
4493 #endif /* not 0 */
4494
4495 #endif /* EMU_NON_COMPILE not defined */