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