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