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