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    Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998,
   4    1999, 2000, 2002 Free Software Foundation, Inc.
   5    Contributed by Stephen L. Moshier (moshier@world.std.com).
   6
   7 This file is part of GCC.
   8
   9 GCC is free software; you can redistribute it and/or modify it under
  10 the terms of the GNU General Public License as published by the Free
  11 Software Foundation; either version 2, or (at your option) any later
  12 version.
  13
  14 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  15 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  17 for more details.
  18
  19 You should have received a copy of the GNU General Public License
  20 along with GCC; see the file COPYING.  If not, write to the Free
  21 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  22 02111-1307, USA.  */
  23
  24 #include "config.h"
  25 #include "system.h"
  26 #include "real.h"
  27 #include "tree.h"
  28 #include "toplev.h"
  29 #include "tm_p.h"
  30
  31 /* To enable support of XFmode extended real floating point, define
  32 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  33
  34 Machine files (tm.h etc) must not contain any code
  35 that tries to use host floating point arithmetic to convert
  36 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  37 etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  38 be intelligible to the host computer's native arithmetic.
  39
  40 The first part of this file interfaces gcc to a floating point
  41 arithmetic suite that was not written with gcc in mind.  Avoid
  42 changing the low-level arithmetic routines unless you have suitable
  43 test programs available.  A special version of the PARANOIA floating
  44 point arithmetic tester, modified for this purpose, can be found on
  45 usc.edu: /pub/C-numanal/ieeetest.zoo.  Other tests, and libraries of
  46 XFmode and TFmode transcendental functions, can be obtained by ftp from
  47 netlib.att.com: netlib/cephes.  */
  48 \f
  49 /* Type of computer arithmetic.
  50    Only one of DEC, IBM, IEEE, C4X, or UNK should get defined.
  51
  52    `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
  53    to big-endian IEEE floating-point data structure.  This definition
  54    should work in SFmode `float' type and DFmode `double' type on
  55    virtually all big-endian IEEE machines.  If LONG_DOUBLE_TYPE_SIZE
  56    has been defined to be 96, then IEEE also invokes the particular
  57    XFmode (`long double' type) data structure used by the Motorola
  58    680x0 series processors.
  59
  60    `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
  61    little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
  62    has been defined to be 96, then IEEE also invokes the particular
  63    XFmode `long double' data structure used by the Intel 80x86 series
  64    processors.
  65
  66    `DEC' refers specifically to the Digital Equipment Corp PDP-11
  67    and VAX floating point data structure.  This model currently
  68    supports no type wider than DFmode.
  69
  70    `IBM' refers specifically to the IBM System/370 and compatible
  71    floating point data structure.  This model currently supports
  72    no type wider than DFmode.  The IBM conversions were contributed by
  73    frank@atom.ansto.gov.au (Frank Crawford).
  74
  75    `C4X' refers specifically to the floating point format used on
  76    Texas Instruments TMS320C3x and TMS320C4x digital signal
  77    processors.  This supports QFmode (32-bit float, double) and HFmode
  78    (40-bit long double) where BITS_PER_BYTE is 32. Unlike IEEE
  79    floats, C4x floats are not rounded to be even. The C4x conversions
  80    were contributed by m.hayes@elec.canterbury.ac.nz (Michael Hayes) and
  81    Haj.Ten.Brugge@net.HCC.nl (Herman ten Brugge).
  82
  83    If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
  84    then `long double' and `double' are both implemented, but they
  85    both mean DFmode.
  86
  87    The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
  88    and may deactivate XFmode since `long double' is used to refer
  89    to both modes.  Defining INTEL_EXTENDED_IEEE_FORMAT to non-zero
  90    at the same time enables 80387-style 80-bit floats in a 128-bit
  91    padded image, as seen on IA-64.
  92
  93    The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
  94    contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
  95    separate the floating point unit's endian-ness from that of
  96    the integer addressing.  This permits one to define a big-endian
  97    FPU on a little-endian machine (e.g., ARM).  An extension to
  98    BYTES_BIG_ENDIAN may be required for some machines in the future.
  99    These optional macros may be defined in tm.h.  In real.h, they
 100    default to WORDS_BIG_ENDIAN, etc., so there is no need to define
 101    them for any normal host or target machine on which the floats
 102    and the integers have the same endian-ness.  */
 103
 104
 105 /* The following converts gcc macros into the ones used by this file.  */
 106
 107 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 108 /* PDP-11, Pro350, VAX: */
 109 #define DEC 1
 110 #else /* it's not VAX */
 111 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
 112 /* IBM System/370 style */
 113 #define IBM 1
 114 #else /* it's also not an IBM */
 115 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
 116 /* TMS320C3x/C4x style */
 117 #define C4X 1
 118 #else /* it's also not a C4X */
 119 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 120 #define IEEE
 121 #else /* it's not IEEE either */
 122 /* UNKnown arithmetic.  We don't support this and can't go on.  */
 123 unknown arithmetic type
 124 #define UNK 1
 125 #endif /* not IEEE */
 126 #endif /* not C4X */
 127 #endif /* not IBM */
 128 #endif /* not VAX */
 129
 130 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
 131
 132 /* Make sure that the endianness is correct for IBM and DEC. */
 133 #if defined(DEC)
 134 #undef LARGEST_EXPONENT_IS_NORMAL
 135 #define LARGEST_EXPONENT_IS_NORMAL(x) 1
 136 #undef REAL_WORDS_BIG_ENDIAN
 137 /* Strangely enough, DEC float most closely resembles big endian IEEE */
 138 #define REAL_WORDS_BIG_ENDIAN 1
 139 /* ... but the halfwords are reversed from IEEE big endian. */
 140 #ifndef VAX_HALFWORD_ORDER
 141 #define VAX_HALFWORD_ORDER 1
 142 #endif
 143 #else
 144 #if defined(IBM) && !REAL_WORDS_BIG_ENDIAN
 145   #error "Little-endian representations are not supported for IBM."
 146 #endif
 147 #endif
 148
 149 #if defined(DEC) && !defined (TARGET_G_FLOAT)
 150 #define TARGET_G_FLOAT 0
 151 #endif
 152
 153 #ifndef VAX_HALFWORD_ORDER
 154 #define VAX_HALFWORD_ORDER 0
 155 #endif
 156
 157 /* Define INFINITY for support of infinity.
 158    Define NANS for support of Not-a-Number's (NaN's).  */
 159 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
 160 #define INFINITY
 161 #define NANS
 162 #endif
 163
 164 /* Support of NaNs requires support of infinity.  */
 165 #ifdef NANS
 166 #ifndef INFINITY
 167 #define INFINITY
 168 #endif
 169 #endif
 170 \f
 171 /* Find a host integer type that is at least 16 bits wide,
 172    and another type at least twice whatever that size is.  */
 173
 174 #if HOST_BITS_PER_CHAR >= 16
 175 #define EMUSHORT char
 176 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
 177 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
 178 #else
 179 #if HOST_BITS_PER_SHORT >= 16
 180 #define EMUSHORT short
 181 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
 182 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
 183 #else
 184 #if HOST_BITS_PER_INT >= 16
 185 #define EMUSHORT int
 186 #define EMUSHORT_SIZE HOST_BITS_PER_INT
 187 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
 188 #else
 189 #if HOST_BITS_PER_LONG >= 16
 190 #define EMUSHORT long
 191 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 192 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 193 #else
 194   #error "You will have to modify this program to have a smaller unit size."
 195 #endif
 196 #endif
 197 #endif
 198 #endif
 199
 200 /* If no 16-bit type has been found and the compiler is GCC, try HImode.  */
 201 #if defined(__GNUC__) && EMUSHORT_SIZE != 16
 202 typedef int HItype __attribute__ ((mode (HI)));
 203 typedef unsigned int UHItype __attribute__ ((mode (HI)));
 204 #undef EMUSHORT
 205 #undef EMUSHORT_SIZE
 206 #undef EMULONG_SIZE
 207 #define EMUSHORT HItype
 208 #define UEMUSHORT UHItype
 209 #define EMUSHORT_SIZE 16
 210 #define EMULONG_SIZE 32
 211 #else
 212 #define UEMUSHORT unsigned EMUSHORT
 213 #endif
 214
 215 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
 216 #define EMULONG short
 217 #else
 218 #if HOST_BITS_PER_INT >= EMULONG_SIZE
 219 #define EMULONG int
 220 #else
 221 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
 222 #define EMULONG long
 223 #else
 224 #if HOST_BITS_PER_LONGLONG >= EMULONG_SIZE
 225 #define EMULONG long long int
 226 #else
 227   #error "You will have to modify this program to have a smaller unit size."
 228 #endif
 229 #endif
 230 #endif
 231 #endif
 232
 233 #if EMUSHORT_SIZE != 16
 234   #error "The host interface doesn't work if no 16-bit size exists."
 235 #endif
 236
 237 /* Calculate the size of the generic "e" type.  This always has
 238    identical in-memory size to REAL_VALUE_TYPE.  The sizes are supposed
 239    to be the same as well, but when REAL_VALUE_TYPE_SIZE is not evenly
 240    divisible by HOST_BITS_PER_WIDE_INT we have some padding in
 241    REAL_VALUE_TYPE.
 242    There are only two supported sizes: ten and six 16-bit words (160
 243    or 96 bits).  */
 244
 245 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && !INTEL_EXTENDED_IEEE_FORMAT
 246 /* TFmode */
 247 # define NE 10
 248 # define MAXDECEXP 4932
 249 # define MINDECEXP -4977
 250 #else
 251 # define NE 6
 252 # define MAXDECEXP 4932
 253 # define MINDECEXP -4956
 254 #endif
 255
 256 /* Fail compilation if 2*NE is not the appropriate size.
 257    If HOST_BITS_PER_WIDE_INT is 64, we're going to have padding
 258    at the end of the array, because neither 96 nor 160 is
 259    evenly divisible by 64.  */
 260 struct compile_test_dummy {
 261   char twice_NE_must_equal_sizeof_REAL_VALUE_TYPE
 262   [(sizeof (REAL_VALUE_TYPE) >= 2*NE) ? 1 : -1];
 263 };
 264
 265 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
 266    In GET_REAL and PUT_REAL, r and e are pointers.
 267    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
 268    in memory, with no holes.  */
 269 #define GET_REAL(r, e)  memcpy ((e), (r), 2*NE)
 270 #define PUT_REAL(e, r)                                          \
 271   do {                                                          \
 272     memcpy (r, e, 2*NE);                                        \
 273     if (2*NE < sizeof (*r))                                     \
 274       memset ((char *) (r) + 2*NE, 0, sizeof (*r) - 2*NE);      \
 275   } while (0)
 276
 277 /* Number of 16 bit words in internal format */
 278 #define NI (NE+3)
 279
 280 /* Array offset to exponent */
 281 #define E 1
 282
 283 /* Array offset to high guard word */
 284 #define M 2
 285
 286 /* Number of bits of precision */
 287 #define NBITS ((NI-4)*16)
 288
 289 /* Maximum number of decimal digits in ASCII conversion
 290  * = NBITS*log10(2)
 291  */
 292 #define NDEC (NBITS*8/27)
 293
 294 /* The exponent of 1.0 */
 295 #define EXONE (0x3fff)
 296
 297 #if defined(HOST_EBCDIC)
 298 /* bit 8 is significant in EBCDIC */
 299 #define CHARMASK 0xff
 300 #else
 301 #define CHARMASK 0x7f
 302 #endif
 303
 304 /* Information about the various IEEE precisions. At the moment, we only
 305    support exponents of 15 bits or less.  */
 306 struct ieee_format
 307 {
 308   /* Precision.  */
 309   int precision;
 310
 311   /* Size of the exponent in bits.  */
 312   int expbits;
 313
 314   /* Overall size of the value in bits.  */
 315   int bits;
 316
 317   /* Mode used for representing the value.  */
 318   enum machine_mode mode;
 319
 320   /* Exponent adjustment for offsets.  */
 321   EMULONG adjustment;
 322 };
 323
 324 /* IEEE float (24 bits).  */
 325 static const struct ieee_format ieee_24 =
 326 {
 327   24,
 328   8,
 329   32,
 330   SFmode,
 331   EXONE - 0x7f
 332 };
 333
 334 /* IEEE double (53 bits).  */
 335 static const struct ieee_format ieee_53 =
 336 {
 337   53,
 338   11,
 339   64,
 340   DFmode,
 341   EXONE - 0x3ff
 342 };
 343
 344 /* IEEE extended double (64 bits).  */
 345 static const struct ieee_format ieee_64 =
 346 {
 347   64,
 348   15,
 349   80,
 350   XFmode,
 351   0
 352 };
 353
 354 /* IEEE long double (113 bits).  */
 355 static const struct ieee_format ieee_113 =
 356 {
 357   113,
 358   15,
 359   128,
 360   TFmode,
 361   0
 362 };
 363
 364 /* DEC F float (24 bits).  */
 365 static const struct ieee_format dec_f =
 366 {
 367   24,
 368   8,
 369   32,
 370   SFmode,
 371   EXONE - 0201
 372 };
 373
 374 /* DEC D float (56 bits).  */
 375 static const struct ieee_format dec_d =
 376 {
 377   56,
 378   8,
 379   64,
 380   DFmode,
 381   EXONE - 0201
 382 };
 383
 384 /* DEC G float (53 bits).  */
 385 static const struct ieee_format dec_g =
 386 {
 387   53,
 388   11,
 389   64,
 390   DFmode,
 391   EXONE - 1025
 392 };
 393
 394 /* DEC H float (113 bits).  (not yet used) */
 395 static const struct ieee_format dec_h =
 396 {
 397   113,
 398   15,
 399   128,
 400   TFmode,
 401   EXONE - 16385
 402 };
 403
 404 extern int extra_warnings;
 405 extern const UEMUSHORT ezero[NE], ehalf[NE], eone[NE], etwo[NE];
 406 extern const UEMUSHORT elog2[NE], esqrt2[NE];
 407
 408 static void endian      PARAMS ((const UEMUSHORT *, long *,
 409                                enum machine_mode));
 410 static void eclear      PARAMS ((UEMUSHORT *));
 411 static void emov        PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 412 #if 0
 413 static void eabs        PARAMS ((UEMUSHORT *));
 414 #endif
 415 static void eneg        PARAMS ((UEMUSHORT *));
 416 static int eisneg       PARAMS ((const UEMUSHORT *));
 417 static int eisinf       PARAMS ((const UEMUSHORT *));
 418 static int eisnan       PARAMS ((const UEMUSHORT *));
 419 static void einfin      PARAMS ((UEMUSHORT *));
 420 #ifdef NANS
 421 static void enan        PARAMS ((UEMUSHORT *, int));
 422 static void einan       PARAMS ((UEMUSHORT *));
 423 static int eiisnan      PARAMS ((const UEMUSHORT *));
 424 static void make_nan    PARAMS ((UEMUSHORT *, int, enum machine_mode));
 425 #endif
 426 static int eiisneg      PARAMS ((const UEMUSHORT *));
 427 static void saturate    PARAMS ((UEMUSHORT *, int, int, int));
 428 static void emovi       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 429 static void emovo       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 430 static void ecleaz      PARAMS ((UEMUSHORT *));
 431 static void ecleazs     PARAMS ((UEMUSHORT *));
 432 static void emovz       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 433 #if 0
 434 static void eiinfin     PARAMS ((UEMUSHORT *));
 435 #endif
 436 #ifdef INFINITY
 437 static int eiisinf      PARAMS ((const UEMUSHORT *));
 438 #endif
 439 static int ecmpm        PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
 440 static void eshdn1      PARAMS ((UEMUSHORT *));
 441 static void eshup1      PARAMS ((UEMUSHORT *));
 442 static void eshdn8      PARAMS ((UEMUSHORT *));
 443 static void eshup8      PARAMS ((UEMUSHORT *));
 444 static void eshup6      PARAMS ((UEMUSHORT *));
 445 static void eshdn6      PARAMS ((UEMUSHORT *));
 446 static void eaddm       PARAMS ((const UEMUSHORT *, UEMUSHORT *));\f
 447 static void esubm       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 448 static void m16m        PARAMS ((unsigned int, const UEMUSHORT *, UEMUSHORT *));
 449 static int edivm        PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 450 static int emulm        PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 451 static void emdnorm     PARAMS ((UEMUSHORT *, int, int, EMULONG, int));
 452 static void esub        PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 453                                  UEMUSHORT *));
 454 static void eadd        PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 455                                  UEMUSHORT *));
 456 static void eadd1       PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 457                                  UEMUSHORT *));
 458 static void ediv        PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 459                                  UEMUSHORT *));
 460 static void emul        PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 461                                  UEMUSHORT *));
 462 static void e53toe      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 463 static void e64toe      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 464 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 465 static void e113toe     PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 466 #endif
 467 static void e24toe      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 468 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 469 static void etoe113     PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 470 static void toe113      PARAMS ((UEMUSHORT *, UEMUSHORT *));
 471 #endif
 472 static void etoe64      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 473 static void toe64       PARAMS ((UEMUSHORT *, UEMUSHORT *));
 474 static void etoe53      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 475 static void toe53       PARAMS ((UEMUSHORT *, UEMUSHORT *));
 476 static void etoe24      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 477 static void toe24       PARAMS ((UEMUSHORT *, UEMUSHORT *));
 478 static void ieeetoe     PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 479                                  const struct ieee_format *));
 480 static void etoieee     PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 481                                  const struct ieee_format *));
 482 static void toieee      PARAMS ((UEMUSHORT *, UEMUSHORT *,
 483                                  const struct ieee_format *));
 484 static int ecmp         PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
 485 #if 0
 486 static void eround      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 487 #endif
 488 static void ltoe        PARAMS ((const HOST_WIDE_INT *, UEMUSHORT *));
 489 static void ultoe       PARAMS ((const unsigned HOST_WIDE_INT *, UEMUSHORT *));
 490 static void eifrac      PARAMS ((const UEMUSHORT *, HOST_WIDE_INT *,
 491                                  UEMUSHORT *));
 492 static void euifrac     PARAMS ((const UEMUSHORT *, unsigned HOST_WIDE_INT *,
 493                                  UEMUSHORT *));
 494 static int eshift       PARAMS ((UEMUSHORT *, int));
 495 static int enormlz      PARAMS ((UEMUSHORT *));
 496 #if 0
 497 static void e24toasc    PARAMS ((const UEMUSHORT *, char *, int));
 498 static void e53toasc    PARAMS ((const UEMUSHORT *, char *, int));
 499 static void e64toasc    PARAMS ((const UEMUSHORT *, char *, int));
 500 static void e113toasc   PARAMS ((const UEMUSHORT *, char *, int));
 501 #endif /* 0 */
 502 static void etoasc      PARAMS ((const UEMUSHORT *, char *, int));
 503 static void asctoe24    PARAMS ((const char *, UEMUSHORT *));
 504 static void asctoe53    PARAMS ((const char *, UEMUSHORT *));
 505 static void asctoe64    PARAMS ((const char *, UEMUSHORT *));
 506 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 507 static void asctoe113   PARAMS ((const char *, UEMUSHORT *));
 508 #endif
 509 static void asctoe      PARAMS ((const char *, UEMUSHORT *));
 510 static void asctoeg     PARAMS ((const char *, UEMUSHORT *, int));
 511 static void efloor      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 512 #if 0
 513 static void efrexp      PARAMS ((const UEMUSHORT *, int *,
 514                                  UEMUSHORT *));
 515 #endif
 516 static void eldexp      PARAMS ((const UEMUSHORT *, int, UEMUSHORT *));
 517 #if 0
 518 static void eremain     PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
 519                                  UEMUSHORT *));
 520 #endif
 521 static void eiremain    PARAMS ((UEMUSHORT *, UEMUSHORT *));
 522 static void mtherr      PARAMS ((const char *, int));
 523 #ifdef DEC
 524 static void dectoe      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 525 static void etodec      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 526 static void todec       PARAMS ((UEMUSHORT *, UEMUSHORT *));
 527 #endif
 528 #ifdef IBM
 529 static void ibmtoe      PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 530                                  enum machine_mode));
 531 static void etoibm      PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 532                                  enum machine_mode));
 533 static void toibm       PARAMS ((UEMUSHORT *, UEMUSHORT *,
 534                                  enum machine_mode));
 535 #endif
 536 #ifdef C4X
 537 static void c4xtoe      PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 538                                  enum machine_mode));
 539 static void etoc4x      PARAMS ((const UEMUSHORT *, UEMUSHORT *,
 540                                  enum machine_mode));
 541 static void toc4x       PARAMS ((UEMUSHORT *, UEMUSHORT *,
 542                                  enum machine_mode));
 543 #endif
 544 #if 0
 545 static void uditoe      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 546 static void ditoe       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 547 static void etoudi      PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 548 static void etodi       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 549 static void esqrt       PARAMS ((const UEMUSHORT *, UEMUSHORT *));
 550 #endif
 551 \f
 552 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
 553    swapping ends if required, into output array of longs.  The
 554    result is normally passed to fprintf by the ASM_OUTPUT_ macros.  */
 555
 556 static void
 557 endian (e, x, mode)
 558      const UEMUSHORT e[];
 559      long x[];
 560      enum machine_mode mode;
 561 {
 562   unsigned long th, t;
 563
 564   if (REAL_WORDS_BIG_ENDIAN && !VAX_HALFWORD_ORDER)
 565     {
 566       switch (mode)
 567         {
 568         case TFmode:
 569 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 570           /* Swap halfwords in the fourth long.  */
 571           th = (unsigned long) e[6] & 0xffff;
 572           t = (unsigned long) e[7] & 0xffff;
 573           t |= th << 16;
 574           x[3] = (long) t;
 575 #else
 576           x[3] = 0;
 577 #endif
 578           /* FALLTHRU */
 579
 580         case XFmode:
 581           /* Swap halfwords in the third long.  */
 582           th = (unsigned long) e[4] & 0xffff;
 583           t = (unsigned long) e[5] & 0xffff;
 584           t |= th << 16;
 585           x[2] = (long) t;
 586           /* FALLTHRU */
 587
 588         case DFmode:
 589           /* Swap halfwords in the second word.  */
 590           th = (unsigned long) e[2] & 0xffff;
 591           t = (unsigned long) e[3] & 0xffff;
 592           t |= th << 16;
 593           x[1] = (long) t;
 594           /* FALLTHRU */
 595
 596         case SFmode:
 597         case HFmode:
 598           /* Swap halfwords in the first word.  */
 599           th = (unsigned long) e[0] & 0xffff;
 600           t = (unsigned long) e[1] & 0xffff;
 601           t |= th << 16;
 602           x[0] = (long) t;
 603           break;
 604
 605         default:
 606           abort ();
 607         }
 608     }
 609   else
 610     {
 611       /* Pack the output array without swapping.  */
 612
 613       switch (mode)
 614         {
 615         case TFmode:
 616 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 617           /* Pack the fourth long.  */
 618           th = (unsigned long) e[7] & 0xffff;
 619           t = (unsigned long) e[6] & 0xffff;
 620           t |= th << 16;
 621           x[3] = (long) t;
 622 #else
 623           x[3] = 0;
 624 #endif
 625           /* FALLTHRU */
 626
 627         case XFmode:
 628           /* Pack the third long.
 629              Each element of the input REAL_VALUE_TYPE array has 16 useful bits
 630              in it.  */
 631           th = (unsigned long) e[5] & 0xffff;
 632           t = (unsigned long) e[4] & 0xffff;
 633           t |= th << 16;
 634           x[2] = (long) t;
 635           /* FALLTHRU */
 636
 637         case DFmode:
 638           /* Pack the second long */
 639           th = (unsigned long) e[3] & 0xffff;
 640           t = (unsigned long) e[2] & 0xffff;
 641           t |= th << 16;
 642           x[1] = (long) t;
 643           /* FALLTHRU */
 644
 645         case SFmode:
 646         case HFmode:
 647           /* Pack the first long */
 648           th = (unsigned long) e[1] & 0xffff;
 649           t = (unsigned long) e[0] & 0xffff;
 650           t |= th << 16;
 651           x[0] = (long) t;
 652           break;
 653
 654         default:
 655           abort ();
 656         }
 657     }
 658 }
 659
 660
 661 /* This is the implementation of the REAL_ARITHMETIC macro.  */
 662
 663 void
 664 earith (value, icode, r1, r2)
 665      REAL_VALUE_TYPE *value;
 666      int icode;
 667      REAL_VALUE_TYPE *r1;
 668      REAL_VALUE_TYPE *r2;
 669 {
 670   UEMUSHORT d1[NE], d2[NE], v[NE];
 671   enum tree_code code;
 672
 673   GET_REAL (r1, d1);
 674   GET_REAL (r2, d2);
 675 #ifdef NANS
 676 /*  Return NaN input back to the caller.  */
 677   if (eisnan (d1))
 678     {
 679       PUT_REAL (d1, value);
 680       return;
 681     }
 682   if (eisnan (d2))
 683     {
 684       PUT_REAL (d2, value);
 685       return;
 686     }
 687 #endif
 688   code = (enum tree_code) icode;
 689   switch (code)
 690     {
 691     case PLUS_EXPR:
 692       eadd (d2, d1, v);
 693       break;
 694
 695     case MINUS_EXPR:
 696       esub (d2, d1, v);         /* d1 - d2 */
 697       break;
 698
 699     case MULT_EXPR:
 700       emul (d2, d1, v);
 701       break;
 702
 703     case RDIV_EXPR:
 704 #ifndef INFINITY
 705       if (ecmp (d2, ezero) == 0)
 706         abort ();
 707 #endif
 708       ediv (d2, d1, v); /* d1/d2 */
 709       break;
 710
 711     case MIN_EXPR:              /* min (d1,d2) */
 712       if (ecmp (d1, d2) < 0)
 713         emov (d1, v);
 714       else
 715         emov (d2, v);
 716       break;
 717
 718     case MAX_EXPR:              /* max (d1,d2) */
 719       if (ecmp (d1, d2) > 0)
 720         emov (d1, v);
 721       else
 722         emov (d2, v);
 723       break;
 724     default:
 725       emov (ezero, v);
 726       break;
 727     }
 728 PUT_REAL (v, value);
 729 }
 730
 731
 732 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
 733    implements REAL_VALUE_RNDZINT (x) (etrunci (x)).  */
 734
 735 REAL_VALUE_TYPE
 736 etrunci (x)
 737      REAL_VALUE_TYPE x;
 738 {
 739   UEMUSHORT f[NE], g[NE];
 740   REAL_VALUE_TYPE r;
 741   HOST_WIDE_INT l;
 742
 743   GET_REAL (&x, g);
 744 #ifdef NANS
 745   if (eisnan (g))
 746     return (x);
 747 #endif
 748   eifrac (g, &l, f);
 749   ltoe (&l, g);
 750   PUT_REAL (g, &r);
 751   return (r);
 752 }
 753
 754
 755 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
 756    implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)).  */
 757
 758 REAL_VALUE_TYPE
 759 etruncui (x)
 760      REAL_VALUE_TYPE x;
 761 {
 762   UEMUSHORT f[NE], g[NE];
 763   REAL_VALUE_TYPE r;
 764   unsigned HOST_WIDE_INT l;
 765
 766   GET_REAL (&x, g);
 767 #ifdef NANS
 768   if (eisnan (g))
 769     return (x);
 770 #endif
 771   euifrac (g, &l, f);
 772   ultoe (&l, g);
 773   PUT_REAL (g, &r);
 774   return (r);
 775 }
 776
 777
 778 /* This is the REAL_VALUE_ATOF function.  It converts a decimal or hexadecimal
 779    string to binary, rounding off as indicated by the machine_mode argument.
 780    Then it promotes the rounded value to REAL_VALUE_TYPE.  */
 781
 782 REAL_VALUE_TYPE
 783 ereal_atof (s, t)
 784      const char *s;
 785      enum machine_mode t;
 786 {
 787   UEMUSHORT tem[NE], e[NE];
 788   REAL_VALUE_TYPE r;
 789
 790   switch (t)
 791     {
 792 #ifdef C4X
 793     case QFmode:
 794     case HFmode:
 795       asctoe53 (s, tem);
 796       e53toe (tem, e);
 797       break;
 798 #else
 799     case HFmode:
 800 #endif
 801
 802     case SFmode:
 803       asctoe24 (s, tem);
 804       e24toe (tem, e);
 805       break;
 806
 807     case DFmode:
 808       asctoe53 (s, tem);
 809       e53toe (tem, e);
 810       break;
 811
 812     case TFmode:
 813 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 814       asctoe113 (s, tem);
 815       e113toe (tem, e);
 816       break;
 817 #endif
 818       /* FALLTHRU */
 819
 820     case XFmode:
 821       asctoe64 (s, tem);
 822       e64toe (tem, e);
 823       break;
 824
 825     default:
 826       asctoe (s, e);
 827     }
 828   PUT_REAL (e, &r);
 829   return (r);
 830 }
 831
 832
 833 /* Expansion of REAL_NEGATE.  */
 834
 835 REAL_VALUE_TYPE
 836 ereal_negate (x)
 837      REAL_VALUE_TYPE x;
 838 {
 839   UEMUSHORT e[NE];
 840   REAL_VALUE_TYPE r;
 841
 842   GET_REAL (&x, e);
 843   eneg (e);
 844   PUT_REAL (e, &r);
 845   return (r);
 846 }
 847
 848
 849 /* Round real toward zero to HOST_WIDE_INT;
 850    implements REAL_VALUE_FIX (x).  */
 851
 852 HOST_WIDE_INT
 853 efixi (x)
 854      REAL_VALUE_TYPE x;
 855 {
 856   UEMUSHORT f[NE], g[NE];
 857   HOST_WIDE_INT l;
 858
 859   GET_REAL (&x, f);
 860 #ifdef NANS
 861   if (eisnan (f))
 862     {
 863       warning ("conversion from NaN to int");
 864       return (-1);
 865     }
 866 #endif
 867   eifrac (f, &l, g);
 868   return l;
 869 }
 870
 871 /* Round real toward zero to unsigned HOST_WIDE_INT
 872    implements  REAL_VALUE_UNSIGNED_FIX (x).
 873    Negative input returns zero.  */
 874
 875 unsigned HOST_WIDE_INT
 876 efixui (x)
 877      REAL_VALUE_TYPE x;
 878 {
 879   UEMUSHORT f[NE], g[NE];
 880   unsigned HOST_WIDE_INT l;
 881
 882   GET_REAL (&x, f);
 883 #ifdef NANS
 884   if (eisnan (f))
 885     {
 886       warning ("conversion from NaN to unsigned int");
 887       return (-1);
 888     }
 889 #endif
 890   euifrac (f, &l, g);
 891   return l;
 892 }
 893
 894
 895 /* REAL_VALUE_FROM_INT macro.  */
 896
 897 void
 898 ereal_from_int (d, i, j, mode)
 899      REAL_VALUE_TYPE *d;
 900      HOST_WIDE_INT i, j;
 901      enum machine_mode mode;
 902 {
 903   UEMUSHORT df[NE], dg[NE];
 904   HOST_WIDE_INT low, high;
 905   int sign;
 906
 907   if (GET_MODE_CLASS (mode) != MODE_FLOAT)
 908     abort ();
 909   sign = 0;
 910   low = i;
 911   if ((high = j) < 0)
 912     {
 913       sign = 1;
 914       /* complement and add 1 */
 915       high = ~high;
 916       if (low)
 917         low = -low;
 918       else
 919         high += 1;
 920     }
 921   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
 922   ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
 923   emul (dg, df, dg);
 924   ultoe ((unsigned HOST_WIDE_INT *) &low, df);
 925   eadd (df, dg, dg);
 926   if (sign)
 927     eneg (dg);
 928
 929   /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
 930      Avoid double-rounding errors later by rounding off now from the
 931      extra-wide internal format to the requested precision.  */
 932   switch (GET_MODE_BITSIZE (mode))
 933     {
 934     case 32:
 935       etoe24 (dg, df);
 936       e24toe (df, dg);
 937       break;
 938
 939     case 64:
 940       etoe53 (dg, df);
 941       e53toe (df, dg);
 942       break;
 943
 944     case 96:
 945       etoe64 (dg, df);
 946       e64toe (df, dg);
 947       break;
 948
 949     case 128:
 950 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
 951       etoe113 (dg, df);
 952       e113toe (df, dg);
 953 #else
 954       etoe64 (dg, df);
 955       e64toe (df, dg);
 956 #endif
 957       break;
 958
 959     default:
 960       abort ();
 961   }
 962
 963   PUT_REAL (dg, d);
 964 }
 965
 966
 967 /* REAL_VALUE_FROM_UNSIGNED_INT macro.  */
 968
 969 void
 970 ereal_from_uint (d, i, j, mode)
 971      REAL_VALUE_TYPE *d;
 972      unsigned HOST_WIDE_INT i, j;
 973      enum machine_mode mode;
 974 {
 975   UEMUSHORT df[NE], dg[NE];
 976   unsigned HOST_WIDE_INT low, high;
 977
 978   if (GET_MODE_CLASS (mode) != MODE_FLOAT)
 979     abort ();
 980   low = i;
 981   high = j;
 982   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
 983   ultoe (&high, dg);
 984   emul (dg, df, dg);
 985   ultoe (&low, df);
 986   eadd (df, dg, dg);
 987
 988   /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
 989      Avoid double-rounding errors later by rounding off now from the
 990      extra-wide internal format to the requested precision.  */
 991   switch (GET_MODE_BITSIZE (mode))
 992     {
 993     case 32:
 994       etoe24 (dg, df);
 995       e24toe (df, dg);
 996       break;
 997
 998     case 64:
 999       etoe53 (dg, df);
1000       e53toe (df, dg);
1001       break;
1002
1003     case 96:
1004       etoe64 (dg, df);
1005       e64toe (df, dg);
1006       break;
1007
1008     case 128:
1009 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1010       etoe113 (dg, df);
1011       e113toe (df, dg);
1012 #else
1013       etoe64 (dg, df);
1014       e64toe (df, dg);
1015 #endif
1016       break;
1017
1018     default:
1019       abort ();
1020   }
1021
1022   PUT_REAL (dg, d);
1023 }
1024
1025
1026 /* REAL_VALUE_TO_INT macro.  */
1027
1028 void
1029 ereal_to_int (low, high, rr)
1030      HOST_WIDE_INT *low, *high;
1031      REAL_VALUE_TYPE rr;
1032 {
1033   UEMUSHORT d[NE], df[NE], dg[NE], dh[NE];
1034   int s;
1035
1036   GET_REAL (&rr, d);
1037 #ifdef NANS
1038   if (eisnan (d))
1039     {
1040       warning ("conversion from NaN to int");
1041       *low = -1;
1042       *high = -1;
1043       return;
1044     }
1045 #endif
1046   /* convert positive value */
1047   s = 0;
1048   if (eisneg (d))
1049     {
1050       eneg (d);
1051       s = 1;
1052     }
1053   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
1054   ediv (df, d, dg);             /* dg = d / 2^32 is the high word */
1055   euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
1056   emul (df, dh, dg);            /* fractional part is the low word */
1057   euifrac (dg, (unsigned HOST_WIDE_INT *) low, dh);
1058   if (s)
1059     {
1060       /* complement and add 1 */
1061       *high = ~(*high);
1062       if (*low)
1063         *low = -(*low);
1064       else
1065         *high += 1;
1066     }
1067 }
1068
1069
1070 /* REAL_VALUE_LDEXP macro.  */
1071
1072 REAL_VALUE_TYPE
1073 ereal_ldexp (x, n)
1074      REAL_VALUE_TYPE x;
1075      int n;
1076 {
1077   UEMUSHORT e[NE], y[NE];
1078   REAL_VALUE_TYPE r;
1079
1080   GET_REAL (&x, e);
1081 #ifdef NANS
1082   if (eisnan (e))
1083     return (x);
1084 #endif
1085   eldexp (e, n, y);
1086   PUT_REAL (y, &r);
1087   return (r);
1088 }
1089
1090 /* Check for infinity in a REAL_VALUE_TYPE.  */
1091
1092 int
1093 target_isinf (x)
1094      REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
1095 {
1096 #ifdef INFINITY
1097   UEMUSHORT e[NE];
1098
1099   GET_REAL (&x, e);
1100   return (eisinf (e));
1101 #else
1102   return 0;
1103 #endif
1104 }
1105
1106 /* Check whether a REAL_VALUE_TYPE item is a NaN.  */
1107
1108 int
1109 target_isnan (x)
1110      REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
1111 {
1112 #ifdef NANS
1113   UEMUSHORT e[NE];
1114
1115   GET_REAL (&x, e);
1116   return (eisnan (e));
1117 #else
1118   return (0);
1119 #endif
1120 }
1121
1122
1123 /* Check for a negative REAL_VALUE_TYPE number.
1124    This just checks the sign bit, so that -0 counts as negative.  */
1125
1126 int
1127 target_negative (x)
1128      REAL_VALUE_TYPE x;
1129 {
1130   return ereal_isneg (x);
1131 }
1132
1133 /* Expansion of REAL_VALUE_TRUNCATE.
1134    The result is in floating point, rounded to nearest or even.  */
1135
1136 REAL_VALUE_TYPE
1137 real_value_truncate (mode, arg)
1138      enum machine_mode mode;
1139      REAL_VALUE_TYPE arg;
1140 {
1141   UEMUSHORT e[NE], t[NE];
1142   REAL_VALUE_TYPE r;
1143
1144   GET_REAL (&arg, e);
1145 #ifdef NANS
1146   if (eisnan (e))
1147     return (arg);
1148 #endif
1149   eclear (t);
1150   switch (mode)
1151     {
1152     case TFmode:
1153 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1154       etoe113 (e, t);
1155       e113toe (t, t);
1156       break;
1157 #endif
1158       /* FALLTHRU */
1159
1160     case XFmode:
1161       etoe64 (e, t);
1162       e64toe (t, t);
1163       break;
1164
1165     case DFmode:
1166       etoe53 (e, t);
1167       e53toe (t, t);
1168       break;
1169
1170     case SFmode:
1171 #ifndef C4X
1172     case HFmode:
1173 #endif
1174       etoe24 (e, t);
1175       e24toe (t, t);
1176       break;
1177
1178 #ifdef C4X
1179     case HFmode:
1180     case QFmode:
1181       etoe53 (e, t);
1182       e53toe (t, t);
1183       break;
1184 #endif
1185
1186     case SImode:
1187       r = etrunci (arg);
1188       return (r);
1189
1190     /* If an unsupported type was requested, presume that
1191        the machine files know something useful to do with
1192        the unmodified value.  */
1193
1194     default:
1195       return (arg);
1196     }
1197   PUT_REAL (t, &r);
1198   return (r);
1199 }
1200
1201 /* Return true if ARG can be represented exactly in MODE.  */
1202
1203 bool
1204 exact_real_truncate (mode, arg)
1205      enum machine_mode mode;
1206      REAL_VALUE_TYPE *arg;
1207 {
1208   REAL_VALUE_TYPE trunc;
1209
1210   if (target_isnan (*arg))
1211     return false;
1212
1213   trunc = real_value_truncate (mode, *arg);
1214   return ereal_cmp (*arg, trunc) == 0;
1215 }
1216
1217 /* Try to change R into its exact multiplicative inverse in machine mode
1218    MODE.  Return nonzero function value if successful.  */
1219
1220 int
1221 exact_real_inverse (mode, r)
1222      enum machine_mode mode;
1223      REAL_VALUE_TYPE *r;
1224 {
1225   UEMUSHORT e[NE], einv[NE];
1226   REAL_VALUE_TYPE rinv;
1227   int i;
1228
1229   GET_REAL (r, e);
1230
1231   /* Test for input in range.  Don't transform IEEE special values.  */
1232   if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
1233     return 0;
1234
1235   /* Test for a power of 2: all significand bits zero except the MSB.
1236      We are assuming the target has binary (or hex) arithmetic.  */
1237   if (e[NE - 2] != 0x8000)
1238     return 0;
1239
1240   for (i = 0; i < NE - 2; i++)
1241     {
1242       if (e[i] != 0)
1243         return 0;
1244     }
1245
1246   /* Compute the inverse and truncate it to the required mode.  */
1247   ediv (e, eone, einv);
1248   PUT_REAL (einv, &rinv);
1249   rinv = real_value_truncate (mode, rinv);
1250
1251 #ifdef CHECK_FLOAT_VALUE
1252   /* This check is not redundant.  It may, for example, flush
1253      a supposedly IEEE denormal value to zero.  */
1254   i = 0;
1255   if (CHECK_FLOAT_VALUE (mode, rinv, i))
1256     return 0;
1257 #endif
1258   GET_REAL (&rinv, einv);
1259
1260   /* Check the bits again, because the truncation might have
1261      generated an arbitrary saturation value on overflow.  */
1262   if (einv[NE - 2] != 0x8000)
1263     return 0;
1264
1265   for (i = 0; i < NE - 2; i++)
1266     {
1267       if (einv[i] != 0)
1268         return 0;
1269     }
1270
1271   /* Fail if the computed inverse is out of range.  */
1272   if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
1273     return 0;
1274
1275   /* Output the reciprocal and return success flag.  */
1276   PUT_REAL (einv, r);
1277   return 1;
1278 }
1279
1280 /* Used for debugging--print the value of R in human-readable format
1281    on stderr.  */
1282
1283 void
1284 debug_real (r)
1285      REAL_VALUE_TYPE r;
1286 {
1287   char dstr[30];
1288
1289   REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1290   fprintf (stderr, "%s", dstr);
1291 }
1292
1293 \f
1294 /* The following routines convert REAL_VALUE_TYPE to the various floating
1295    point formats that are meaningful to supported computers.
1296
1297    The results are returned in 32-bit pieces, each piece stored in a `long'.
1298    This is so they can be printed by statements like
1299
1300       fprintf (file, "%lx, %lx", L[0],  L[1]);
1301
1302    that will work on both narrow- and wide-word host computers.  */
1303
1304 /* Convert R to a 128-bit long double precision value.  The output array L
1305    contains four 32-bit pieces of the result, in the order they would appear
1306    in memory.  */
1307
1308 void
1309 etartdouble (r, l)
1310      REAL_VALUE_TYPE r;
1311      long l[];
1312 {
1313   UEMUSHORT e[NE];
1314
1315   GET_REAL (&r, e);
1316 #if INTEL_EXTENDED_IEEE_FORMAT == 0
1317   etoe113 (e, e);
1318 #else
1319   etoe64 (e, e);
1320 #endif
1321   endian (e, l, TFmode);
1322 }
1323
1324 /* Convert R to a double extended precision value.  The output array L
1325    contains three 32-bit pieces of the result, in the order they would
1326    appear in memory.  */
1327
1328 void
1329 etarldouble (r, l)
1330      REAL_VALUE_TYPE r;
1331      long l[];
1332 {
1333   UEMUSHORT e[NE];
1334
1335   GET_REAL (&r, e);
1336   etoe64 (e, e);
1337   endian (e, l, XFmode);
1338 }
1339
1340 /* Convert R to a double precision value.  The output array L contains two
1341    32-bit pieces of the result, in the order they would appear in memory.  */
1342
1343 void
1344 etardouble (r, l)
1345      REAL_VALUE_TYPE r;
1346      long l[];
1347 {
1348   UEMUSHORT e[NE];
1349
1350   GET_REAL (&r, e);
1351   etoe53 (e, e);
1352   endian (e, l, DFmode);
1353 }
1354
1355 /* Convert R to a single precision float value stored in the least-significant
1356    bits of a `long'.  */
1357
1358 long
1359 etarsingle (r)
1360      REAL_VALUE_TYPE r;
1361 {
1362   UEMUSHORT e[NE];
1363   long l;
1364
1365   GET_REAL (&r, e);
1366   etoe24 (e, e);
1367   endian (e, &l, SFmode);
1368   return ((long) l);
1369 }
1370
1371 /* Convert X to a decimal ASCII string S for output to an assembly
1372    language file.  Note, there is no standard way to spell infinity or
1373    a NaN, so these values may require special treatment in the tm.h
1374    macros.  */
1375
1376 void
1377 ereal_to_decimal (x, s)
1378      REAL_VALUE_TYPE x;
1379      char *s;
1380 {
1381   UEMUSHORT e[NE];
1382
1383   GET_REAL (&x, e);
1384   etoasc (e, s, 20);
1385 }
1386
1387 /* Compare X and Y.  Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1388    or -2 if either is a NaN.  */
1389
1390 int
1391 ereal_cmp (x, y)
1392      REAL_VALUE_TYPE x, y;
1393 {
1394   UEMUSHORT ex[NE], ey[NE];
1395
1396   GET_REAL (&x, ex);
1397   GET_REAL (&y, ey);
1398   return (ecmp (ex, ey));
1399 }
1400
1401 /*  Return 1 if the sign bit of X is set, else return 0.  */
1402
1403 int
1404 ereal_isneg (x)
1405      REAL_VALUE_TYPE x;
1406 {
1407   UEMUSHORT ex[NE];
1408
1409   GET_REAL (&x, ex);
1410   return (eisneg (ex));
1411 }
1412
1413 \f
1414 /*
1415   Extended precision IEEE binary floating point arithmetic routines
1416
1417   Numbers are stored in C language as arrays of 16-bit unsigned
1418   short integers.  The arguments of the routines are pointers to
1419   the arrays.
1420
1421   External e type data structure, similar to Intel 8087 chip
1422   temporary real format but possibly with a larger significand:
1423
1424         NE-1 significand words  (least significant word first,
1425                                  most significant bit is normally set)
1426         exponent                (value = EXONE for 1.0,
1427                                 top bit is the sign)
1428
1429
1430   Internal exploded e-type data structure of a number (a "word" is 16 bits):
1431
1432   ei[0] sign word       (0 for positive, 0xffff for negative)
1433   ei[1] biased exponent (value = EXONE for the number 1.0)
1434   ei[2] high guard word (always zero after normalization)
1435   ei[3]
1436   to ei[NI-2]   significand     (NI-4 significand words,
1437                                  most significant word first,
1438                                  most significant bit is set)
1439   ei[NI-1]      low guard word  (0x8000 bit is rounding place)
1440
1441
1442
1443                 Routines for external format e-type numbers
1444
1445         asctoe (string, e)      ASCII string to extended double e type
1446         asctoe64 (string, &d)   ASCII string to long double
1447         asctoe53 (string, &d)   ASCII string to double
1448         asctoe24 (string, &f)   ASCII string to single
1449         asctoeg (string, e, prec) ASCII string to specified precision
1450         e24toe (&f, e)          IEEE single precision to e type
1451         e53toe (&d, e)          IEEE double precision to e type
1452         e64toe (&d, e)          IEEE long double precision to e type
1453         e113toe (&d, e)         128-bit long double precision to e type
1454 #if 0
1455         eabs (e)                        absolute value
1456 #endif
1457         eadd (a, b, c)          c = b + a
1458         eclear (e)              e = 0
1459         ecmp (a, b)             Returns 1 if a > b, 0 if a == b,
1460                                 -1 if a < b, -2 if either a or b is a NaN.
1461         ediv (a, b, c)          c = b / a
1462         efloor (a, b)           truncate to integer, toward -infinity
1463         efrexp (a, exp, s)      extract exponent and significand
1464         eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
1465         euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
1466         einfin (e)              set e to infinity, leaving its sign alone
1467         eldexp (a, n, b)        multiply by 2**n
1468         emov (a, b)             b = a
1469         emul (a, b, c)          c = b * a
1470         eneg (e)                        e = -e
1471 #if 0
1472         eround (a, b)           b = nearest integer value to a
1473 #endif
1474         esub (a, b, c)          c = b - a
1475 #if 0
1476         e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
1477         e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
1478         e64toasc (&d, str, n)   80-bit long double to ASCII string
1479         e113toasc (&d, str, n)  128-bit long double to ASCII string
1480 #endif
1481         etoasc (e, str, n)      e to ASCII string, n digits after decimal
1482         etoe24 (e, &f)          convert e type to IEEE single precision
1483         etoe53 (e, &d)          convert e type to IEEE double precision
1484         etoe64 (e, &d)          convert e type to IEEE long double precision
1485         ltoe (&l, e)            HOST_WIDE_INT to e type
1486         ultoe (&l, e)           unsigned HOST_WIDE_INT to e type
1487         eisneg (e)              1 if sign bit of e != 0, else 0
1488         eisinf (e)              1 if e has maximum exponent (non-IEEE)
1489                                 or is infinite (IEEE)
1490         eisnan (e)              1 if e is a NaN
1491
1492
1493                 Routines for internal format exploded e-type numbers
1494
1495         eaddm (ai, bi)          add significands, bi = bi + ai
1496         ecleaz (ei)             ei = 0
1497         ecleazs (ei)            set ei = 0 but leave its sign alone
1498         ecmpm (ai, bi)          compare significands, return 1, 0, or -1
1499         edivm (ai, bi)          divide  significands, bi = bi / ai
1500         emdnorm (ai,l,s,exp)    normalize and round off
1501         emovi (a, ai)           convert external a to internal ai
1502         emovo (ai, a)           convert internal ai to external a
1503         emovz (ai, bi)          bi = ai, low guard word of bi = 0
1504         emulm (ai, bi)          multiply significands, bi = bi * ai
1505         enormlz (ei)            left-justify the significand
1506         eshdn1 (ai)             shift significand and guards down 1 bit
1507         eshdn8 (ai)             shift down 8 bits
1508         eshdn6 (ai)             shift down 16 bits
1509         eshift (ai, n)          shift ai n bits up (or down if n < 0)
1510         eshup1 (ai)             shift significand and guards up 1 bit
1511         eshup8 (ai)             shift up 8 bits
1512         eshup6 (ai)             shift up 16 bits
1513         esubm (ai, bi)          subtract significands, bi = bi - ai
1514         eiisinf (ai)            1 if infinite
1515         eiisnan (ai)            1 if a NaN
1516         eiisneg (ai)            1 if sign bit of ai != 0, else 0
1517         einan (ai)              set ai = NaN
1518 #if 0
1519         eiinfin (ai)            set ai = infinity
1520 #endif
1521
1522   The result is always normalized and rounded to NI-4 word precision
1523   after each arithmetic operation.
1524
1525   Exception flags are NOT fully supported.
1526
1527   Signaling NaN's are NOT supported; they are treated the same
1528   as quiet NaN's.
1529
1530   Define INFINITY for support of infinity; otherwise a
1531   saturation arithmetic is implemented.
1532
1533   Define NANS for support of Not-a-Number items; otherwise the
1534   arithmetic will never produce a NaN output, and might be confused
1535   by a NaN input.
1536   If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1537   either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1538   may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1539   if in doubt.
1540
1541   Denormals are always supported here where appropriate (e.g., not
1542   for conversion to DEC numbers).  */
1543
1544 /* Definitions for error codes that are passed to the common error handling
1545    routine mtherr.
1546
1547    For Digital Equipment PDP-11 and VAX computers, certain
1548   IBM systems, and others that use numbers with a 56-bit
1549   significand, the symbol DEC should be defined.  In this
1550   mode, most floating point constants are given as arrays
1551   of octal integers to eliminate decimal to binary conversion
1552   errors that might be introduced by the compiler.
1553
1554   For computers, such as IBM PC, that follow the IEEE
1555   Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1556   Std 754-1985), the symbol IEEE should be defined.
1557   These numbers have 53-bit significands.  In this mode, constants
1558   are provided as arrays of hexadecimal 16 bit integers.
1559   The endian-ness of generated values is controlled by
1560   REAL_WORDS_BIG_ENDIAN.
1561
1562   To accommodate other types of computer arithmetic, all
1563   constants are also provided in a normal decimal radix
1564   which one can hope are correctly converted to a suitable
1565   format by the available C language compiler.  To invoke
1566   this mode, the symbol UNK is defined.
1567
1568   An important difference among these modes is a predefined
1569   set of machine arithmetic constants for each.  The numbers
1570   MACHEP (the machine roundoff error), MAXNUM (largest number
1571   represented), and several other parameters are preset by
1572   the configuration symbol.  Check the file const.c to
1573   ensure that these values are correct for your computer.
1574
1575   For ANSI C compatibility, define ANSIC equal to 1.  Currently
1576   this affects only the atan2 function and others that use it.  */
1577
1578 /* Constant definitions for math error conditions.  */
1579
1580 #define DOMAIN          1       /* argument domain error */
1581 #define SING            2       /* argument singularity */
1582 #define OVERFLOW        3       /* overflow range error */
1583 #define UNDERFLOW       4       /* underflow range error */
1584 #define TLOSS           5       /* total loss of precision */
1585 #define PLOSS           6       /* partial loss of precision */
1586 #define INVALID         7       /* NaN-producing operation */
1587
1588 /*  e type constants used by high precision check routines */
1589
1590 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1591 /* 0.0 */
1592 const UEMUSHORT ezero[NE] =
1593  {0x0000, 0x0000, 0x0000, 0x0000,
1594   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1595
1596 /* 5.0E-1 */
1597 const UEMUSHORT ehalf[NE] =
1598  {0x0000, 0x0000, 0x0000, 0x0000,
1599   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1600
1601 /* 1.0E0 */
1602 const UEMUSHORT eone[NE] =
1603  {0x0000, 0x0000, 0x0000, 0x0000,
1604   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1605
1606 /* 2.0E0 */
1607 const UEMUSHORT etwo[NE] =
1608  {0x0000, 0x0000, 0x0000, 0x0000,
1609   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1610
1611 /* 3.2E1 */
1612 const UEMUSHORT e32[NE] =
1613  {0x0000, 0x0000, 0x0000, 0x0000,
1614   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1615
1616 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1617 const UEMUSHORT elog2[NE] =
1618  {0x40f3, 0xf6af, 0x03f2, 0xb398,
1619   0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1620
1621 /* 1.41421356237309504880168872420969807856967187537695E0 */
1622 const UEMUSHORT esqrt2[NE] =
1623  {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1624   0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1625
1626 /* 3.14159265358979323846264338327950288419716939937511E0 */
1627 const UEMUSHORT epi[NE] =
1628  {0x2902, 0x1cd1, 0x80dc, 0x628b,
1629   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1630
1631 #else
1632 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1633 const UEMUSHORT ezero[NE] =
1634  {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1635 const UEMUSHORT ehalf[NE] =
1636  {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1637 const UEMUSHORT eone[NE] =
1638  {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1639 const UEMUSHORT etwo[NE] =
1640  {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1641 const UEMUSHORT e32[NE] =
1642  {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1643 const UEMUSHORT elog2[NE] =
1644  {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1645 const UEMUSHORT esqrt2[NE] =
1646  {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1647 const UEMUSHORT epi[NE] =
1648  {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1649 #endif
1650
1651 /* Control register for rounding precision.
1652    This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */
1653
1654 int rndprc = NBITS;
1655 extern int rndprc;
1656
1657 /*  Clear out entire e-type number X.  */
1658
1659 static void
1660 eclear (x)
1661      UEMUSHORT *x;
1662 {
1663   int i;
1664
1665   for (i = 0; i < NE; i++)
1666     *x++ = 0;
1667 }
1668
1669 /* Move e-type number from A to B.  */
1670
1671 static void
1672 emov (a, b)
1673      const UEMUSHORT *a;
1674      UEMUSHORT *b;
1675 {
1676   int i;
1677
1678   for (i = 0; i < NE; i++)
1679     *b++ = *a++;
1680 }
1681
1682
1683 #if 0
1684 /* Absolute value of e-type X.  */
1685
1686 static void
1687 eabs (x)
1688      UEMUSHORT x[];
1689 {
1690   /* sign is top bit of last word of external format */
1691   x[NE - 1] &= 0x7fff;
1692 }
1693 #endif /* 0 */
1694
1695 /* Negate the e-type number X.  */
1696
1697 static void
1698 eneg (x)
1699      UEMUSHORT x[];
1700 {
1701
1702   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1703 }
1704
1705 /* Return 1 if sign bit of e-type number X is nonzero, else zero.  */
1706
1707 static int
1708 eisneg (x)
1709      const UEMUSHORT x[];
1710 {
1711
1712   if (x[NE - 1] & 0x8000)
1713     return (1);
1714   else
1715     return (0);
1716 }
1717
1718 /* Return 1 if e-type number X is infinity, else return zero.  */
1719
1720 static int
1721 eisinf (x)
1722      const UEMUSHORT x[];
1723 {
1724
1725 #ifdef NANS
1726   if (eisnan (x))
1727     return (0);
1728 #endif
1729   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1730     return (1);
1731   else
1732     return (0);
1733 }
1734
1735 /* Check if e-type number is not a number.  The bit pattern is one that we
1736    defined, so we know for sure how to detect it.  */
1737
1738 static int
1739 eisnan (x)
1740      const UEMUSHORT x[] ATTRIBUTE_UNUSED;
1741 {
1742 #ifdef NANS
1743   int i;
1744
1745   /* NaN has maximum exponent */
1746   if ((x[NE - 1] & 0x7fff) != 0x7fff)
1747     return (0);
1748   /* ... and non-zero significand field.  */
1749   for (i = 0; i < NE - 1; i++)
1750     {
1751       if (*x++ != 0)
1752         return (1);
1753     }
1754 #endif
1755
1756   return (0);
1757 }
1758
1759 /*  Fill e-type number X with infinity pattern (IEEE)
1760     or largest possible number (non-IEEE).  */
1761
1762 static void
1763 einfin (x)
1764      UEMUSHORT *x;
1765 {
1766   int i;
1767
1768 #ifdef INFINITY
1769   for (i = 0; i < NE - 1; i++)
1770     *x++ = 0;
1771   *x |= 32767;
1772 #else
1773   for (i = 0; i < NE - 1; i++)
1774     *x++ = 0xffff;
1775   *x |= 32766;
1776   if (rndprc < NBITS)
1777     {
1778       if (rndprc == 113)
1779         {
1780           *(x - 9) = 0;
1781           *(x - 8) = 0;
1782         }
1783       if (rndprc == 64)
1784         {
1785           *(x - 5) = 0;
1786         }
1787       if (rndprc == 53)
1788         {
1789           *(x - 4) = 0xf800;
1790         }
1791       else
1792         {
1793           *(x - 4) = 0;
1794           *(x - 3) = 0;
1795           *(x - 2) = 0xff00;
1796         }
1797     }
1798 #endif
1799 }
1800
1801 /* Output an e-type NaN.
1802    This generates Intel's quiet NaN pattern for extended real.
1803    The exponent is 7fff, the leading mantissa word is c000.  */
1804
1805 #ifdef NANS
1806 static void
1807 enan (x, sign)
1808      UEMUSHORT *x;
1809      int sign;
1810 {
1811   int i;
1812
1813   for (i = 0; i < NE - 2; i++)
1814     *x++ = 0;
1815   *x++ = 0xc000;
1816   *x = (sign << 15) | 0x7fff;
1817 }
1818 #endif /* NANS */
1819
1820 /* Move in an e-type number A, converting it to exploded e-type B.  */
1821
1822 static void
1823 emovi (a, b)
1824      const UEMUSHORT *a;
1825      UEMUSHORT *b;
1826 {
1827   const UEMUSHORT *p;
1828   UEMUSHORT *q;
1829   int i;
1830
1831   q = b;
1832   p = a + (NE - 1);             /* point to last word of external number */
1833   /* get the sign bit */
1834   if (*p & 0x8000)
1835     *q++ = 0xffff;
1836   else
1837     *q++ = 0;
1838   /* get the exponent */
1839   *q = *p--;
1840   *q++ &= 0x7fff;               /* delete the sign bit */
1841 #ifdef INFINITY
1842   if ((*(q - 1) & 0x7fff) == 0x7fff)
1843     {
1844 #ifdef NANS
1845       if (eisnan (a))
1846         {
1847           *q++ = 0;
1848           for (i = 3; i < NI; i++)
1849             *q++ = *p--;
1850           return;
1851         }
1852 #endif
1853
1854       for (i = 2; i < NI; i++)
1855         *q++ = 0;
1856       return;
1857     }
1858 #endif
1859
1860   /* clear high guard word */
1861   *q++ = 0;
1862   /* move in the significand */
1863   for (i = 0; i < NE - 1; i++)
1864     *q++ = *p--;
1865   /* clear low guard word */
1866   *q = 0;
1867 }
1868
1869 /* Move out exploded e-type number A, converting it to e type B.  */
1870
1871 static void
1872 emovo (a, b)
1873      const UEMUSHORT *a;
1874      UEMUSHORT *b;
1875 {
1876   const UEMUSHORT *p;
1877   UEMUSHORT *q;
1878   UEMUSHORT i;
1879   int j;
1880
1881   p = a;
1882   q = b + (NE - 1);             /* point to output exponent */
1883   /* combine sign and exponent */
1884   i = *p++;
1885   if (i)
1886     *q-- = *p++ | 0x8000;
1887   else
1888     *q-- = *p++;
1889 #ifdef INFINITY
1890   if (*(p - 1) == 0x7fff)
1891     {
1892 #ifdef NANS
1893       if (eiisnan (a))
1894         {
1895           enan (b, eiisneg (a));
1896           return;
1897         }
1898 #endif
1899       einfin (b);
1900         return;
1901     }
1902 #endif
1903   /* skip over guard word */
1904   ++p;
1905   /* move the significand */
1906   for (j = 0; j < NE - 1; j++)
1907     *q-- = *p++;
1908 }
1909
1910 /* Clear out exploded e-type number XI.  */
1911
1912 static void
1913 ecleaz (xi)
1914      UEMUSHORT *xi;
1915 {
1916   int i;
1917
1918   for (i = 0; i < NI; i++)
1919     *xi++ = 0;
1920 }
1921
1922 /* Clear out exploded e-type XI, but don't touch the sign.  */
1923
1924 static void
1925 ecleazs (xi)
1926      UEMUSHORT *xi;
1927 {
1928   int i;
1929
1930   ++xi;
1931   for (i = 0; i < NI - 1; i++)
1932     *xi++ = 0;
1933 }
1934
1935 /* Move exploded e-type number from A to B.  */
1936
1937 static void
1938 emovz (a, b)
1939      const UEMUSHORT *a;
1940      UEMUSHORT *b;
1941 {
1942   int i;
1943
1944   for (i = 0; i < NI - 1; i++)
1945     *b++ = *a++;
1946   /* clear low guard word */
1947   *b = 0;
1948 }
1949
1950 /* Generate exploded e-type NaN.
1951    The explicit pattern for this is maximum exponent and
1952    top two significant bits set.  */
1953
1954 #ifdef NANS
1955 static void
1956 einan (x)
1957      UEMUSHORT x[];
1958 {
1959
1960   ecleaz (x);
1961   x[E] = 0x7fff;
1962   x[M + 1] = 0xc000;
1963 }
1964 #endif /* NANS */
1965
1966 /* Return nonzero if exploded e-type X is a NaN.  */
1967
1968 #ifdef NANS
1969 static int
1970 eiisnan (x)
1971      const UEMUSHORT x[];
1972 {
1973   int i;
1974
1975   if ((x[E] & 0x7fff) == 0x7fff)
1976     {
1977       for (i = M + 1; i < NI; i++)
1978         {
1979           if (x[i] != 0)
1980             return (1);
1981         }
1982     }
1983   return (0);
1984 }
1985 #endif /* NANS */
1986
1987 /* Return nonzero if sign of exploded e-type X is nonzero.  */
1988
1989 static int
1990 eiisneg (x)
1991      const UEMUSHORT x[];
1992 {
1993
1994   return x[0] != 0;
1995 }
1996
1997 #if 0
1998 /* Fill exploded e-type X with infinity pattern.
1999    This has maximum exponent and significand all zeros.  */
2000
2001 static void
2002 eiinfin (x)
2003      UEMUSHORT x[];
2004 {
2005
2006   ecleaz (x);
2007   x[E] = 0x7fff;
2008 }
2009 #endif /* 0 */
2010
2011 /* Return nonzero if exploded e-type X is infinite.  */
2012
2013 #ifdef INFINITY
2014 static int
2015 eiisinf (x)
2016      const UEMUSHORT x[];
2017 {
2018
2019 #ifdef NANS
2020   if (eiisnan (x))
2021     return (0);
2022 #endif
2023   if ((x[E] & 0x7fff) == 0x7fff)
2024     return (1);
2025   return (0);
2026 }
2027 #endif /* INFINITY */
2028
2029 /* Compare significands of numbers in internal exploded e-type format.
2030    Guard words are included in the comparison.
2031
2032    Returns      +1 if a > b
2033                  0 if a == b
2034                 -1 if a < b   */
2035
2036 static int
2037 ecmpm (a, b)
2038      const UEMUSHORT *a, *b;
2039 {
2040   int i;
2041
2042   a += M;                       /* skip up to significand area */
2043   b += M;
2044   for (i = M; i < NI; i++)
2045     {
2046       if (*a++ != *b++)
2047         goto difrnt;
2048     }
2049   return (0);
2050
2051  difrnt:
2052   if (*(--a) > *(--b))
2053     return (1);
2054   else
2055     return (-1);
2056 }
2057
2058 /* Shift significand of exploded e-type X down by 1 bit.  */
2059
2060 static void
2061 eshdn1 (x)
2062      UEMUSHORT *x;
2063 {
2064   UEMUSHORT bits;
2065   int i;
2066
2067   x += M;                       /* point to significand area */
2068
2069   bits = 0;
2070   for (i = M; i < NI; i++)
2071     {
2072       if (*x & 1)
2073         bits |= 1;
2074       *x >>= 1;
2075       if (bits & 2)
2076         *x |= 0x8000;
2077       bits <<= 1;
2078       ++x;
2079     }
2080 }
2081
2082 /* Shift significand of exploded e-type X up by 1 bit.  */
2083
2084 static void
2085 eshup1 (x)
2086      UEMUSHORT *x;
2087 {
2088   UEMUSHORT bits;
2089   int i;
2090
2091   x += NI - 1;
2092   bits = 0;
2093
2094   for (i = M; i < NI; i++)
2095     {
2096       if (*x & 0x8000)
2097         bits |= 1;
2098       *x <<= 1;
2099       if (bits & 2)
2100         *x |= 1;
2101       bits <<= 1;
2102       --x;
2103     }
2104 }
2105
2106
2107 /* Shift significand of exploded e-type X down by 8 bits.  */
2108
2109 static void
2110 eshdn8 (x)
2111      UEMUSHORT *x;
2112 {
2113   UEMUSHORT newbyt, oldbyt;
2114   int i;
2115
2116   x += M;
2117   oldbyt = 0;
2118   for (i = M; i < NI; i++)
2119     {
2120       newbyt = *x << 8;
2121       *x >>= 8;
2122       *x |= oldbyt;
2123       oldbyt = newbyt;
2124       ++x;
2125     }
2126 }
2127
2128 /* Shift significand of exploded e-type X up by 8 bits.  */
2129
2130 static void
2131 eshup8 (x)
2132      UEMUSHORT *x;
2133 {
2134   int i;
2135   UEMUSHORT newbyt, oldbyt;
2136
2137   x += NI - 1;
2138   oldbyt = 0;
2139
2140   for (i = M; i < NI; i++)
2141     {
2142       newbyt = *x >> 8;
2143       *x <<= 8;
2144       *x |= oldbyt;
2145       oldbyt = newbyt;
2146       --x;
2147     }
2148 }
2149
2150 /* Shift significand of exploded e-type X up by 16 bits.  */
2151
2152 static void
2153 eshup6 (x)
2154      UEMUSHORT *x;
2155 {
2156   int i;
2157   UEMUSHORT *p;
2158
2159   p = x + M;
2160   x += M + 1;
2161
2162   for (i = M; i < NI - 1; i++)
2163     *p++ = *x++;
2164
2165   *p = 0;
2166 }
2167
2168 /* Shift significand of exploded e-type X down by 16 bits.  */
2169
2170 static void
2171 eshdn6 (x)
2172      UEMUSHORT *x;
2173 {
2174   int i;
2175   UEMUSHORT *p;
2176
2177   x += NI - 1;
2178   p = x + 1;
2179
2180   for (i = M; i < NI - 1; i++)
2181     *(--p) = *(--x);
2182
2183   *(--p) = 0;
2184 }
2185
2186 /* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */
2187
2188 static void
2189 eaddm (x, y)
2190      const UEMUSHORT *x;
2191      UEMUSHORT *y;
2192 {
2193   unsigned EMULONG a;
2194   int i;
2195   unsigned int carry;
2196
2197   x += NI - 1;
2198   y += NI - 1;
2199   carry = 0;
2200   for (i = M; i < NI; i++)
2201     {
2202       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
2203       if (a & 0x10000)
2204         carry = 1;
2205       else
2206         carry = 0;
2207       *y = (UEMUSHORT) a;
2208       --x;
2209       --y;
2210     }
2211 }
2212
2213 /* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */
2214
2215 static void
2216 esubm (x, y)
2217      const UEMUSHORT *x;
2218      UEMUSHORT *y;
2219 {
2220   unsigned EMULONG a;
2221   int i;
2222   unsigned int carry;
2223
2224   x += NI - 1;
2225   y += NI - 1;
2226   carry = 0;
2227   for (i = M; i < NI; i++)
2228     {
2229       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
2230       if (a & 0x10000)
2231         carry = 1;
2232       else
2233         carry = 0;
2234       *y = (UEMUSHORT) a;
2235       --x;
2236       --y;
2237     }
2238 }
2239
2240
2241 static UEMUSHORT equot[NI];
2242
2243
2244 #if 0
2245 /* Radix 2 shift-and-add versions of multiply and divide  */
2246
2247
2248 /* Divide significands */
2249
2250 int
2251 edivm (den, num)
2252      UEMUSHORT den[], num[];
2253 {
2254   int i;
2255   UEMUSHORT *p, *q;
2256   UEMUSHORT j;
2257
2258   p = &equot[0];
2259   *p++ = num[0];
2260   *p++ = num[1];
2261
2262   for (i = M; i < NI; i++)
2263     {
2264       *p++ = 0;
2265     }
2266
2267   /* Use faster compare and subtraction if denominator has only 15 bits of
2268      significance.  */
2269
2270   p = &den[M + 2];
2271   if (*p++ == 0)
2272     {
2273       for (i = M + 3; i < NI; i++)
2274         {
2275           if (*p++ != 0)
2276             goto fulldiv;
2277         }
2278       if ((den[M + 1] & 1) != 0)
2279         goto fulldiv;
2280       eshdn1 (num);
2281       eshdn1 (den);
2282
2283       p = &den[M + 1];
2284       q = &num[M + 1];
2285
2286       for (i = 0; i < NBITS + 2; i++)
2287         {
2288           if (*p <= *q)
2289             {
2290               *q -= *p;
2291               j = 1;
2292             }
2293           else
2294             {
2295               j = 0;
2296             }
2297           eshup1 (equot);
2298           equot[NI - 2] |= j;
2299           eshup1 (num);
2300         }
2301       goto divdon;
2302     }
2303
2304   /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2305      bit + 1 roundoff bit.  */
2306
2307  fulldiv:
2308
2309   p = &equot[NI - 2];
2310   for (i = 0; i < NBITS + 2; i++)
2311     {
2312       if (ecmpm (den, num) <= 0)
2313         {
2314           esubm (den, num);
2315           j = 1;                /* quotient bit = 1 */
2316         }
2317       else
2318         j = 0;
2319       eshup1 (equot);
2320       *p |= j;
2321       eshup1 (num);
2322     }
2323
2324  divdon:
2325
2326   eshdn1 (equot);
2327   eshdn1 (equot);
2328
2329   /* test for nonzero remainder after roundoff bit */
2330   p = &num[M];
2331   j = 0;
2332   for (i = M; i < NI; i++)
2333     {
2334       j |= *p++;
2335     }
2336   if (j)
2337     j = 1;
2338
2339
2340   for (i = 0; i < NI; i++)
2341     num[i] = equot[i];
2342   return ((int) j);
2343 }
2344
2345
2346 /* Multiply significands */
2347
2348 int
2349 emulm (a, b)
2350      UEMUSHORT a[], b[];
2351 {
2352   UEMUSHORT *p, *q;
2353   int i, j, k;
2354
2355   equot[0] = b[0];
2356   equot[1] = b[1];
2357   for (i = M; i < NI; i++)
2358     equot[i] = 0;
2359
2360   p = &a[NI - 2];
2361   k = NBITS;
2362   while (*p == 0)               /* significand is not supposed to be zero */
2363     {
2364       eshdn6 (a);
2365       k -= 16;
2366     }
2367   if ((*p & 0xff) == 0)
2368     {
2369       eshdn8 (a);
2370       k -= 8;
2371     }
2372
2373   q = &equot[NI - 1];
2374   j = 0;
2375   for (i = 0; i < k; i++)
2376     {
2377       if (*p & 1)
2378         eaddm (b, equot);
2379       /* remember if there were any nonzero bits shifted out */
2380       if (*q & 1)
2381         j |= 1;
2382       eshdn1 (a);
2383       eshdn1 (equot);
2384     }
2385
2386   for (i = 0; i < NI; i++)
2387     b[i] = equot[i];
2388
2389   /* return flag for lost nonzero bits */
2390   return (j);
2391 }
2392
2393 #else
2394
2395 /* Radix 65536 versions of multiply and divide.  */
2396
2397 /* Multiply significand of e-type number B
2398    by 16-bit quantity A, return e-type result to C.  */
2399
2400 static void
2401 m16m (a, b, c)
2402      unsigned int a;
2403      const UEMUSHORT b[];
2404      UEMUSHORT c[];
2405 {
2406   UEMUSHORT *pp;
2407   unsigned EMULONG carry;
2408   const UEMUSHORT *ps;
2409   UEMUSHORT p[NI];
2410   unsigned EMULONG aa, m;
2411   int i;
2412
2413   aa = a;
2414   pp = &p[NI-2];
2415   *pp++ = 0;
2416   *pp = 0;
2417   ps = &b[NI-1];
2418
2419   for (i=M+1; i<NI; i++)
2420     {
2421       if (*ps == 0)
2422         {
2423           --ps;
2424           --pp;
2425           *(pp-1) = 0;
2426         }
2427       else
2428         {
2429           m = (unsigned EMULONG) aa * *ps--;
2430           carry = (m & 0xffff) + *pp;
2431           *pp-- = (UEMUSHORT) carry;
2432           carry = (carry >> 16) + (m >> 16) + *pp;
2433           *pp = (UEMUSHORT) carry;
2434           *(pp-1) = carry >> 16;
2435         }
2436     }
2437   for (i=M; i<NI; i++)
2438     c[i] = p[i];
2439 }
2440
2441 /* Divide significands of exploded e-types NUM / DEN.  Neither the
2442    numerator NUM nor the denominator DEN is permitted to have its high guard
2443    word nonzero.  */
2444
2445 static int
2446 edivm (den, num)
2447      const UEMUSHORT den[];
2448      UEMUSHORT num[];
2449 {
2450   int i;
2451   UEMUSHORT *p;
2452   unsigned EMULONG tnum;
2453   UEMUSHORT j, tdenm, tquot;
2454   UEMUSHORT tprod[NI+1];
2455
2456   p = &equot[0];
2457   *p++ = num[0];
2458   *p++ = num[1];
2459
2460   for (i=M; i<NI; i++)
2461     {
2462       *p++ = 0;
2463     }
2464   eshdn1 (num);
2465   tdenm = den[M+1];
2466   for (i=M; i<NI; i++)
2467     {
2468       /* Find trial quotient digit (the radix is 65536).  */
2469       tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2470
2471       /* Do not execute the divide instruction if it will overflow.  */
2472       if ((tdenm * (unsigned long) 0xffff) < tnum)
2473         tquot = 0xffff;
2474       else
2475         tquot = tnum / tdenm;
2476       /* Multiply denominator by trial quotient digit.  */
2477       m16m ((unsigned int) tquot, den, tprod);
2478       /* The quotient digit may have been overestimated.  */
2479       if (ecmpm (tprod, num) > 0)
2480         {
2481           tquot -= 1;
2482           esubm (den, tprod);
2483           if (ecmpm (tprod, num) > 0)
2484             {
2485               tquot -= 1;
2486               esubm (den, tprod);
2487             }
2488         }
2489       esubm (tprod, num);
2490       equot[i] = tquot;
2491       eshup6 (num);
2492     }
2493   /* test for nonzero remainder after roundoff bit */
2494   p = &num[M];
2495   j = 0;
2496   for (i=M; i<NI; i++)
2497     {
2498       j |= *p++;
2499     }
2500   if (j)
2501     j = 1;
2502
2503   for (i=0; i<NI; i++)
2504     num[i] = equot[i];
2505
2506   return ((int) j);
2507 }
2508
2509 /* Multiply significands of exploded e-type A and B, result in B.  */
2510
2511 static int
2512 emulm (a, b)
2513      const UEMUSHORT a[];
2514      UEMUSHORT b[];
2515 {
2516   const UEMUSHORT *p;
2517   UEMUSHORT *q;
2518   UEMUSHORT pprod[NI];
2519   UEMUSHORT j;
2520   int i;
2521
2522   equot[0] = b[0];
2523   equot[1] = b[1];
2524   for (i=M; i<NI; i++)
2525     equot[i] = 0;
2526
2527   j = 0;
2528   p = &a[NI-1];
2529   q = &equot[NI-1];
2530   for (i=M+1; i<NI; i++)
2531     {
2532       if (*p == 0)
2533         {
2534           --p;
2535         }
2536       else
2537         {
2538           m16m ((unsigned int) *p--, b, pprod);
2539           eaddm (pprod, equot);
2540         }
2541       j |= *q;
2542       eshdn6 (equot);
2543     }
2544
2545   for (i=0; i<NI; i++)
2546     b[i] = equot[i];
2547
2548   /* return flag for lost nonzero bits */
2549   return ((int) j);
2550 }
2551 #endif
2552
2553
2554 /* Normalize and round off.
2555
2556   The internal format number to be rounded is S.
2557   Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
2558
2559   Input SUBFLG indicates whether the number was obtained
2560   by a subtraction operation.  In that case if LOST is nonzero
2561   then the number is slightly smaller than indicated.
2562
2563   Input EXP is the biased exponent, which may be negative.
2564   the exponent field of S is ignored but is replaced by
2565   EXP as adjusted by normalization and rounding.
2566
2567   Input RCNTRL is the rounding control.  If it is nonzero, the
2568   returned value will be rounded to RNDPRC bits.
2569
2570   For future reference:  In order for emdnorm to round off denormal
2571    significands at the right point, the input exponent must be
2572    adjusted to be the actual value it would have after conversion to
2573    the final floating point type.  This adjustment has been
2574    implemented for all type conversions (etoe53, etc.) and decimal
2575    conversions, but not for the arithmetic functions (eadd, etc.).
2576    Data types having standard 15-bit exponents are not affected by
2577    this, but SFmode and DFmode are affected. For example, ediv with
2578    rndprc = 24 will not round correctly to 24-bit precision if the
2579    result is denormal.  */
2580
2581 static int rlast = -1;
2582 static int rw = 0;
2583 static UEMUSHORT rmsk = 0;
2584 static UEMUSHORT rmbit = 0;
2585 static UEMUSHORT rebit = 0;
2586 static int re = 0;
2587 static UEMUSHORT rbit[NI];
2588
2589 static void
2590 emdnorm (s, lost, subflg, exp, rcntrl)
2591      UEMUSHORT s[];
2592      int lost;
2593      int subflg;
2594      EMULONG exp;
2595      int rcntrl;
2596 {
2597   int i, j;
2598   UEMUSHORT r;
2599
2600   /* Normalize */
2601   j = enormlz (s);
2602
2603   /* a blank significand could mean either zero or infinity.  */
2604 #ifndef INFINITY
2605   if (j > NBITS)
2606     {
2607       ecleazs (s);
2608       return;
2609     }
2610 #endif
2611   exp -= j;
2612 #ifndef INFINITY
2613   if (exp >= 32767L)
2614     goto overf;
2615 #else
2616   if ((j > NBITS) && (exp < 32767))
2617     {
2618       ecleazs (s);
2619       return;
2620     }
2621 #endif
2622   if (exp < 0L)
2623     {
2624       if (exp > (EMULONG) (-NBITS - 1))
2625         {
2626           j = (int) exp;
2627           i = eshift (s, j);
2628           if (i)
2629             lost = 1;
2630         }
2631       else
2632         {
2633           ecleazs (s);
2634           return;
2635         }
2636     }
2637   /* Round off, unless told not to by rcntrl.  */
2638   if (rcntrl == 0)
2639     goto mdfin;
2640   /* Set up rounding parameters if the control register changed.  */
2641   if (rndprc != rlast)
2642     {
2643       ecleaz (rbit);
2644       switch (rndprc)
2645         {
2646         default:
2647         case NBITS:
2648           rw = NI - 1;          /* low guard word */
2649           rmsk = 0xffff;
2650           rmbit = 0x8000;
2651           re = rw - 1;
2652           rebit = 1;
2653           break;
2654
2655         case 113:
2656           rw = 10;
2657           rmsk = 0x7fff;
2658           rmbit = 0x4000;
2659           rebit = 0x8000;
2660           re = rw;
2661           break;
2662
2663         case 64:
2664           rw = 7;
2665           rmsk = 0xffff;
2666           rmbit = 0x8000;
2667           re = rw - 1;
2668           rebit = 1;
2669           break;
2670
2671           /* For DEC or IBM arithmetic */
2672         case 56:
2673           rw = 6;
2674           rmsk = 0xff;
2675           rmbit = 0x80;
2676           rebit = 0x100;
2677           re = rw;
2678           break;
2679
2680         case 53:
2681           rw = 6;
2682           rmsk = 0x7ff;
2683           rmbit = 0x0400;
2684           rebit = 0x800;
2685           re = rw;
2686           break;
2687
2688           /* For C4x arithmetic */
2689         case 32:
2690           rw = 5;
2691           rmsk = 0xffff;
2692           rmbit = 0x8000;
2693           rebit = 1;
2694           re = rw - 1;
2695           break;
2696
2697         case 24:
2698           rw = 4;
2699           rmsk = 0xff;
2700           rmbit = 0x80;
2701           rebit = 0x100;
2702           re = rw;
2703           break;
2704         }
2705       rbit[re] = rebit;
2706       rlast = rndprc;
2707     }
2708
2709   /* Shift down 1 temporarily if the data structure has an implied
2710      most significant bit and the number is denormal.
2711      Intel long double denormals also lose one bit of precision.  */
2712   if ((exp <= 0) && (rndprc != NBITS)
2713       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2714     {
2715       lost |= s[NI - 1] & 1;
2716       eshdn1 (s);
2717     }
2718   /* Clear out all bits below the rounding bit,
2719      remembering in r if any were nonzero.  */
2720   r = s[rw] & rmsk;
2721   if (rndprc < NBITS)
2722     {
2723       i = rw + 1;
2724       while (i < NI)
2725         {
2726           if (s[i])
2727             r |= 1;
2728           s[i] = 0;
2729           ++i;
2730         }
2731     }
2732   s[rw] &= ~rmsk;
2733   if ((r & rmbit) != 0)
2734     {
2735 #ifndef C4X
2736       if (r == rmbit)
2737         {
2738           if (lost == 0)
2739             {                   /* round to even */
2740               if ((s[re] & rebit) == 0)
2741                 goto mddone;
2742             }
2743           else
2744             {
2745               if (subflg != 0)
2746                 goto mddone;
2747             }
2748         }
2749 #endif
2750       eaddm (rbit, s);
2751     }
2752  mddone:
2753 /* Undo the temporary shift for denormal values.  */
2754   if ((exp <= 0) && (rndprc != NBITS)
2755       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2756     {
2757       eshup1 (s);
2758     }
2759   if (s[2] != 0)
2760     {                           /* overflow on roundoff */
2761       eshdn1 (s);
2762       exp += 1;
2763     }
2764  mdfin:
2765   s[NI - 1] = 0;
2766   if (exp >= 32767L)
2767     {
2768 #ifndef INFINITY
2769     overf:
2770 #endif
2771 #ifdef INFINITY
2772       s[1] = 32767;
2773       for (i = 2; i < NI - 1; i++)
2774         s[i] = 0;
2775       if (extra_warnings)
2776         warning ("floating point overflow");
2777 #else
2778       s[1] = 32766;
2779       s[2] = 0;
2780       for (i = M + 1; i < NI - 1; i++)
2781         s[i] = 0xffff;
2782       s[NI - 1] = 0;
2783       if ((rndprc < 64) || (rndprc == 113))
2784         {
2785           s[rw] &= ~rmsk;
2786           if (rndprc == 24)
2787             {
2788               s[5] = 0;
2789               s[6] = 0;
2790             }
2791         }
2792 #endif
2793       return;
2794     }
2795   if (exp < 0)
2796     s[1] = 0;
2797   else
2798     s[1] = (UEMUSHORT) exp;
2799 }
2800
2801 /*  Subtract.  C = B - A, all e type numbers.  */
2802
2803 static int subflg = 0;
2804
2805 static void
2806 esub (a, b, c)
2807      const UEMUSHORT *a, *b;
2808      UEMUSHORT *c;
2809 {
2810
2811 #ifdef NANS
2812   if (eisnan (a))
2813     {
2814       emov (a, c);
2815       return;
2816     }
2817   if (eisnan (b))
2818     {
2819       emov (b, c);
2820       return;
2821     }
2822 /* Infinity minus infinity is a NaN.
2823    Test for subtracting infinities of the same sign.  */
2824   if (eisinf (a) && eisinf (b)
2825       && ((eisneg (a) ^ eisneg (b)) == 0))
2826     {
2827       mtherr ("esub", INVALID);
2828       enan (c, 0);
2829       return;
2830     }
2831 #endif
2832   subflg = 1;
2833   eadd1 (a, b, c);
2834 }
2835
2836 /* Add.  C = A + B, all e type.  */
2837
2838 static void
2839 eadd (a, b, c)
2840      const UEMUSHORT *a, *b;
2841      UEMUSHORT *c;
2842 {
2843
2844 #ifdef NANS
2845 /* NaN plus anything is a NaN.  */
2846   if (eisnan (a))
2847     {
2848       emov (a, c);
2849       return;
2850     }
2851   if (eisnan (b))
2852     {
2853       emov (b, c);
2854       return;
2855     }
2856 /* Infinity minus infinity is a NaN.
2857    Test for adding infinities of opposite signs.  */
2858   if (eisinf (a) && eisinf (b)
2859       && ((eisneg (a) ^ eisneg (b)) != 0))
2860     {
2861       mtherr ("esub", INVALID);
2862       enan (c, 0);
2863       return;
2864     }
2865 #endif
2866   subflg = 0;
2867   eadd1 (a, b, c);
2868 }
2869
2870 /* Arithmetic common to both addition and subtraction.  */
2871
2872 static void
2873 eadd1 (a, b, c)
2874      const UEMUSHORT *a, *b;
2875      UEMUSHORT *c;
2876 {
2877   UEMUSHORT ai[NI], bi[NI], ci[NI];
2878   int i, lost, j, k;
2879   EMULONG lt, lta, ltb;
2880
2881 #ifdef INFINITY
2882   if (eisinf (a))
2883     {
2884       emov (a, c);
2885       if (subflg)
2886         eneg (c);
2887       return;
2888     }
2889   if (eisinf (b))
2890     {
2891       emov (b, c);
2892       return;
2893     }
2894 #endif
2895   emovi (a, ai);
2896   emovi (b, bi);
2897   if (subflg)
2898     ai[0] = ~ai[0];
2899
2900   /* compare exponents */
2901   lta = ai[E];
2902   ltb = bi[E];
2903   lt = lta - ltb;
2904   if (lt > 0L)
2905     {                           /* put the larger number in bi */
2906       emovz (bi, ci);
2907       emovz (ai, bi);
2908       emovz (ci, ai);
2909       ltb = bi[E];
2910       lt = -lt;
2911     }
2912   lost = 0;
2913   if (lt != 0L)
2914     {
2915       if (lt < (EMULONG) (-NBITS - 1))
2916         goto done;              /* answer same as larger addend */
2917       k = (int) lt;
2918       lost = eshift (ai, k);    /* shift the smaller number down */
2919     }
2920   else
2921     {
2922       /* exponents were the same, so must compare significands */
2923       i = ecmpm (ai, bi);
2924       if (i == 0)
2925         {                       /* the numbers are identical in magnitude */
2926           /* if different signs, result is zero */
2927           if (ai[0] != bi[0])
2928             {
2929               eclear (c);
2930               return;
2931             }
2932           /* if same sign, result is double */
2933           /* double denormalized tiny number */
2934           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2935             {
2936               eshup1 (bi);
2937               goto done;
2938             }
2939           /* add 1 to exponent unless both are zero! */
2940           for (j = 1; j < NI - 1; j++)
2941             {
2942               if (bi[j] != 0)
2943                 {
2944                   ltb += 1;
2945                   if (ltb >= 0x7fff)
2946                     {
2947                       eclear (c);
2948                       if (ai[0] != 0)
2949                         eneg (c);
2950                       einfin (c);
2951                       return;
2952                     }
2953                   break;
2954                 }
2955             }
2956           bi[E] = (UEMUSHORT) ltb;
2957           goto done;
2958         }
2959       if (i > 0)
2960         {                       /* put the larger number in bi */
2961           emovz (bi, ci);
2962           emovz (ai, bi);
2963           emovz (ci, ai);
2964         }
2965     }
2966   if (ai[0] == bi[0])
2967     {
2968       eaddm (ai, bi);
2969       subflg = 0;
2970     }
2971   else
2972     {
2973       esubm (ai, bi);
2974       subflg = 1;
2975     }
2976   emdnorm (bi, lost, subflg, ltb, !ROUND_TOWARDS_ZERO);
2977
2978  done:
2979   emovo (bi, c);
2980 }
2981
2982 /* Divide: C = B/A, all e type.  */
2983
2984 static void
2985 ediv (a, b, c)
2986      const UEMUSHORT *a, *b;
2987      UEMUSHORT *c;
2988 {
2989   UEMUSHORT ai[NI], bi[NI];
2990   int i, sign;
2991   EMULONG lt, lta, ltb;
2992
2993 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2994    operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
2995   sign = eisneg (a) ^ eisneg (b);
2996
2997 #ifdef NANS
2998 /* Return any NaN input.  */
2999   if (eisnan (a))
3000     {
3001     emov (a, c);
3002     return;
3003     }
3004   if (eisnan (b))
3005     {
3006     emov (b, c);
3007     return;
3008     }
3009 /* Zero over zero, or infinity over infinity, is a NaN.  */
3010   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
3011       || (eisinf (a) && eisinf (b)))
3012     {
3013     mtherr ("ediv", INVALID);
3014     enan (c, sign);
3015     return;
3016     }
3017 #endif
3018 /* Infinity over anything else is infinity.  */
3019 #ifdef INFINITY
3020   if (eisinf (b))
3021     {
3022       einfin (c);
3023       goto divsign;
3024     }
3025 /* Anything else over infinity is zero.  */
3026   if (eisinf (a))
3027     {
3028       eclear (c);
3029       goto divsign;
3030     }
3031 #endif
3032   emovi (a, ai);
3033   emovi (b, bi);
3034   lta = ai[E];
3035   ltb = bi[E];
3036   if (bi[E] == 0)
3037     {                           /* See if numerator is zero.  */
3038       for (i = 1; i < NI - 1; i++)
3039         {
3040           if (bi[i] != 0)
3041             {
3042               ltb -= enormlz (bi);
3043               goto dnzro1;
3044             }
3045         }
3046       eclear (c);
3047       goto divsign;
3048     }
3049  dnzro1:
3050
3051   if (ai[E] == 0)
3052     {                           /* possible divide by zero */
3053       for (i = 1; i < NI - 1; i++)
3054         {
3055           if (ai[i] != 0)
3056             {
3057               lta -= enormlz (ai);
3058               goto dnzro2;
3059             }
3060         }
3061 /* Divide by zero is not an invalid operation.
3062    It is a divide-by-zero operation!   */
3063       einfin (c);
3064       mtherr ("ediv", SING);
3065       goto divsign;
3066     }
3067  dnzro2:
3068
3069   i = edivm (ai, bi);
3070   /* calculate exponent */
3071   lt = ltb - lta + EXONE;
3072   emdnorm (bi, i, 0, lt, !ROUND_TOWARDS_ZERO);
3073   emovo (bi, c);
3074
3075  divsign:
3076
3077   if (sign
3078 #ifndef IEEE
3079       && (ecmp (c, ezero) != 0)
3080 #endif
3081       )
3082      *(c+(NE-1)) |= 0x8000;
3083   else
3084      *(c+(NE-1)) &= ~0x8000;
3085 }
3086
3087 /* Multiply e-types A and B, return e-type product C.  */
3088
3089 static void
3090 emul (a, b, c)
3091      const UEMUSHORT *a, *b;
3092      UEMUSHORT *c;
3093 {
3094   UEMUSHORT ai[NI], bi[NI];
3095   int i, j, sign;
3096   EMULONG lt, lta, ltb;
3097
3098 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3099    operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
3100   sign = eisneg (a) ^ eisneg (b);
3101
3102 #ifdef NANS
3103 /* NaN times anything is the same NaN.  */
3104   if (eisnan (a))
3105     {
3106     emov (a, c);
3107     return;
3108     }
3109   if (eisnan (b))
3110     {
3111     emov (b, c);
3112     return;
3113     }
3114 /* Zero times infinity is a NaN.  */
3115   if ((eisinf (a) && (ecmp (b, ezero) == 0))
3116       || (eisinf (b) && (ecmp (a, ezero) == 0)))
3117     {
3118     mtherr ("emul", INVALID);
3119     enan (c, sign);
3120     return;
3121     }
3122 #endif
3123 /* Infinity times anything else is infinity.  */
3124 #ifdef INFINITY
3125   if (eisinf (a) || eisinf (b))
3126     {
3127       einfin (c);
3128       goto mulsign;
3129     }
3130 #endif
3131   emovi (a, ai);
3132   emovi (b, bi);
3133   lta = ai[E];
3134   ltb = bi[E];
3135   if (ai[E] == 0)
3136     {
3137       for (i = 1; i < NI - 1; i++)
3138         {
3139           if (ai[i] != 0)
3140             {
3141               lta -= enormlz (ai);
3142               goto mnzer1;
3143             }
3144         }
3145       eclear (c);
3146       goto mulsign;
3147     }
3148  mnzer1:
3149
3150   if (bi[E] == 0)
3151     {
3152       for (i = 1; i < NI - 1; i++)
3153         {
3154           if (bi[i] != 0)
3155             {
3156               ltb -= enormlz (bi);
3157               goto mnzer2;
3158             }
3159         }
3160       eclear (c);
3161       goto mulsign;
3162     }
3163  mnzer2:
3164
3165   /* Multiply significands */
3166   j = emulm (ai, bi);
3167   /* calculate exponent */
3168   lt = lta + ltb - (EXONE - 1);
3169   emdnorm (bi, j, 0, lt, !ROUND_TOWARDS_ZERO);
3170   emovo (bi, c);
3171
3172  mulsign:
3173
3174   if (sign
3175 #ifndef IEEE
3176       && (ecmp (c, ezero) != 0)
3177 #endif
3178       )
3179      *(c+(NE-1)) |= 0x8000;
3180   else
3181      *(c+(NE-1)) &= ~0x8000;
3182 }
3183
3184 /* Convert double precision PE to e-type Y.  */
3185
3186 static void
3187 e53toe (pe, y)
3188      const UEMUSHORT *pe;
3189      UEMUSHORT *y;
3190 {
3191 #ifdef DEC
3192
3193   dectoe (pe, y);
3194
3195 #else
3196 #ifdef IBM
3197
3198   ibmtoe (pe, y, DFmode);
3199
3200 #else
3201 #ifdef C4X
3202
3203   c4xtoe (pe, y, HFmode);
3204
3205 #else
3206
3207   ieeetoe (pe, y, &ieee_53);
3208
3209 #endif /* not C4X */
3210 #endif /* not IBM */
3211 #endif /* not DEC */
3212 }
3213
3214 /* Convert double extended precision float PE to e type Y.  */
3215
3216 static void
3217 e64toe (pe, y)
3218      const UEMUSHORT *pe;
3219      UEMUSHORT *y;
3220 {
3221   UEMUSHORT yy[NI];
3222   const UEMUSHORT *e;
3223   UEMUSHORT *p, *q;
3224   int i;
3225
3226   e = pe;
3227   p = yy;
3228   for (i = 0; i < NE - 5; i++)
3229     *p++ = 0;
3230 #ifndef C4X
3231   /* REAL_WORDS_BIG_ENDIAN is always 0 for DEC and 1 for IBM.
3232      This precision is not ordinarily supported on DEC or IBM.  */
3233   if (! REAL_WORDS_BIG_ENDIAN)
3234     {
3235       for (i = 0; i < 5; i++)
3236         *p++ = *e++;
3237
3238 #ifdef IEEE
3239       /* For denormal long double Intel format, shift significand up one
3240          -- but only if the top significand bit is zero.  A top bit of 1
3241          is "pseudodenormal" when the exponent is zero.  */
3242       if ((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3243         {
3244           UEMUSHORT temp[NI];
3245
3246           emovi (yy, temp);
3247           eshup1 (temp);
3248           emovo (temp,y);
3249           return;
3250         }
3251 #endif /* IEEE */
3252     }
3253   else
3254     {
3255       p = &yy[0] + (NE - 1);
3256       *p-- = *e++;
3257       ++e;
3258       for (i = 0; i < 4; i++)
3259         *p-- = *e++;
3260     }
3261 #endif  /* not C4X */
3262 #ifdef INFINITY
3263   /* Point to the exponent field and check max exponent cases.  */
3264   p = &yy[NE - 1];
3265   if ((*p & 0x7fff) == 0x7fff)
3266     {
3267 #ifdef NANS
3268       if (! REAL_WORDS_BIG_ENDIAN)
3269         {
3270           for (i = 0; i < 4; i++)
3271             {
3272               if ((i != 3 && pe[i] != 0)
3273                   /* Anything but 0x8000 here, including 0, is a NaN.  */
3274                   || (i == 3 && pe[i] != 0x8000))
3275                 {
3276                   enan (y, (*p & 0x8000) != 0);
3277                   return;
3278                 }
3279             }
3280         }
3281       else
3282         {
3283           /* In Motorola extended precision format, the most significant
3284              bit of an infinity mantissa could be either 1 or 0.  It is
3285              the lower order bits that tell whether the value is a NaN.  */
3286           if ((pe[2] & 0x7fff) != 0)
3287             goto bigend_nan;
3288
3289           for (i = 3; i <= 5; i++)
3290             {
3291               if (pe[i] != 0)
3292                 {
3293 bigend_nan:
3294                   enan (y, (*p & 0x8000) != 0);
3295                   return;
3296                 }
3297             }
3298         }
3299 #endif /* NANS */
3300       eclear (y);
3301       einfin (y);
3302       if (*p & 0x8000)
3303         eneg (y);
3304       return;
3305     }
3306 #endif  /* INFINITY */
3307   p = yy;
3308   q = y;
3309   for (i = 0; i < NE; i++)
3310     *q++ = *p++;
3311 }
3312
3313 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3314 /* Convert 128-bit long double precision float PE to e type Y.  */
3315
3316 static void
3317 e113toe (pe, y)
3318      const UEMUSHORT *pe;
3319      UEMUSHORT *y;
3320 {
3321   ieeetoe (pe, y, &ieee_113);
3322 }
3323 #endif  /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3324
3325 /* Convert single precision float PE to e type Y.  */
3326
3327 static void
3328 e24toe (pe, y)
3329      const UEMUSHORT *pe;
3330      UEMUSHORT *y;
3331 {
3332 #ifdef IBM
3333
3334   ibmtoe (pe, y, SFmode);
3335
3336 #else
3337
3338 #ifdef C4X
3339
3340   c4xtoe (pe, y, QFmode);
3341
3342 #else
3343 #ifdef DEC
3344
3345   ieeetoe (pe, y, &dec_f);
3346
3347 #else
3348
3349   ieeetoe (pe, y, &ieee_24);
3350
3351 #endif /* not DEC */
3352 #endif /* not C4X */
3353 #endif /* not IBM */
3354 }
3355
3356 /* Convert machine format float of specified format PE to e type Y.  */
3357
3358 static void
3359 ieeetoe (pe, y, fmt)
3360      const UEMUSHORT *pe;
3361      UEMUSHORT *y;
3362      const struct ieee_format *fmt;
3363 {
3364   UEMUSHORT r;
3365   const UEMUSHORT *e;
3366   UEMUSHORT *p;
3367   UEMUSHORT yy[NI];
3368   int denorm, i, k;
3369   int shortsm1 = fmt->bits / 16 - 1;
3370 #ifdef INFINITY
3371   int expmask = (1 << fmt->expbits) - 1;
3372 #endif
3373 int expshift = (fmt->precision - 1) & 0x0f;
3374   int highbit = 1 << expshift;
3375
3376   e = pe;
3377   denorm = 0;
3378   ecleaz (yy);
3379   if (! REAL_WORDS_BIG_ENDIAN)
3380     e += shortsm1;
3381   r = *e;
3382   yy[0] = 0;
3383   if (r & 0x8000)
3384     yy[0] = 0xffff;
3385   yy[M] = (r & (highbit - 1)) | highbit;
3386   r = (r & 0x7fff) >> expshift;
3387 #ifdef INFINITY
3388   if (!LARGEST_EXPONENT_IS_NORMAL (fmt->precision) && r == expmask)
3389     {
3390 #ifdef NANS
3391       /* First check the word where high order mantissa and exponent live */
3392       if ((*e & (highbit - 1)) != 0)
3393         {
3394           enan (y, yy[0] != 0);
3395           return;
3396         }
3397       if (! REAL_WORDS_BIG_ENDIAN)
3398         {
3399           for (i = 0; i < shortsm1; i++)
3400             {
3401               if (pe[i] != 0)
3402                 {
3403                   enan (y, yy[0] != 0);
3404                   return;
3405                 }
3406             }
3407         }
3408       else
3409         {
3410           for (i = 1; i < shortsm1 + 1; i++)
3411             {
3412               if (pe[i] != 0)
3413                 {
3414                   enan (y, yy[0] != 0);
3415                   return;
3416                 }
3417             }
3418         }
3419 #endif /* NANS */
3420       eclear (y);
3421       einfin (y);
3422       if (yy[0])
3423         eneg (y);
3424       return;
3425     }
3426 #endif  /* INFINITY */
3427   /* If zero exponent, then the significand is denormalized.
3428      So take back the understood high significand bit.  */
3429   if (r == 0)
3430     {
3431       denorm = 1;
3432       yy[M] &= ~highbit;
3433     }
3434   r += fmt->adjustment;
3435   yy[E] = r;
3436   p = &yy[M + 1];
3437   if (! REAL_WORDS_BIG_ENDIAN)
3438     {
3439       for (i = 0; i < shortsm1; i++)
3440         *p++ = *(--e);
3441     }
3442   else
3443     {
3444       ++e;
3445       for (i = 0; i < shortsm1; i++)
3446         *p++ = *e++;
3447     }
3448   if (fmt->precision == 113)
3449     {
3450       /* denorm is left alone in 113 bit format */
3451       if (!denorm)
3452         eshift (yy, -1);
3453     }
3454   else
3455     {
3456       eshift (yy, -(expshift + 1));
3457       if (denorm)
3458         { /* if zero exponent, then normalize the significand */
3459           if ((k = enormlz (yy)) > NBITS)
3460             ecleazs (yy);
3461           else
3462             yy[E] -= (UEMUSHORT) (k - 1);
3463         }
3464     }
3465   emovo (yy, y);
3466 }
3467
3468 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3469 /* Convert e-type X to IEEE 128-bit long double format E.  */
3470
3471 static void
3472 etoe113 (x, e)
3473      const UEMUSHORT *x;
3474      UEMUSHORT *e;
3475 {
3476   etoieee (x, e, &ieee_113);
3477 }
3478
3479 /* Convert exploded e-type X, that has already been rounded to
3480    113-bit precision, to IEEE 128-bit long double format Y.  */
3481
3482 static void
3483 toe113 (x, y)
3484      UEMUSHORT *x, *y;
3485 {
3486   toieee (x, y, &ieee_113);
3487 }
3488
3489 #endif  /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3490
3491 /* Convert e-type X to IEEE double extended format E.  */
3492
3493 static void
3494 etoe64 (x, e)
3495      const UEMUSHORT *x;
3496      UEMUSHORT *e;
3497 {
3498   etoieee (x, e, &ieee_64);
3499 }
3500
3501 /* Convert exploded e-type X, that has already been rounded to
3502    64-bit precision, to IEEE double extended format Y.  */
3503
3504 static void
3505 toe64 (x, y)
3506      UEMUSHORT *x, *y;
3507 {
3508   toieee (x, y, &ieee_64);
3509 }
3510
3511 /* e type to double precision.  */
3512
3513 #ifdef DEC
3514 /* Convert e-type X to DEC-format double E.  */
3515
3516 static void
3517 etoe53 (x, e)
3518      const UEMUSHORT *x;
3519      UEMUSHORT *e;
3520 {
3521   etodec (x, e);
3522 }
3523
3524 /* Convert exploded e-type X, that has already been rounded to
3525    56-bit double precision, to DEC double Y.  */
3526
3527 static void
3528 toe53 (x, y)
3529      UEMUSHORT *x, *y;
3530 {
3531   todec (x, y);
3532 }
3533
3534 #else
3535 #ifdef IBM
3536 /* Convert e-type X to IBM 370-format double E.  */
3537
3538 static void
3539 etoe53 (x, e)
3540      const UEMUSHORT *x;
3541      UEMUSHORT *e;
3542 {
3543   etoibm (x, e, DFmode);
3544 }
3545
3546 /* Convert exploded e-type X, that has already been rounded to
3547    56-bit precision, to IBM 370 double Y.  */
3548
3549 static void
3550 toe53 (x, y)
3551      UEMUSHORT *x, *y;
3552 {
3553   toibm (x, y, DFmode);
3554 }
3555
3556 #else /* it's neither DEC nor IBM */
3557 #ifdef C4X
3558 /* Convert e-type X to C4X-format long double E.  */
3559
3560 static void
3561 etoe53 (x, e)
3562      const UEMUSHORT *x;
3563      UEMUSHORT *e;
3564 {
3565   etoc4x (x, e, HFmode);
3566 }
3567
3568 /* Convert exploded e-type X, that has already been rounded to
3569    56-bit precision, to IBM 370 double Y.  */
3570
3571 static void
3572 toe53 (x, y)
3573      UEMUSHORT *x, *y;
3574 {
3575   toc4x (x, y, HFmode);
3576 }
3577
3578 #else  /* it's neither DEC nor IBM nor C4X */
3579
3580 /* Convert e-type X to IEEE double E.  */
3581
3582 static void
3583 etoe53 (x, e)
3584       const UEMUSHORT *x;
3585       UEMUSHORT *e;
3586 {
3587   etoieee (x, e, &ieee_53);
3588 }
3589
3590 /* Convert exploded e-type X, that has already been rounded to
3591    53-bit precision, to IEEE double Y.  */
3592
3593 static void
3594 toe53 (x, y)
3595      UEMUSHORT *x, *y;
3596 {
3597   toieee (x, y, &ieee_53);
3598 }
3599
3600 #endif /* not C4X */
3601 #endif /* not IBM */
3602 #endif /* not DEC */
3603
3604
3605
3606 /* e type to single precision.  */
3607
3608 #ifdef IBM
3609 /* Convert e-type X to IBM 370 float E.  */
3610
3611 static void
3612 etoe24 (x, e)
3613      const UEMUSHORT *x;
3614      UEMUSHORT *e;
3615 {
3616   etoibm (x, e, SFmode);
3617 }
3618
3619 /* Convert exploded e-type X, that has already been rounded to
3620    float precision, to IBM 370 float Y.  */
3621
3622 static void
3623 toe24 (x, y)
3624      UEMUSHORT *x, *y;
3625 {
3626   toibm (x, y, SFmode);
3627 }
3628
3629 #else /* it's not IBM */
3630
3631 #ifdef C4X
3632 /* Convert e-type X to C4X float E.  */
3633
3634 static void
3635 etoe24 (x, e)
3636      const UEMUSHORT *x;
3637      UEMUSHORT *e;
3638 {
3639   etoc4x (x, e, QFmode);
3640 }
3641
3642 /* Convert exploded e-type X, that has already been rounded to
3643    float precision, to IBM 370 float Y.  */
3644
3645 static void
3646 toe24 (x, y)
3647      UEMUSHORT *x, *y;
3648 {
3649   toc4x (x, y, QFmode);
3650 }
3651
3652 #else /* it's neither IBM nor C4X */
3653
3654 #ifdef DEC
3655
3656 /* Convert e-type X to DEC F-float E.  */
3657
3658 static void
3659 etoe24 (x, e)
3660       const UEMUSHORT *x;
3661       UEMUSHORT *e;
3662 {
3663   etoieee (x, e, &dec_f);
3664 }
3665
3666 /* Convert exploded e-type X, that has already been rounded to
3667    float precision, to DEC F-float Y.  */
3668
3669 static void
3670 toe24 (x, y)
3671       UEMUSHORT *x, *y;
3672 {
3673   toieee (x, y, &dec_f);
3674 }
3675
3676 #else
3677
3678 /* Convert e-type X to IEEE float E.  */
3679
3680 static void
3681 etoe24 (x, e)
3682       const UEMUSHORT *x;
3683       UEMUSHORT *e;
3684 {
3685   etoieee (x, e, &ieee_24);
3686 }
3687
3688 /* Convert exploded e-type X, that has already been rounded to
3689    float precision, to IEEE float Y.  */
3690
3691 static void
3692 toe24 (x, y)
3693       UEMUSHORT *x, *y;
3694 {
3695   toieee (x, y, &ieee_24);
3696 }
3697
3698 #endif  /* not DEC */
3699 #endif  /* not C4X */
3700 #endif  /* not IBM */
3701
3702
3703 /* Convert e-type X to the IEEE format described by FMT.  */
3704
3705 static void
3706 etoieee (x, e, fmt)
3707       const UEMUSHORT *x;
3708       UEMUSHORT *e;
3709       const struct ieee_format *fmt;
3710 {
3711   UEMUSHORT xi[NI];
3712   EMULONG exp;
3713   int rndsav;
3714
3715 #ifdef NANS
3716   if (eisnan (x))
3717     {
3718       make_nan (e, eisneg (x), fmt->mode);
3719       return;
3720     }
3721 #endif
3722
3723   emovi (x, xi);
3724
3725 #ifdef INFINITY
3726   if (eisinf (x))
3727     goto nonorm;
3728 #endif
3729   /* Adjust exponent for offset.  */
3730   exp = (EMULONG) xi[E] - fmt->adjustment;
3731
3732   /* Round off to nearest or even.  */
3733   rndsav = rndprc;
3734   rndprc = fmt->precision;
3735   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3736   rndprc = rndsav;
3737 #ifdef INFINITY
3738  nonorm:
3739 #endif
3740   toieee (xi, e, fmt);
3741 }
3742
3743 /* Convert exploded e-type X, that has already been rounded to
3744    the necessary precision, to the IEEE format described by FMT.  */
3745
3746 static void
3747 toieee (x, y, fmt)
3748       UEMUSHORT *x, *y;
3749       const struct ieee_format *fmt;
3750 {
3751   UEMUSHORT maxexp;
3752   UEMUSHORT *q;
3753   int words;
3754   int i;
3755
3756   maxexp = (1 << fmt->expbits) - 1;
3757   words = (fmt->bits - fmt->expbits) / EMUSHORT_SIZE;
3758
3759 #ifdef NANS
3760   if (eiisnan (x))
3761     {
3762       make_nan (y, eiisneg (x), fmt->mode);
3763       return;
3764     }
3765 #endif
3766
3767   if (fmt->expbits < 15
3768       && LARGEST_EXPONENT_IS_NORMAL (fmt->bits)
3769       && x[E] > maxexp)
3770     {
3771       saturate (y, eiisneg (x), fmt->bits, 1);
3772       return;
3773     }
3774
3775   /* Point to the exponent.  */
3776   if (REAL_WORDS_BIG_ENDIAN)
3777     q = y;
3778   else
3779     q = y + words;
3780
3781   /* Copy the sign.  */
3782   if (x[0])
3783     *q = 0x8000;
3784   else
3785     *q = 0;
3786
3787   if (fmt->expbits < 15
3788       && !LARGEST_EXPONENT_IS_NORMAL (fmt->bits)
3789       && x[E] >= maxexp)
3790     {
3791       /* Saturate at largest number less that infinity.  */
3792       UEMUSHORT fill;
3793 #ifdef INFINITY
3794       *q |= maxexp << (15 - fmt->expbits);
3795       fill = 0;
3796 #else
3797       *q |= (maxexp << (15 - fmt->expbits)) - 1;
3798       fill = 0xffff;
3799 #endif
3800
3801       if (!REAL_WORDS_BIG_ENDIAN)
3802         {
3803           for (i = 0; i < words; i++)
3804             *(--q) = fill;
3805         }
3806       else
3807         {
3808           for (i = 0; i < words; i++)
3809             *(++q) = fill;
3810         }
3811 #if defined(INFINITY) && defined(ERANGE)
3812       errno = ERANGE;
3813 #endif
3814       return;
3815     }
3816
3817   /* If denormal and DEC float, return zero (DEC has no denormals) */
3818 #ifdef DEC
3819   if (x[E] == 0)
3820     {
3821       for (i = 0; i < fmt->bits / EMUSHORT_SIZE ; i++)
3822         q[i] = 0;
3823       return;
3824     }
3825 #endif /* DEC */
3826
3827   /* Delete the implied bit unless denormal, except for
3828      64-bit precision.  */
3829   if (fmt->precision != 64 && x[E] != 0)
3830     {
3831       eshup1 (x);
3832     }
3833
3834   /* Shift denormal double extended Intel format significand down
3835      one bit. */
3836   if (fmt->precision == 64 && x[E] == 0 && ! REAL_WORDS_BIG_ENDIAN)
3837     eshdn1 (x);
3838
3839   if (fmt->expbits < 15)
3840     {
3841       /* Shift the significand.  */
3842       eshift (x, 15 - fmt->expbits);
3843
3844       /* Combine the exponent and upper bits of the significand.  */
3845       *q |= x[E] << (15 - fmt->expbits);
3846       *q |= x[M] & (UEMUSHORT) ~((maxexp << (15 - fmt->expbits)) | 0x8000);
3847     }
3848   else
3849     {
3850       /* Copy the exponent.  */
3851       *q |= x[E];
3852     }
3853
3854   /* Add padding after the exponent. At the moment, this is only necessary for
3855      64-bit precision; in this case, the padding is 16 bits.  */
3856   if (fmt->precision == 64)
3857     {
3858       *(q + 1) = 0;
3859
3860       /* Skip padding.  */
3861       if (REAL_WORDS_BIG_ENDIAN)
3862         ++q;
3863     }
3864
3865   /* Copy the significand.  */
3866   if (REAL_WORDS_BIG_ENDIAN)
3867     {
3868       for (i = 0; i < words; i++)
3869         *(++q) = x[i + M + 1];
3870     }
3871 #ifdef INFINITY
3872   else if (fmt->precision == 64 && eiisinf (x))
3873     {
3874       /* Intel double extended infinity significand.  */
3875       *(--q) = 0x8000;
3876       *(--q) = 0;
3877       *(--q) = 0;
3878       *(--q) = 0;
3879     }
3880 #endif
3881   else
3882     {
3883       for (i = 0; i < words; i++)
3884         *(--q) = x[i + M + 1];
3885     }
3886 }
3887
3888
3889 /* Compare two e type numbers.
3890    Return +1 if a > b
3891            0 if a == b
3892           -1 if a < b
3893           -2 if either a or b is a NaN.  */
3894
3895 static int
3896 ecmp (a, b)
3897      const UEMUSHORT *a, *b;
3898 {
3899   UEMUSHORT ai[NI], bi[NI];
3900   UEMUSHORT *p, *q;
3901   int i;
3902   int msign;
3903
3904 #ifdef NANS
3905   if (eisnan (a)  || eisnan (b))
3906       return (-2);
3907 #endif
3908   emovi (a, ai);
3909   p = ai;
3910   emovi (b, bi);
3911   q = bi;
3912
3913   if (*p != *q)
3914     {                           /* the signs are different */
3915       /* -0 equals + 0 */
3916       for (i = 1; i < NI - 1; i++)
3917         {
3918           if (ai[i] != 0)
3919             goto nzro;
3920           if (bi[i] != 0)
3921             goto nzro;
3922         }
3923       return (0);
3924     nzro:
3925       if (*p == 0)
3926         return (1);
3927       else
3928         return (-1);
3929     }
3930   /* both are the same sign */
3931   if (*p == 0)
3932     msign = 1;
3933   else
3934     msign = -1;
3935   i = NI - 1;
3936   do
3937     {
3938       if (*p++ != *q++)
3939         {
3940           goto diff;
3941         }
3942     }
3943   while (--i > 0);
3944
3945   return (0);                   /* equality */
3946
3947  diff:
3948
3949   if (*(--p) > *(--q))
3950     return (msign);             /* p is bigger */
3951   else
3952     return (-msign);            /* p is littler */
3953 }
3954
3955 #if 0
3956 /* Find e-type nearest integer to X, as floor (X + 0.5).  */
3957
3958 static void
3959 eround (x, y)
3960      const UEMUSHORT *x;
3961      UEMUSHORT *y;
3962 {
3963   eadd (ehalf, x, y);
3964   efloor (y, y);
3965 }
3966 #endif /* 0 */
3967
3968 /* Convert HOST_WIDE_INT LP to e type Y.  */
3969
3970 static void
3971 ltoe (lp, y)
3972      const HOST_WIDE_INT *lp;
3973      UEMUSHORT *y;
3974 {
3975   UEMUSHORT yi[NI];
3976   unsigned HOST_WIDE_INT ll;
3977   int k;
3978
3979   ecleaz (yi);
3980   if (*lp < 0)
3981     {
3982       /* make it positive */
3983       ll = (unsigned HOST_WIDE_INT) (-(*lp));
3984       yi[0] = 0xffff;           /* put correct sign in the e type number */
3985     }
3986   else
3987     {
3988       ll = (unsigned HOST_WIDE_INT) (*lp);
3989     }
3990   /* move the long integer to yi significand area */
3991 #if HOST_BITS_PER_WIDE_INT == 64
3992   yi[M] = (UEMUSHORT) (ll >> 48);
3993   yi[M + 1] = (UEMUSHORT) (ll >> 32);
3994   yi[M + 2] = (UEMUSHORT) (ll >> 16);
3995   yi[M + 3] = (UEMUSHORT) ll;
3996   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
3997 #else
3998   yi[M] = (UEMUSHORT) (ll >> 16);
3999   yi[M + 1] = (UEMUSHORT) ll;
4000   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4001 #endif
4002
4003   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4004     ecleaz (yi);                /* it was zero */
4005   else
4006     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
4007   emovo (yi, y);                /* output the answer */
4008 }
4009
4010 /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
4011
4012 static void
4013 ultoe (lp, y)
4014      const unsigned HOST_WIDE_INT *lp;
4015      UEMUSHORT *y;
4016 {
4017   UEMUSHORT yi[NI];
4018   unsigned HOST_WIDE_INT ll;
4019   int k;
4020
4021   ecleaz (yi);
4022   ll = *lp;
4023
4024   /* move the long integer to ayi significand area */
4025 #if HOST_BITS_PER_WIDE_INT == 64
4026   yi[M] = (UEMUSHORT) (ll >> 48);
4027   yi[M + 1] = (UEMUSHORT) (ll >> 32);
4028   yi[M + 2] = (UEMUSHORT) (ll >> 16);
4029   yi[M + 3] = (UEMUSHORT) ll;
4030   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
4031 #else
4032   yi[M] = (UEMUSHORT) (ll >> 16);
4033   yi[M + 1] = (UEMUSHORT) ll;
4034   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4035 #endif
4036
4037   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4038     ecleaz (yi);                /* it was zero */
4039   else
4040     yi[E] -= (UEMUSHORT) k;  /* subtract shift count from exponent */
4041   emovo (yi, y);                /* output the answer */
4042 }
4043
4044
4045 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4046    part FRAC of e-type (packed internal format) floating point input X.
4047    The integer output I has the sign of the input, except that
4048    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4049    The output e-type fraction FRAC is the positive fractional
4050    part of abs (X).  */
4051
4052 static void
4053 eifrac (x, i, frac)
4054      const UEMUSHORT *x;
4055      HOST_WIDE_INT *i;
4056      UEMUSHORT *frac;
4057 {
4058   UEMUSHORT xi[NI];
4059   int j, k;
4060   unsigned HOST_WIDE_INT ll;
4061
4062   emovi (x, xi);
4063   k = (int) xi[E] - (EXONE - 1);
4064   if (k <= 0)
4065     {
4066       /* if exponent <= 0, integer = 0 and real output is fraction */
4067       *i = 0L;
4068       emovo (xi, frac);
4069       return;
4070     }
4071   if (k > (HOST_BITS_PER_WIDE_INT - 1))
4072     {
4073       /* long integer overflow: output large integer
4074          and correct fraction  */
4075       if (xi[0])
4076         *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4077       else
4078         {
4079 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4080           /* In this case, let it overflow and convert as if unsigned.  */
4081           euifrac (x, &ll, frac);
4082           *i = (HOST_WIDE_INT) ll;
4083           return;
4084 #else
4085           /* In other cases, return the largest positive integer.  */
4086           *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4087 #endif
4088         }
4089       eshift (xi, k);
4090       if (extra_warnings)
4091         warning ("overflow on truncation to integer");
4092     }
4093   else if (k > 16)
4094     {
4095       /* Shift more than 16 bits: first shift up k-16 mod 16,
4096          then shift up by 16's.  */
4097       j = k - ((k >> 4) << 4);
4098       eshift (xi, j);
4099       ll = xi[M];
4100       k -= j;
4101       do
4102         {
4103           eshup6 (xi);
4104           ll = (ll << 16) | xi[M];
4105         }
4106       while ((k -= 16) > 0);
4107       *i = ll;
4108       if (xi[0])
4109         *i = -(*i);
4110     }
4111   else
4112     {
4113       /* shift not more than 16 bits */
4114       eshift (xi, k);
4115       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4116       if (xi[0])
4117         *i = -(*i);
4118     }
4119   xi[0] = 0;
4120   xi[E] = EXONE - 1;
4121   xi[M] = 0;
4122   if ((k = enormlz (xi)) > NBITS)
4123     ecleaz (xi);
4124   else
4125     xi[E] -= (UEMUSHORT) k;
4126
4127   emovo (xi, frac);
4128 }
4129
4130
4131 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4132    FRAC of e-type X.  A negative input yields integer output = 0 but
4133    correct fraction.  */
4134
4135 static void
4136 euifrac (x, i, frac)
4137      const UEMUSHORT *x;
4138      unsigned HOST_WIDE_INT *i;
4139      UEMUSHORT *frac;
4140 {
4141   unsigned HOST_WIDE_INT ll;
4142   UEMUSHORT xi[NI];
4143   int j, k;
4144
4145   emovi (x, xi);
4146   k = (int) xi[E] - (EXONE - 1);
4147   if (k <= 0)
4148     {
4149       /* if exponent <= 0, integer = 0 and argument is fraction */
4150       *i = 0L;
4151       emovo (xi, frac);
4152       return;
4153     }
4154   if (k > HOST_BITS_PER_WIDE_INT)
4155     {
4156       /* Long integer overflow: output large integer
4157          and correct fraction.
4158          Note, the BSD MicroVAX compiler says that ~(0UL)
4159          is a syntax error.  */
4160       *i = ~(0L);
4161       eshift (xi, k);
4162       if (extra_warnings)
4163         warning ("overflow on truncation to unsigned integer");
4164     }
4165   else if (k > 16)
4166     {
4167       /* Shift more than 16 bits: first shift up k-16 mod 16,
4168          then shift up by 16's.  */
4169       j = k - ((k >> 4) << 4);
4170       eshift (xi, j);
4171       ll = xi[M];
4172       k -= j;
4173       do
4174         {
4175           eshup6 (xi);
4176           ll = (ll << 16) | xi[M];
4177         }
4178       while ((k -= 16) > 0);
4179       *i = ll;
4180     }
4181   else
4182     {
4183       /* shift not more than 16 bits */
4184       eshift (xi, k);
4185       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4186     }
4187
4188   if (xi[0])  /* A negative value yields unsigned integer 0.  */
4189     *i = 0L;
4190
4191   xi[0] = 0;
4192   xi[E] = EXONE - 1;
4193   xi[M] = 0;
4194   if ((k = enormlz (xi)) > NBITS)
4195     ecleaz (xi);
4196   else
4197     xi[E] -= (UEMUSHORT) k;
4198
4199   emovo (xi, frac);
4200 }
4201
4202 /* Shift the significand of exploded e-type X up or down by SC bits.  */
4203
4204 static int
4205 eshift (x, sc)
4206      UEMUSHORT *x;
4207      int sc;
4208 {
4209   UEMUSHORT lost;
4210   UEMUSHORT *p;
4211
4212   if (sc == 0)
4213     return (0);
4214
4215   lost = 0;
4216   p = x + NI - 1;
4217
4218   if (sc < 0)
4219     {
4220       sc = -sc;
4221       while (sc >= 16)
4222         {
4223           lost |= *p;           /* remember lost bits */
4224           eshdn6 (x);
4225           sc -= 16;
4226         }
4227
4228       while (sc >= 8)
4229         {
4230           lost |= *p & 0xff;
4231           eshdn8 (x);
4232           sc -= 8;
4233         }
4234
4235       while (sc > 0)
4236         {
4237           lost |= *p & 1;
4238           eshdn1 (x);
4239           sc -= 1;
4240         }
4241     }
4242   else
4243     {
4244       while (sc >= 16)
4245         {
4246           eshup6 (x);
4247           sc -= 16;
4248         }
4249
4250       while (sc >= 8)
4251         {
4252           eshup8 (x);
4253           sc -= 8;
4254         }
4255
4256       while (sc > 0)
4257         {
4258           eshup1 (x);
4259           sc -= 1;
4260         }
4261     }
4262   if (lost)
4263     lost = 1;
4264   return ((int) lost);
4265 }
4266
4267 /* Shift normalize the significand area of exploded e-type X.
4268    Return the shift count (up = positive).  */
4269
4270 static int
4271 enormlz (x)
4272      UEMUSHORT x[];
4273 {
4274   UEMUSHORT *p;
4275   int sc;
4276
4277   sc = 0;
4278   p = &x[M];
4279   if (*p != 0)
4280     goto normdn;
4281   ++p;
4282   if (*p & 0x8000)
4283     return (0);                 /* already normalized */
4284   while (*p == 0)
4285     {
4286       eshup6 (x);
4287       sc += 16;
4288
4289       /* With guard word, there are NBITS+16 bits available.
4290        Return true if all are zero.  */
4291       if (sc > NBITS)
4292         return (sc);
4293     }
4294   /* see if high byte is zero */
4295   while ((*p & 0xff00) == 0)
4296     {
4297       eshup8 (x);
4298       sc += 8;
4299     }
4300   /* now shift 1 bit at a time */
4301   while ((*p & 0x8000) == 0)
4302     {
4303       eshup1 (x);
4304       sc += 1;
4305       if (sc > NBITS)
4306         {
4307           mtherr ("enormlz", UNDERFLOW);
4308           return (sc);
4309         }
4310     }
4311   return (sc);
4312
4313   /* Normalize by shifting down out of the high guard word
4314      of the significand */
4315  normdn:
4316
4317   if (*p & 0xff00)
4318     {
4319       eshdn8 (x);
4320       sc -= 8;
4321     }
4322   while (*p != 0)
4323     {
4324       eshdn1 (x);
4325       sc -= 1;
4326
4327       if (sc < -NBITS)
4328         {
4329           mtherr ("enormlz", OVERFLOW);
4330           return (sc);
4331         }
4332     }
4333   return (sc);
4334 }
4335
4336 /* Powers of ten used in decimal <-> binary conversions.  */
4337
4338 #define NTEN 12
4339 #define MAXP 4096
4340
4341 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4342 static const UEMUSHORT etens[NTEN + 1][NE] =
4343 {
4344   {0x6576, 0x4a92, 0x804a, 0x153f,
4345    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4346   {0x6a32, 0xce52, 0x329a, 0x28ce,
4347    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4348   {0x526c, 0x50ce, 0xf18b, 0x3d28,
4349    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4350   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4351    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4352   {0x851e, 0xeab7, 0x98fe, 0x901b,
4353    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4354   {0x0235, 0x0137, 0x36b1, 0x336c,
4355    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4356   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4357    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4358   {0x0000, 0x0000, 0x0000, 0x0000,
4359    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4360   {0x0000, 0x0000, 0x0000, 0x0000,
4361    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4362   {0x0000, 0x0000, 0x0000, 0x0000,
4363    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4364   {0x0000, 0x0000, 0x0000, 0x0000,
4365    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4366   {0x0000, 0x0000, 0x0000, 0x0000,
4367    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4368   {0x0000, 0x0000, 0x0000, 0x0000,
4369    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4370 };
4371
4372 static const UEMUSHORT emtens[NTEN + 1][NE] =
4373 {
4374   {0x2030, 0xcffc, 0xa1c3, 0x8123,
4375    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4376   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4377    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4378   {0xf53f, 0xf698, 0x6bd3, 0x0158,
4379    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4380   {0xe731, 0x04d4, 0xe3f2, 0xd332,
4381    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4382   {0xa23e, 0x5308, 0xfefb, 0x1155,
4383    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4384   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4385    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4386   {0x2a20, 0x6224, 0x47b3, 0x98d7,
4387    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4388   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4389    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4390   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4391    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4392   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4393    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4394   {0xc155, 0xa4a8, 0x404e, 0x6113,
4395    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4396   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4397    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4398   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4399    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4400 };
4401 #else
4402 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4403 static const UEMUSHORT etens[NTEN + 1][NE] =
4404 {
4405   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4406   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4407   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4408   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4409   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4410   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4411   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4412   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4413   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4414   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4415   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4416   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4417   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4418 };
4419
4420 static const UEMUSHORT emtens[NTEN + 1][NE] =
4421 {
4422   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4423   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4424   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4425   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4426   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4427   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4428   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4429   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4430   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4431   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4432   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4433   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4434   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4435 };
4436 #endif
4437
4438 #if 0
4439 /* Convert float value X to ASCII string STRING with NDIG digits after
4440    the decimal point.  */
4441
4442 static void
4443 e24toasc (x, string, ndigs)
4444      const UEMUSHORT x[];
4445      char *string;
4446      int ndigs;
4447 {
4448   UEMUSHORT w[NI];
4449
4450   e24toe (x, w);
4451   etoasc (w, string, ndigs);
4452 }
4453
4454 /* Convert double value X to ASCII string STRING with NDIG digits after
4455    the decimal point.  */
4456
4457 static void
4458 e53toasc (x, string, ndigs)
4459      const UEMUSHORT x[];
4460      char *string;
4461      int ndigs;
4462 {
4463   UEMUSHORT w[NI];
4464
4465   e53toe (x, w);
4466   etoasc (w, string, ndigs);
4467 }
4468
4469 /* Convert double extended value X to ASCII string STRING with NDIG digits
4470    after the decimal point.  */
4471
4472 static void
4473 e64toasc (x, string, ndigs)
4474      const UEMUSHORT x[];
4475      char *string;
4476      int ndigs;
4477 {
4478   UEMUSHORT w[NI];
4479
4480   e64toe (x, w);
4481   etoasc (w, string, ndigs);
4482 }
4483
4484 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4485    after the decimal point.  */
4486
4487 static void
4488 e113toasc (x, string, ndigs)
4489      const UEMUSHORT x[];
4490      char *string;
4491      int ndigs;
4492 {
4493   UEMUSHORT w[NI];
4494
4495   e113toe (x, w);
4496   etoasc (w, string, ndigs);
4497 }
4498 #endif /* 0 */
4499
4500 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4501    the decimal point.  */
4502
4503 static char wstring[80];        /* working storage for ASCII output */
4504
4505 static void
4506 etoasc (x, string, ndigs)
4507      const UEMUSHORT x[];
4508      char *string;
4509      int ndigs;
4510 {
4511   EMUSHORT digit;
4512   UEMUSHORT y[NI], t[NI], u[NI], w[NI];
4513   const UEMUSHORT *p, *r, *ten;
4514   UEMUSHORT sign;
4515   int i, j, k, expon, rndsav;
4516   char *s, *ss;
4517   UEMUSHORT m;
4518
4519
4520   rndsav = rndprc;
4521   ss = string;
4522   s = wstring;
4523   *ss = '\0';
4524   *s = '\0';
4525 #ifdef NANS
4526   if (eisnan (x))
4527     {
4528       sprintf (wstring, " NaN ");
4529       goto bxit;
4530     }
4531 #endif
4532   rndprc = NBITS;               /* set to full precision */
4533   emov (x, y);                  /* retain external format */
4534   if (y[NE - 1] & 0x8000)
4535     {
4536       sign = 0xffff;
4537       y[NE - 1] &= 0x7fff;
4538     }
4539   else
4540     {
4541       sign = 0;
4542     }
4543   expon = 0;
4544   ten = &etens[NTEN][0];
4545   emov (eone, t);
4546   /* Test for zero exponent */
4547   if (y[NE - 1] == 0)
4548     {
4549       for (k = 0; k < NE - 1; k++)
4550         {
4551           if (y[k] != 0)
4552             goto tnzro;         /* denormalized number */
4553         }
4554       goto isone;               /* valid all zeros */
4555     }
4556  tnzro:
4557
4558   /* Test for infinity.  */
4559   if (y[NE - 1] == 0x7fff)
4560     {
4561       if (sign)
4562         sprintf (wstring, " -Infinity ");
4563       else
4564         sprintf (wstring, " Infinity ");
4565       goto bxit;
4566     }
4567
4568   /* Test for exponent nonzero but significand denormalized.
4569    * This is an error condition.
4570    */
4571   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4572     {
4573       mtherr ("etoasc", DOMAIN);
4574       sprintf (wstring, "NaN");
4575       goto bxit;
4576     }
4577
4578   /* Compare to 1.0 */
4579   i = ecmp (eone, y);
4580   if (i == 0)
4581     goto isone;
4582
4583   if (i == -2)
4584     abort ();
4585
4586   if (i < 0)
4587     {                           /* Number is greater than 1 */
4588       /* Convert significand to an integer and strip trailing decimal zeros.  */
4589       emov (y, u);
4590       u[NE - 1] = EXONE + NBITS - 1;
4591
4592       p = &etens[NTEN - 4][0];
4593       m = 16;
4594       do
4595         {
4596           ediv (p, u, t);
4597           efloor (t, w);
4598           for (j = 0; j < NE - 1; j++)
4599             {
4600               if (t[j] != w[j])
4601                 goto noint;
4602             }
4603           emov (t, u);
4604           expon += (int) m;
4605         noint:
4606           p += NE;
4607           m >>= 1;
4608         }
4609       while (m != 0);
4610
4611       /* Rescale from integer significand */
4612       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4613       emov (u, y);
4614       /* Find power of 10 */
4615       emov (eone, t);
4616       m = MAXP;
4617       p = &etens[0][0];
4618       /* An unordered compare result shouldn't happen here.  */
4619       while (ecmp (ten, u) <= 0)
4620         {
4621           if (ecmp (p, u) <= 0)
4622             {
4623               ediv (p, u, u);
4624               emul (p, t, t);
4625               expon += (int) m;
4626             }
4627           m >>= 1;
4628           if (m == 0)
4629             break;
4630           p += NE;
4631         }
4632     }
4633   else
4634     {                           /* Number is less than 1.0 */
4635       /* Pad significand with trailing decimal zeros.  */
4636       if (y[NE - 1] == 0)
4637         {
4638           while ((y[NE - 2] & 0x8000) == 0)
4639             {
4640               emul (ten, y, y);
4641               expon -= 1;
4642             }
4643         }
4644       else
4645         {
4646           emovi (y, w);
4647           for (i = 0; i < NDEC + 1; i++)
4648             {
4649               if ((w[NI - 1] & 0x7) != 0)
4650                 break;
4651               /* multiply by 10 */
4652               emovz (w, u);
4653               eshdn1 (u);
4654               eshdn1 (u);
4655               eaddm (w, u);
4656               u[1] += 3;
4657               while (u[2] != 0)
4658                 {
4659                   eshdn1 (u);
4660                   u[1] += 1;
4661                 }
4662               if (u[NI - 1] != 0)
4663                 break;
4664               if (eone[NE - 1] <= u[1])
4665                 break;
4666               emovz (u, w);
4667               expon -= 1;
4668             }
4669           emovo (w, y);
4670         }
4671       k = -MAXP;
4672       p = &emtens[0][0];
4673       r = &etens[0][0];
4674       emov (y, w);
4675       emov (eone, t);
4676       while (ecmp (eone, w) > 0)
4677         {
4678           if (ecmp (p, w) >= 0)
4679             {
4680               emul (r, w, w);
4681               emul (r, t, t);
4682               expon += k;
4683             }
4684           k /= 2;
4685           if (k == 0)
4686             break;
4687           p += NE;
4688           r += NE;
4689         }
4690       ediv (t, eone, t);
4691     }
4692  isone:
4693   /* Find the first (leading) digit.  */
4694   emovi (t, w);
4695   emovz (w, t);
4696   emovi (y, w);
4697   emovz (w, y);
4698   eiremain (t, y);
4699   digit = equot[NI - 1];
4700   while ((digit == 0) && (ecmp (y, ezero) != 0))
4701     {
4702       eshup1 (y);
4703       emovz (y, u);
4704       eshup1 (u);
4705       eshup1 (u);
4706       eaddm (u, y);
4707       eiremain (t, y);
4708       digit = equot[NI - 1];
4709       expon -= 1;
4710     }
4711   s = wstring;
4712   if (sign)
4713     *s++ = '-';
4714   else
4715     *s++ = ' ';
4716   /* Examine number of digits requested by caller.  */
4717   if (ndigs < 0)
4718     ndigs = 0;
4719   if (ndigs > NDEC)
4720     ndigs = NDEC;
4721   if (digit == 10)
4722     {
4723       *s++ = '1';
4724       *s++ = '.';
4725       if (ndigs > 0)
4726         {
4727           *s++ = '0';
4728           ndigs -= 1;
4729         }
4730       expon += 1;
4731     }
4732   else
4733     {
4734       *s++ = (char) digit + '0';
4735       *s++ = '.';
4736     }
4737   /* Generate digits after the decimal point.  */
4738   for (k = 0; k <= ndigs; k++)
4739     {
4740       /* multiply current number by 10, without normalizing */
4741       eshup1 (y);
4742       emovz (y, u);
4743       eshup1 (u);
4744       eshup1 (u);
4745       eaddm (u, y);
4746       eiremain (t, y);
4747       *s++ = (char) equot[NI - 1] + '0';
4748     }
4749   digit = equot[NI - 1];
4750   --s;
4751   ss = s;
4752   /* round off the ASCII string */
4753   if (digit > 4)
4754     {
4755       /* Test for critical rounding case in ASCII output.  */
4756       if (digit == 5)
4757         {
4758           emovo (y, t);
4759           if (ecmp (t, ezero) != 0)
4760             goto roun;          /* round to nearest */
4761 #ifndef C4X
4762           if ((*(s - 1) & 1) == 0)
4763             goto doexp;         /* round to even */
4764 #endif
4765         }
4766       /* Round up and propagate carry-outs */
4767     roun:
4768       --s;
4769       k = *s & CHARMASK;
4770       /* Carry out to most significant digit? */
4771       if (k == '.')
4772         {
4773           --s;
4774           k = *s;
4775           k += 1;
4776           *s = (char) k;
4777           /* Most significant digit carries to 10? */
4778           if (k > '9')
4779             {
4780               expon += 1;
4781               *s = '1';
4782             }
4783           goto doexp;
4784         }
4785       /* Round up and carry out from less significant digits */
4786       k += 1;
4787       *s = (char) k;
4788       if (k > '9')
4789         {
4790           *s = '0';
4791           goto roun;
4792         }
4793     }
4794  doexp:
4795   /* Strip trailing zeros, but leave at least one.  */
4796   while (ss[-1] == '0' && ss[-2] != '.')
4797     --ss;
4798   sprintf (ss, "e%d", expon);
4799  bxit:
4800   rndprc = rndsav;
4801   /* copy out the working string */
4802   s = string;
4803   ss = wstring;
4804   while (*ss == ' ')            /* strip possible leading space */
4805     ++ss;
4806   while ((*s++ = *ss++) != '\0')
4807     ;
4808 }
4809
4810
4811 /* Convert ASCII string to floating point.
4812
4813    Numeric input is a free format decimal number of any length, with
4814    or without decimal point.  Entering E after the number followed by an
4815    integer number causes the second number to be interpreted as a power of
4816    10 to be multiplied by the first number (i.e., "scientific" notation).  */
4817
4818 /* Convert ASCII string S to single precision float value Y.  */
4819
4820 static void
4821 asctoe24 (s, y)
4822      const char *s;
4823      UEMUSHORT *y;
4824 {
4825   asctoeg (s, y, 24);
4826 }
4827
4828
4829 /* Convert ASCII string S to double precision value Y.  */
4830
4831 static void
4832 asctoe53 (s, y)
4833      const char *s;
4834      UEMUSHORT *y;
4835 {
4836 #if defined(DEC) || defined(IBM)
4837   asctoeg (s, y, 56);
4838 #else
4839 #if defined(C4X)
4840   asctoeg (s, y, 32);
4841 #else
4842   asctoeg (s, y, 53);
4843 #endif
4844 #endif
4845 }
4846
4847
4848 /* Convert ASCII string S to double extended value Y.  */
4849
4850 static void
4851 asctoe64 (s, y)
4852      const char *s;
4853      UEMUSHORT *y;
4854 {
4855   asctoeg (s, y, 64);
4856 }
4857
4858 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
4859 /* Convert ASCII string S to 128-bit long double Y.  */
4860
4861 static void
4862 asctoe113 (s, y)
4863      const char *s;
4864      UEMUSHORT *y;
4865 {
4866   asctoeg (s, y, 113);
4867 }
4868 #endif
4869
4870 /* Convert ASCII string S to e type Y.  */
4871
4872 static void
4873 asctoe (s, y)
4874      const char *s;
4875      UEMUSHORT *y;
4876 {
4877   asctoeg (s, y, NBITS);
4878 }
4879
4880 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4881    of OPREC bits.  BASE is 16 for C99 hexadecimal floating constants.  */
4882
4883 static void
4884 asctoeg (ss, y, oprec)
4885      const char *ss;
4886      UEMUSHORT *y;
4887      int oprec;
4888 {
4889   UEMUSHORT yy[NI], xt[NI], tt[NI];
4890   int esign, decflg, sgnflg, nexp, exp, prec, lost;
4891   int i, k, trail, c, rndsav;
4892   EMULONG lexp;
4893   UEMUSHORT nsign;
4894   char *sp, *s, *lstr;
4895   int base = 10;
4896
4897   /* Copy the input string.  */
4898   lstr = (char *) alloca (strlen (ss) + 1);
4899
4900   while (*ss == ' ')            /* skip leading spaces */
4901     ++ss;
4902
4903   sp = lstr;
4904   while ((*sp++ = *ss++) != '\0')
4905     ;
4906   s = lstr;
4907
4908   if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
4909     {
4910       base = 16;
4911       s += 2;
4912     }
4913
4914   rndsav = rndprc;
4915   rndprc = NBITS;               /* Set to full precision */
4916   lost = 0;
4917   nsign = 0;
4918   decflg = 0;
4919   sgnflg = 0;
4920   nexp = 0;
4921   exp = 0;
4922   prec = 0;
4923   ecleaz (yy);
4924   trail = 0;
4925
4926  nxtcom:
4927   k = hex_value (*s);
4928   if ((k >= 0) && (k < base))
4929     {
4930       /* Ignore leading zeros */
4931       if ((prec == 0) && (decflg == 0) && (k == 0))
4932         goto donchr;
4933       /* Identify and strip trailing zeros after the decimal point.  */
4934       if ((trail == 0) && (decflg != 0))
4935         {
4936           sp = s;
4937           while (ISDIGIT (*sp) || (base == 16 && ISXDIGIT (*sp)))
4938             ++sp;
4939           /* Check for syntax error */
4940           c = *sp & CHARMASK;
4941           if ((base != 10 || ((c != 'e') && (c != 'E')))
4942               && (base != 16 || ((c != 'p') && (c != 'P')))
4943               && (c != '\0')
4944               && (c != '\n') && (c != '\r') && (c != ' ')
4945               && (c != ','))
4946             goto unexpected_char_error;
4947           --sp;
4948           while (*sp == '0')
4949             *sp-- = 'z';
4950           trail = 1;
4951           if (*s == 'z')
4952             goto donchr;
4953         }
4954
4955       /* If enough digits were given to more than fill up the yy register,
4956          continuing until overflow into the high guard word yy[2]
4957          guarantees that there will be a roundoff bit at the top
4958          of the low guard word after normalization.  */
4959
4960       if (yy[2] == 0)
4961         {
4962           if (base == 16)
4963             {
4964               if (decflg)
4965                 nexp += 4;      /* count digits after decimal point */
4966
4967               eshup1 (yy);      /* multiply current number by 16 */
4968               eshup1 (yy);
4969               eshup1 (yy);
4970               eshup1 (yy);
4971             }
4972           else
4973             {
4974               if (decflg)
4975                 nexp += 1;              /* count digits after decimal point */
4976
4977               eshup1 (yy);              /* multiply current number by 10 */
4978               emovz (yy, xt);
4979               eshup1 (xt);
4980               eshup1 (xt);
4981               eaddm (xt, yy);
4982             }
4983           /* Insert the current digit.  */
4984           ecleaz (xt);
4985           xt[NI - 2] = (UEMUSHORT) k;
4986           eaddm (xt, yy);
4987         }
4988       else
4989         {
4990           /* Mark any lost non-zero digit.  */
4991           lost |= k;
4992           /* Count lost digits before the decimal point.  */
4993           if (decflg == 0)
4994             {
4995               if (base == 10)
4996                 nexp -= 1;
4997               else
4998                 nexp -= 4;
4999             }
5000         }
5001       prec += 1;
5002       goto donchr;
5003     }
5004
5005   switch (*s)
5006     {
5007     case 'z':
5008       break;
5009     case 'E':
5010     case 'e':
5011     case 'P':
5012     case 'p':
5013       goto expnt;
5014     case '.':                   /* decimal point */
5015       if (decflg)
5016         goto unexpected_char_error;
5017       ++decflg;
5018       break;
5019     case '-':
5020       nsign = 0xffff;
5021       if (sgnflg)
5022         goto unexpected_char_error;
5023       ++sgnflg;
5024       break;
5025     case '+':
5026       if (sgnflg)
5027         goto unexpected_char_error;
5028       ++sgnflg;
5029       break;
5030     case ',':
5031     case ' ':
5032     case '\0':
5033     case '\n':
5034     case '\r':
5035       goto daldone;
5036     case 'i':
5037     case 'I':
5038       goto infinite;
5039     default:
5040     unexpected_char_error:
5041 #ifdef NANS
5042       einan (yy);
5043 #else
5044       mtherr ("asctoe", DOMAIN);
5045       eclear (yy);
5046 #endif
5047       goto aexit;
5048     }
5049  donchr:
5050   ++s;
5051   goto nxtcom;
5052
5053   /* Exponent interpretation */
5054  expnt:
5055   /* 0.0eXXX is zero, regardless of XXX.  Check for the 0.0.  */
5056   for (k = 0; k < NI; k++)
5057     {
5058       if (yy[k] != 0)
5059         goto read_expnt;
5060     }
5061   goto aexit;
5062
5063 read_expnt:
5064   esign = 1;
5065   exp = 0;
5066   ++s;
5067   /* check for + or - */
5068   if (*s == '-')
5069     {
5070       esign = -1;
5071       ++s;
5072     }
5073   if (*s == '+')
5074     ++s;
5075   while (ISDIGIT (*s))
5076     {
5077       exp *= 10;
5078       exp += *s++ - '0';
5079       if (exp > 999999)
5080         break;
5081     }
5082   if (esign < 0)
5083     exp = -exp;
5084   if ((exp > MAXDECEXP) && (base == 10))
5085     {
5086  infinite:
5087       ecleaz (yy);
5088       yy[E] = 0x7fff;           /* infinity */
5089       goto aexit;
5090     }
5091   if ((exp < MINDECEXP) && (base == 10))
5092     {
5093  zero:
5094       ecleaz (yy);
5095       goto aexit;
5096     }
5097
5098  daldone:
5099   if (base == 16)
5100     {
5101       /* Base 16 hexadecimal floating constant.  */
5102       if ((k = enormlz (yy)) > NBITS)
5103         {
5104           ecleaz (yy);
5105           goto aexit;
5106         }
5107       /* Adjust the exponent.  NEXP is the number of hex digits,
5108          EXP is a power of 2.  */
5109       lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5110       if (lexp > 0x7fff)
5111         goto infinite;
5112       if (lexp < 0)
5113         goto zero;
5114       yy[E] = lexp;
5115       goto expdon;
5116     }
5117
5118   nexp = exp - nexp;
5119   /* Pad trailing zeros to minimize power of 10, per IEEE spec.  */
5120   while ((nexp > 0) && (yy[2] == 0))
5121     {
5122       emovz (yy, xt);
5123       eshup1 (xt);
5124       eshup1 (xt);
5125       eaddm (yy, xt);
5126       eshup1 (xt);
5127       if (xt[2] != 0)
5128         break;
5129       nexp -= 1;
5130       emovz (xt, yy);
5131     }
5132   if ((k = enormlz (yy)) > NBITS)
5133     {
5134       ecleaz (yy);
5135       goto aexit;
5136     }
5137   lexp = (EXONE - 1 + NBITS) - k;
5138   emdnorm (yy, lost, 0, lexp, 64);
5139   lost = 0;
5140
5141   /* Convert to external format:
5142
5143      Multiply by 10**nexp.  If precision is 64 bits,
5144      the maximum relative error incurred in forming 10**n
5145      for 0 <= n <= 324 is 8.2e-20, at 10**180.
5146      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5147      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
5148
5149   lexp = yy[E];
5150   if (nexp == 0)
5151     {
5152       k = 0;
5153       goto expdon;
5154     }
5155   esign = 1;
5156   if (nexp < 0)
5157     {
5158       nexp = -nexp;
5159       esign = -1;
5160       if (nexp > 4096)
5161         {
5162           /* Punt.  Can't handle this without 2 divides.  */
5163           emovi (etens[0], tt);
5164           lexp -= tt[E];
5165           k = edivm (tt, yy);
5166           lexp += EXONE;
5167           nexp -= 4096;
5168         }
5169     }
5170   emov (eone, xt);
5171   exp = 1;
5172   i = NTEN;
5173   do
5174     {
5175       if (exp & nexp)
5176         emul (etens[i], xt, xt);
5177       i--;
5178       exp = exp + exp;
5179     }
5180   while (exp <= MAXP);
5181
5182   emovi (xt, tt);
5183   if (esign < 0)
5184     {
5185       lexp -= tt[E];
5186       k = edivm (tt, yy);
5187       lexp += EXONE;
5188     }
5189   else
5190     {
5191       lexp += tt[E];
5192       k = emulm (tt, yy);
5193       lexp -= EXONE - 1;
5194     }
5195   lost = k;
5196
5197  expdon:
5198
5199   /* Round and convert directly to the destination type */
5200   if (oprec == 53)
5201     lexp -= EXONE - 0x3ff;
5202 #ifdef C4X
5203   else if (oprec == 24 || oprec == 32)
5204     lexp -= (EXONE - 0x7f);
5205 #else
5206 #ifdef IBM
5207   else if (oprec == 24 || oprec == 56)
5208     lexp -= EXONE - (0x41 << 2);
5209 #else
5210 #ifdef DEC
5211   else if (oprec == 24)
5212     lexp -= dec_f.adjustment;
5213   else if (oprec == 56)
5214     {
5215       if (TARGET_G_FLOAT)
5216         lexp -= dec_g.adjustment;
5217       else
5218         lexp -= dec_d.adjustment;
5219     }
5220 #else
5221 else if (oprec == 24)
5222     lexp -= EXONE - 0177;
5223 #endif /* DEC */
5224 #endif /* IBM */
5225 #endif /* C4X */
5226   rndprc = oprec;
5227   emdnorm (yy, lost, 0, lexp, 64);
5228
5229  aexit:
5230
5231   rndprc = rndsav;
5232   yy[0] = nsign;
5233   switch (oprec)
5234     {
5235 #ifdef DEC
5236     case 56:
5237       todec (yy, y);
5238       break;
5239 #endif
5240 #ifdef IBM
5241     case 56:
5242       toibm (yy, y, DFmode);
5243       break;
5244 #endif
5245 #ifdef C4X
5246     case 32:
5247       toc4x (yy, y, HFmode);
5248       break;
5249 #endif
5250
5251     case 53:
5252       toe53 (yy, y);
5253       break;
5254     case 24:
5255       toe24 (yy, y);
5256       break;
5257     case 64:
5258       toe64 (yy, y);
5259       break;
5260 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5261     case 113:
5262       toe113 (yy, y);
5263       break;
5264 #endif
5265     case NBITS:
5266       emovo (yy, y);
5267       break;
5268     }
5269 }
5270
5271
5272
5273 /* Return Y = largest integer not greater than X (truncated toward minus
5274    infinity).  */
5275
5276 static const UEMUSHORT bmask[] =
5277 {
5278   0xffff,
5279   0xfffe,
5280   0xfffc,
5281   0xfff8,
5282   0xfff0,
5283   0xffe0,
5284   0xffc0,
5285   0xff80,
5286   0xff00,
5287   0xfe00,
5288   0xfc00,
5289   0xf800,
5290   0xf000,
5291   0xe000,
5292   0xc000,
5293   0x8000,
5294   0x0000,
5295 };
5296
5297 static void
5298 efloor (x, y)
5299      const UEMUSHORT x[];
5300      UEMUSHORT y[];
5301 {
5302   UEMUSHORT *p;
5303   int e, expon, i;
5304   UEMUSHORT f[NE];
5305
5306   emov (x, f);                  /* leave in external format */
5307   expon = (int) f[NE - 1];
5308   e = (expon & 0x7fff) - (EXONE - 1);
5309   if (e <= 0)
5310     {
5311       eclear (y);
5312       goto isitneg;
5313     }
5314   /* number of bits to clear out */
5315   e = NBITS - e;
5316   emov (f, y);
5317   if (e <= 0)
5318     return;
5319
5320   p = &y[0];
5321   while (e >= 16)
5322     {
5323       *p++ = 0;
5324       e -= 16;
5325     }
5326   /* clear the remaining bits */
5327   *p &= bmask[e];
5328   /* truncate negatives toward minus infinity */
5329  isitneg:
5330
5331   if ((UEMUSHORT) expon & (UEMUSHORT) 0x8000)
5332     {
5333       for (i = 0; i < NE - 1; i++)
5334         {
5335           if (f[i] != y[i])
5336             {
5337               esub (eone, y, y);
5338               break;
5339             }
5340         }
5341     }
5342 }
5343
5344
5345 #if 0
5346 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5347    For example, 1.1 = 0.55 * 2^1.  */
5348
5349 static void
5350 efrexp (x, exp, s)
5351      const UEMUSHORT x[];
5352      int *exp;
5353      UEMUSHORT s[];
5354 {
5355   UEMUSHORT xi[NI];
5356   EMULONG li;
5357
5358   emovi (x, xi);
5359   /*  Handle denormalized numbers properly using long integer exponent.  */
5360   li = (EMULONG) ((EMUSHORT) xi[1]);
5361
5362   if (li == 0)
5363     {
5364       li -= enormlz (xi);
5365     }
5366   xi[1] = 0x3ffe;
5367   emovo (xi, s);
5368   *exp = (int) (li - 0x3ffe);
5369 }
5370 #endif
5371
5372 /* Return e type Y = X * 2^PWR2.  */
5373
5374 static void
5375 eldexp (x, pwr2, y)
5376      const UEMUSHORT x[];
5377      int pwr2;
5378      UEMUSHORT y[];
5379 {
5380   UEMUSHORT xi[NI];
5381   EMULONG li;
5382   int i;
5383
5384   emovi (x, xi);
5385   li = xi[1];
5386   li += pwr2;
5387   i = 0;
5388   emdnorm (xi, i, i, li, !ROUND_TOWARDS_ZERO);
5389   emovo (xi, y);
5390 }
5391
5392
5393 #if 0
5394 /* C = remainder after dividing B by A, all e type values.
5395    Least significant integer quotient bits left in EQUOT.  */
5396
5397 static void
5398 eremain (a, b, c)
5399      const UEMUSHORT a[], b[];
5400      UEMUSHORT c[];
5401 {
5402   UEMUSHORT den[NI], num[NI];
5403
5404 #ifdef NANS
5405   if (eisinf (b)
5406       || (ecmp (a, ezero) == 0)
5407       || eisnan (a)
5408       || eisnan (b))
5409     {
5410       enan (c, 0);
5411       return;
5412     }
5413 #endif
5414   if (ecmp (a, ezero) == 0)
5415     {
5416       mtherr ("eremain", SING);
5417       eclear (c);
5418       return;
5419     }
5420   emovi (a, den);
5421   emovi (b, num);
5422   eiremain (den, num);
5423   /* Sign of remainder = sign of quotient */
5424   if (a[0] == b[0])
5425     num[0] = 0;
5426   else
5427     num[0] = 0xffff;
5428   emovo (num, c);
5429 }
5430 #endif
5431
5432 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5433     remainder in NUM.  */
5434
5435 static void
5436 eiremain (den, num)
5437      UEMUSHORT den[], num[];
5438 {
5439   EMULONG ld, ln;
5440   UEMUSHORT j;
5441
5442   ld = den[E];
5443   ld -= enormlz (den);
5444   ln = num[E];
5445   ln -= enormlz (num);
5446   ecleaz (equot);
5447   while (ln >= ld)
5448     {
5449       if (ecmpm (den, num) <= 0)
5450         {
5451           esubm (den, num);
5452           j = 1;
5453         }
5454       else
5455           j = 0;
5456       eshup1 (equot);
5457       equot[NI - 1] |= j;
5458       eshup1 (num);
5459       ln -= 1;
5460     }
5461   emdnorm (num, 0, 0, ln, 0);
5462 }
5463
5464 /* Report an error condition CODE encountered in function NAME.
5465
5466     Mnemonic        Value          Significance
5467
5468      DOMAIN            1       argument domain error
5469      SING              2       function singularity
5470      OVERFLOW          3       overflow range error
5471      UNDERFLOW         4       underflow range error
5472      TLOSS             5       total loss of precision
5473      PLOSS             6       partial loss of precision
5474      INVALID           7       NaN - producing operation
5475      EDOM             33       Unix domain error code
5476      ERANGE           34       Unix range error code
5477
5478    The order of appearance of the following messages is bound to the
5479    error codes defined above.  */
5480
5481 int merror = 0;
5482 extern int merror;
5483
5484 static void
5485 mtherr (name, code)
5486      const char *name;
5487      int code;
5488 {
5489   /* The string passed by the calling program is supposed to be the
5490      name of the function in which the error occurred.
5491      The code argument selects which error message string will be printed.  */
5492
5493   if (strcmp (name, "esub") == 0)
5494     name = "subtraction";
5495   else if (strcmp (name, "ediv") == 0)
5496     name = "division";
5497   else if (strcmp (name, "emul") == 0)
5498     name = "multiplication";
5499   else if (strcmp (name, "enormlz") == 0)
5500     name = "normalization";
5501   else if (strcmp (name, "etoasc") == 0)
5502     name = "conversion to text";
5503   else if (strcmp (name, "asctoe") == 0)
5504     name = "parsing";
5505   else if (strcmp (name, "eremain") == 0)
5506     name = "modulus";
5507   else if (strcmp (name, "esqrt") == 0)
5508     name = "square root";
5509   if (extra_warnings)
5510     {
5511       switch (code)
5512         {
5513         case DOMAIN:    warning ("%s: argument domain error"    , name); break;
5514         case SING:      warning ("%s: function singularity"     , name); break;
5515         case OVERFLOW:  warning ("%s: overflow range error"     , name); break;
5516         case UNDERFLOW: warning ("%s: underflow range error"    , name); break;
5517         case TLOSS:     warning ("%s: total loss of precision"  , name); break;
5518         case PLOSS:     warning ("%s: partial loss of precision", name); break;
5519         case INVALID:   warning ("%s: NaN - producing operation", name); break;
5520         default:        abort ();
5521         }
5522     }
5523
5524   /* Set global error message word */
5525   merror = code + 1;
5526 }
5527
5528 #ifdef DEC
5529 /* Convert DEC double precision D to e type E.  */
5530
5531 static void
5532 dectoe (d, e)
5533      const UEMUSHORT *d;
5534      UEMUSHORT *e;
5535 {
5536   if (TARGET_G_FLOAT)
5537     ieeetoe (d, e, &dec_g);
5538   else
5539     ieeetoe (d, e, &dec_d);
5540 }
5541
5542 /* Convert e type X to DEC double precision D.  */
5543
5544 static void
5545 etodec (x, d)
5546      const UEMUSHORT *x;
5547      UEMUSHORT *d;
5548 {
5549   UEMUSHORT xi[NI];
5550   EMULONG exp;
5551   int rndsav;
5552   const struct ieee_format *fmt;
5553
5554   if (TARGET_G_FLOAT)
5555     fmt = &dec_g;
5556   else
5557     fmt = &dec_d;
5558
5559   emovi (x, xi);
5560   /* Adjust exponent for offsets.  */
5561   exp = (EMULONG) xi[E] - fmt->adjustment;
5562   /* Round off to nearest or even.  */
5563   rndsav = rndprc;
5564   rndprc = fmt->precision;
5565   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5566   rndprc = rndsav;
5567   todec (xi, d);
5568 }
5569
5570 /* Convert exploded e-type X, that has already been rounded to
5571    56-bit precision, to DEC format double Y.  */
5572
5573 static void
5574 todec (x, y)
5575      UEMUSHORT *x, *y;
5576 {
5577   if (TARGET_G_FLOAT)
5578     toieee (x, y, &dec_g);
5579   else
5580     toieee (x, y, &dec_d);
5581 }
5582 #endif /* DEC */
5583
5584 #ifdef IBM
5585 /* Convert IBM single/double precision to e type.  */
5586
5587 static void
5588 ibmtoe (d, e, mode)
5589      const UEMUSHORT *d;
5590      UEMUSHORT *e;
5591      enum machine_mode mode;
5592 {
5593   UEMUSHORT y[NI];
5594   UEMUSHORT r, *p;
5595
5596   ecleaz (y);                   /* start with a zero */
5597   p = y;                        /* point to our number */
5598   r = *d;                       /* get IBM exponent word */
5599   if (*d & (unsigned int) 0x8000)
5600     *p = 0xffff;                /* fill in our sign */
5601   ++p;                          /* bump pointer to our exponent word */
5602   r &= 0x7f00;                  /* strip the sign bit */
5603   r >>= 6;                      /* shift exponent word down 6 bits */
5604                                 /* in fact shift by 8 right and 2 left */
5605   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5606                                 /* add our e type exponent offset */
5607   *p++ = r;                     /* to form our exponent */
5608
5609   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5610                                 /* strip off the IBM exponent and sign bits */
5611   if (mode != SFmode)           /* there are only 2 words in SFmode */
5612     {
5613       *p++ = *d++;              /* fill in the rest of our mantissa */
5614       *p++ = *d++;
5615     }
5616   *p = *d;
5617
5618   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5619     y[0] = y[E] = 0;
5620   else
5621     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5622                               /* handle change in RADIX */
5623   emovo (y, e);
5624 }
5625
5626
5627
5628 /* Convert e type to IBM single/double precision.  */
5629
5630 static void
5631 etoibm (x, d, mode)
5632      const UEMUSHORT *x;
5633      UEMUSHORT *d;
5634      enum machine_mode mode;
5635 {
5636   UEMUSHORT xi[NI];
5637   EMULONG exp;
5638   int rndsav;
5639
5640   emovi (x, xi);
5641   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5642                                                         /* round off to nearest or even */
5643   rndsav = rndprc;
5644   rndprc = 56;
5645   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5646   rndprc = rndsav;
5647   toibm (xi, d, mode);
5648 }
5649
5650 static void
5651 toibm (x, y, mode)
5652      UEMUSHORT *x, *y;
5653      enum machine_mode mode;
5654 {
5655   UEMUSHORT i;
5656   UEMUSHORT *p;
5657   int r;
5658
5659   p = x;
5660   *y = 0;
5661   if (*p++)
5662     *y = 0x8000;
5663   i = *p++;
5664   if (i == 0)
5665     {
5666       *y++ = 0;
5667       *y++ = 0;
5668       if (mode != SFmode)
5669         {
5670           *y++ = 0;
5671           *y++ = 0;
5672         }
5673       return;
5674     }
5675   r = i & 0x3;
5676   i >>= 2;
5677   if (i > 0x7f)
5678     {
5679       *y++ |= 0x7fff;
5680       *y++ = 0xffff;
5681       if (mode != SFmode)
5682         {
5683           *y++ = 0xffff;
5684           *y++ = 0xffff;
5685         }
5686 #ifdef ERANGE
5687       errno = ERANGE;
5688 #endif
5689       return;
5690     }
5691   i &= 0x7f;
5692   *y |= (i << 8);
5693   eshift (x, r + 5);
5694   *y++ |= x[M];
5695   *y++ = x[M + 1];
5696   if (mode != SFmode)
5697     {
5698       *y++ = x[M + 2];
5699       *y++ = x[M + 3];
5700     }
5701 }
5702 #endif /* IBM */
5703
5704
5705 #ifdef C4X
5706 /* Convert C4X single/double precision to e type.  */
5707
5708 static void
5709 c4xtoe (d, e, mode)
5710      const UEMUSHORT *d;
5711      UEMUSHORT *e;
5712      enum machine_mode mode;
5713 {
5714   UEMUSHORT y[NI];
5715   UEMUSHORT dn[4];
5716   int r;
5717   int isnegative;
5718   int size;
5719   int i;
5720   int carry;
5721
5722   dn[0] = d[0];
5723   dn[1] = d[1];
5724   if (mode != QFmode)
5725     {
5726       dn[2] = d[3] << 8;
5727       dn[3] = 0;
5728     }
5729
5730   /* Short-circuit the zero case.  */
5731   if ((dn[0] == 0x8000)
5732       && (dn[1] == 0x0000)
5733       && ((mode == QFmode) || ((dn[2] == 0x0000) && (dn[3] == 0x0000))))
5734     {
5735       e[0] = 0;
5736       e[1] = 0;
5737       e[2] = 0;
5738       e[3] = 0;
5739       e[4] = 0;
5740       e[5] = 0;
5741       return;
5742     }
5743
5744   ecleaz (y);                   /* start with a zero */
5745   r = dn[0];                    /* get sign/exponent part */
5746   if (r & (unsigned int) 0x0080)
5747     {
5748       y[0] = 0xffff;            /* fill in our sign */
5749       isnegative = TRUE;
5750     }
5751   else
5752     isnegative = FALSE;
5753
5754   r >>= 8;                      /* Shift exponent word down 8 bits.  */
5755   if (r & 0x80)                 /* Make the exponent negative if it is.  */
5756     r = r | (~0 & ~0xff);
5757
5758   if (isnegative)
5759     {
5760       /* Now do the high order mantissa.  We don't "or" on the high bit
5761          because it is 2 (not 1) and is handled a little differently
5762          below.  */
5763       y[M] = dn[0] & 0x7f;
5764
5765       y[M+1] = dn[1];
5766       if (mode != QFmode)       /* There are only 2 words in QFmode.  */
5767         {
5768           y[M+2] = dn[2];       /* Fill in the rest of our mantissa.  */
5769           y[M+3] = dn[3];
5770           size = 4;
5771         }
5772       else
5773         size = 2;
5774       eshift (y, -8);
5775
5776       /* Now do the two's complement on the data.  */
5777
5778       carry = 1;        /* Initially add 1 for the two's complement.  */
5779       for (i=size + M; i > M; i--)
5780         {
5781           if (carry && (y[i] == 0x0000))
5782             /* We overflowed into the next word, carry is the same.  */
5783             y[i] = carry ? 0x0000 : 0xffff;
5784           else
5785             {
5786               /* No overflow, just invert and add carry.  */
5787               y[i] = ((~y[i]) + carry) & 0xffff;
5788               carry = 0;
5789             }
5790         }
5791
5792       if (carry)
5793         {
5794           eshift (y, -1);
5795           y[M+1] |= 0x8000;
5796           r++;
5797         }
5798        y[1] = r + EXONE;
5799     }
5800   else
5801     {
5802       /* Add our e type exponent offset to form our exponent.  */
5803       r += EXONE;
5804       y[1] = r;
5805
5806      /* Now do the high order mantissa strip off the exponent and sign
5807         bits and add the high 1 bit.  */
5808      y[M] = (dn[0] & 0x7f) | 0x80;
5809
5810      y[M+1] = dn[1];
5811      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
5812        {
5813          y[M+2] = dn[2];        /* Fill in the rest of our mantissa.  */
5814          y[M+3] = dn[3];
5815        }
5816      eshift (y, -8);
5817     }
5818
5819   emovo (y, e);
5820 }
5821
5822
5823 /* Convert e type to C4X single/double precision.  */
5824
5825 static void
5826 etoc4x (x, d, mode)
5827      const UEMUSHORT *x;
5828      UEMUSHORT *d;
5829      enum machine_mode mode;
5830 {
5831   UEMUSHORT xi[NI];
5832   EMULONG exp;
5833   int rndsav;
5834
5835   emovi (x, xi);
5836
5837   /* Adjust exponent for offsets.  */
5838   exp = (EMULONG) xi[E] - (EXONE - 0x7f);
5839
5840   /* Round off to nearest or even.  */
5841   rndsav = rndprc;
5842   rndprc = mode == QFmode ? 24 : 32;
5843   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5844   rndprc = rndsav;
5845   toc4x (xi, d, mode);
5846 }
5847
5848 static void
5849 toc4x (x, y, mode)
5850      UEMUSHORT *x, *y;
5851      enum machine_mode mode;
5852 {
5853   int i;
5854   int v;
5855   int carry;
5856
5857   /* Short-circuit the zero case */
5858   if ((x[0] == 0)       /* Zero exponent and sign */
5859       && (x[1] == 0)
5860       && (x[M] == 0)    /* The rest is for zero mantissa */
5861       && (x[M+1] == 0)
5862       /* Only check for double if necessary */
5863       && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
5864     {
5865       /* We have a zero.  Put it into the output and return.  */
5866       *y++ = 0x8000;
5867       *y++ = 0x0000;
5868       if (mode != QFmode)
5869         {
5870           *y++ = 0x0000;
5871           *y++ = 0x0000;
5872         }
5873       return;
5874     }
5875
5876   *y = 0;
5877
5878   /* Negative number require a two's complement conversion of the
5879      mantissa.  */
5880   if (x[0])
5881     {
5882       *y = 0x0080;
5883
5884       i = ((int) x[1]) - 0x7f;
5885
5886       /* Now add 1 to the inverted data to do the two's complement.  */
5887       if (mode != QFmode)
5888         v = 4 + M;
5889       else
5890         v = 2 + M;
5891       carry = 1;
5892       while (v > M)
5893         {
5894           if (x[v] == 0x0000)
5895             x[v] = carry ? 0x0000 : 0xffff;
5896           else
5897             {
5898               x[v] = ((~x[v]) + carry) & 0xffff;
5899               carry = 0;
5900             }
5901           v--;
5902         }
5903
5904       /* The following is a special case.  The C4X negative float requires
5905          a zero in the high bit (because the format is (2 - x) x 2^m), so
5906          if a one is in that bit, we have to shift left one to get rid
5907          of it.  This only occurs if the number is -1 x 2^m.  */
5908       if (x[M+1] & 0x8000)
5909         {
5910           /* This is the case of -1 x 2^m, we have to rid ourselves of the
5911              high sign bit and shift the exponent.  */
5912           eshift (x, 1);
5913           i--;
5914         }
5915     }
5916   else
5917     i = ((int) x[1]) - 0x7f;
5918
5919   if ((i < -128) || (i > 127))
5920     {
5921       y[0] |= 0xff7f;
5922       y[1] = 0xffff;
5923       if (mode != QFmode)
5924         {
5925           y[2] = 0xffff;
5926           y[3] = 0xffff;
5927           y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
5928           y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
5929         }
5930 #ifdef ERANGE
5931       errno = ERANGE;
5932 #endif
5933       return;
5934     }
5935
5936   y[0] |= ((i & 0xff) << 8);
5937
5938   eshift (x, 8);
5939
5940   y[0] |= x[M] & 0x7f;
5941   y[1] = x[M + 1];
5942   if (mode != QFmode)
5943     {
5944       y[2] = x[M + 2];
5945       y[3] = x[M + 3];
5946       y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
5947       y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
5948     }
5949 }
5950 #endif /* C4X */
5951
5952 /* Output a binary NaN bit pattern in the target machine's format.  */
5953
5954 /* If special NaN bit patterns are required, define them in tm.h
5955    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
5956    patterns here.  */
5957 #ifdef TFMODE_NAN
5958 TFMODE_NAN;
5959 #else
5960 #ifdef IEEE
5961 static const UEMUSHORT TFbignan[8] =
5962  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5963 static const UEMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5964 #endif
5965 #endif
5966
5967 #ifdef XFMODE_NAN
5968 XFMODE_NAN;
5969 #else
5970 #ifdef IEEE
5971 static const UEMUSHORT XFbignan[6] =
5972  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5973 static const UEMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5974 #endif
5975 #endif
5976
5977 #ifdef DFMODE_NAN
5978 DFMODE_NAN;
5979 #else
5980 #ifdef IEEE
5981 static const UEMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5982 static const UEMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5983 #endif
5984 #endif
5985
5986 #ifdef SFMODE_NAN
5987 SFMODE_NAN;
5988 #else
5989 #ifdef IEEE
5990 static const UEMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5991 static const UEMUSHORT SFlittlenan[2] = {0, 0xffc0};
5992 #endif
5993 #endif
5994
5995
5996 #ifdef NANS
5997 static void
5998 make_nan (nan, sign, mode)
5999      UEMUSHORT *nan;
6000      int sign;
6001      enum machine_mode mode;
6002 {
6003   int n;
6004   const UEMUSHORT *p;
6005   int size;
6006
6007   size = GET_MODE_BITSIZE (mode);
6008   if (LARGEST_EXPONENT_IS_NORMAL (size))
6009     {
6010       warning ("%d-bit floats cannot hold NaNs", size);
6011       saturate (nan, sign, size, 0);
6012       return;
6013     }
6014   switch (mode)
6015     {
6016 /* Possibly the `reserved operand' patterns on a VAX can be
6017    used like NaN's, but probably not in the same way as IEEE.  */
6018 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6019     case TFmode:
6020 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6021       n = 8;
6022       if (REAL_WORDS_BIG_ENDIAN)
6023         p = TFbignan;
6024       else
6025         p = TFlittlenan;
6026       break;
6027 #endif
6028       /* FALLTHRU */
6029
6030     case XFmode:
6031       n = 6;
6032       if (REAL_WORDS_BIG_ENDIAN)
6033         p = XFbignan;
6034       else
6035         p = XFlittlenan;
6036       break;
6037
6038     case DFmode:
6039       n = 4;
6040       if (REAL_WORDS_BIG_ENDIAN)
6041         p = DFbignan;
6042       else
6043         p = DFlittlenan;
6044       break;
6045
6046     case SFmode:
6047     case HFmode:
6048       n = 2;
6049       if (REAL_WORDS_BIG_ENDIAN)
6050         p = SFbignan;
6051       else
6052         p = SFlittlenan;
6053       break;
6054 #endif
6055
6056     default:
6057       abort ();
6058     }
6059   if (REAL_WORDS_BIG_ENDIAN)
6060     *nan++ = (sign << 15) | (*p++ & 0x7fff);
6061   while (--n != 0)
6062     *nan++ = *p++;
6063   if (! REAL_WORDS_BIG_ENDIAN)
6064     *nan = (sign << 15) | (*p & 0x7fff);
6065 }
6066 #endif /* NANS */
6067
6068
6069 /* Create a saturation value for a SIZE-bit float, assuming that
6070    LARGEST_EXPONENT_IS_NORMAL (SIZE).
6071
6072    If SIGN is true, fill X with the most negative value, otherwise fill
6073    it with the most positive value.  WARN is true if the function should
6074    warn about overflow.  */
6075
6076 static void
6077 saturate (x, sign, size, warn)
6078      UEMUSHORT *x;
6079      int sign, size, warn;
6080 {
6081   int i;
6082
6083   if (warn && extra_warnings)
6084     warning ("value exceeds the range of a %d-bit float", size);
6085
6086   /* Create the most negative value.  */
6087   for (i = 0; i < size / EMUSHORT_SIZE; i++)
6088     x[i] = 0xffff;
6089
6090   /* Make it positive, if necessary.  */
6091   if (!sign)
6092     x[REAL_WORDS_BIG_ENDIAN? 0 : i - 1] = 0x7fff;
6093 }
6094
6095
6096 /* This is the inverse of the function `etarsingle' invoked by
6097    REAL_VALUE_TO_TARGET_SINGLE.  */
6098
6099 REAL_VALUE_TYPE
6100 ereal_unto_float (f)
6101      long f;
6102 {
6103   REAL_VALUE_TYPE r;
6104   UEMUSHORT s[2];
6105   UEMUSHORT e[NE];
6106
6107   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6108    This is the inverse operation to what the function `endian' does.  */
6109   if (REAL_WORDS_BIG_ENDIAN)
6110     {
6111       s[0] = (UEMUSHORT) (f >> 16);
6112       s[1] = (UEMUSHORT) f;
6113     }
6114   else
6115     {
6116       s[0] = (UEMUSHORT) f;
6117       s[1] = (UEMUSHORT) (f >> 16);
6118     }
6119   /* Convert and promote the target float to E-type.  */
6120   e24toe (s, e);
6121   /* Output E-type to REAL_VALUE_TYPE.  */
6122   PUT_REAL (e, &r);
6123   return r;
6124 }
6125
6126
6127 /* This is the inverse of the function `etardouble' invoked by
6128    REAL_VALUE_TO_TARGET_DOUBLE.  */
6129
6130 REAL_VALUE_TYPE
6131 ereal_unto_double (d)
6132      long d[];
6133 {
6134   REAL_VALUE_TYPE r;
6135   UEMUSHORT s[4];
6136   UEMUSHORT e[NE];
6137
6138   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6139   if (REAL_WORDS_BIG_ENDIAN)
6140     {
6141       s[0] = (UEMUSHORT) (d[0] >> 16);
6142       s[1] = (UEMUSHORT) d[0];
6143       s[2] = (UEMUSHORT) (d[1] >> 16);
6144       s[3] = (UEMUSHORT) d[1];
6145     }
6146   else
6147     {
6148       /* Target float words are little-endian.  */
6149       s[0] = (UEMUSHORT) d[0];
6150       s[1] = (UEMUSHORT) (d[0] >> 16);
6151       s[2] = (UEMUSHORT) d[1];
6152       s[3] = (UEMUSHORT) (d[1] >> 16);
6153     }
6154   /* Convert target double to E-type.  */
6155   e53toe (s, e);
6156   /* Output E-type to REAL_VALUE_TYPE.  */
6157   PUT_REAL (e, &r);
6158   return r;
6159 }
6160
6161
6162 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6163    This is somewhat like ereal_unto_float, but the input types
6164    for these are different.  */
6165
6166 REAL_VALUE_TYPE
6167 ereal_from_float (f)
6168      HOST_WIDE_INT f;
6169 {
6170   REAL_VALUE_TYPE r;
6171   UEMUSHORT s[2];
6172   UEMUSHORT e[NE];
6173
6174   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6175    This is the inverse operation to what the function `endian' does.  */
6176   if (REAL_WORDS_BIG_ENDIAN)
6177     {
6178       s[0] = (UEMUSHORT) (f >> 16);
6179       s[1] = (UEMUSHORT) f;
6180     }
6181   else
6182     {
6183       s[0] = (UEMUSHORT) f;
6184       s[1] = (UEMUSHORT) (f >> 16);
6185     }
6186   /* Convert and promote the target float to E-type.  */
6187   e24toe (s, e);
6188   /* Output E-type to REAL_VALUE_TYPE.  */
6189   PUT_REAL (e, &r);
6190   return r;
6191 }
6192
6193
6194 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6195    This is somewhat like ereal_unto_double, but the input types
6196    for these are different.
6197
6198    The DFmode is stored as an array of HOST_WIDE_INT in the target's
6199    data format, with no holes in the bit packing.  The first element
6200    of the input array holds the bits that would come first in the
6201    target computer's memory.  */
6202
6203 REAL_VALUE_TYPE
6204 ereal_from_double (d)
6205      HOST_WIDE_INT d[];
6206 {
6207   REAL_VALUE_TYPE r;
6208   UEMUSHORT s[4];
6209   UEMUSHORT e[NE];
6210
6211   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6212   if (REAL_WORDS_BIG_ENDIAN)
6213     {
6214 #if HOST_BITS_PER_WIDE_INT == 32
6215       s[0] = (UEMUSHORT) (d[0] >> 16);
6216       s[1] = (UEMUSHORT) d[0];
6217       s[2] = (UEMUSHORT) (d[1] >> 16);
6218       s[3] = (UEMUSHORT) d[1];
6219 #else
6220       /* In this case the entire target double is contained in the
6221          first array element.  The second element of the input is
6222          ignored.  */
6223       s[0] = (UEMUSHORT) (d[0] >> 48);
6224       s[1] = (UEMUSHORT) (d[0] >> 32);
6225       s[2] = (UEMUSHORT) (d[0] >> 16);
6226       s[3] = (UEMUSHORT) d[0];
6227 #endif
6228     }
6229   else
6230     {
6231       /* Target float words are little-endian.  */
6232       s[0] = (UEMUSHORT) d[0];
6233       s[1] = (UEMUSHORT) (d[0] >> 16);
6234 #if HOST_BITS_PER_WIDE_INT == 32
6235       s[2] = (UEMUSHORT) d[1];
6236       s[3] = (UEMUSHORT) (d[1] >> 16);
6237 #else
6238       s[2] = (UEMUSHORT) (d[0] >> 32);
6239       s[3] = (UEMUSHORT) (d[0] >> 48);
6240 #endif
6241     }
6242   /* Convert target double to E-type.  */
6243   e53toe (s, e);
6244   /* Output E-type to REAL_VALUE_TYPE.  */
6245   PUT_REAL (e, &r);
6246   return r;
6247 }
6248
6249
6250 #if 0
6251 /* Convert target computer unsigned 64-bit integer to e-type.
6252    The endian-ness of DImode follows the convention for integers,
6253    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
6254
6255 static void
6256 uditoe (di, e)
6257      const UEMUSHORT *di;  /* Address of the 64-bit int.  */
6258      UEMUSHORT *e;
6259 {
6260   UEMUSHORT yi[NI];
6261   int k;
6262
6263   ecleaz (yi);
6264   if (WORDS_BIG_ENDIAN)
6265     {
6266       for (k = M; k < M + 4; k++)
6267         yi[k] = *di++;
6268     }
6269   else
6270     {
6271       for (k = M + 3; k >= M; k--)
6272         yi[k] = *di++;
6273     }
6274   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6275   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6276     ecleaz (yi);                /* it was zero */
6277   else
6278     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6279   emovo (yi, e);
6280 }
6281
6282 /* Convert target computer signed 64-bit integer to e-type.  */
6283
6284 static void
6285 ditoe (di, e)
6286      const UEMUSHORT *di;  /* Address of the 64-bit int.  */
6287      UEMUSHORT *e;
6288 {
6289   unsigned EMULONG acc;
6290   UEMUSHORT yi[NI];
6291   UEMUSHORT carry;
6292   int k, sign;
6293
6294   ecleaz (yi);
6295   if (WORDS_BIG_ENDIAN)
6296     {
6297       for (k = M; k < M + 4; k++)
6298         yi[k] = *di++;
6299     }
6300   else
6301     {
6302       for (k = M + 3; k >= M; k--)
6303         yi[k] = *di++;
6304     }
6305   /* Take absolute value */
6306   sign = 0;
6307   if (yi[M] & 0x8000)
6308     {
6309       sign = 1;
6310       carry = 0;
6311       for (k = M + 3; k >= M; k--)
6312         {
6313           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6314           yi[k] = acc;
6315           carry = 0;
6316           if (acc & 0x10000)
6317             carry = 1;
6318         }
6319     }
6320   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6321   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6322     ecleaz (yi);                /* it was zero */
6323   else
6324     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6325   emovo (yi, e);
6326   if (sign)
6327         eneg (e);
6328 }
6329
6330
6331 /* Convert e-type to unsigned 64-bit int.  */
6332
6333 static void
6334 etoudi (x, i)
6335      const UEMUSHORT *x;
6336      UEMUSHORT *i;
6337 {
6338   UEMUSHORT xi[NI];
6339   int j, k;
6340
6341   emovi (x, xi);
6342   if (xi[0])
6343     {
6344       xi[M] = 0;
6345       goto noshift;
6346     }
6347   k = (int) xi[E] - (EXONE - 1);
6348   if (k <= 0)
6349     {
6350       for (j = 0; j < 4; j++)
6351         *i++ = 0;
6352       return;
6353     }
6354   if (k > 64)
6355     {
6356       for (j = 0; j < 4; j++)
6357         *i++ = 0xffff;
6358       if (extra_warnings)
6359         warning ("overflow on truncation to integer");
6360       return;
6361     }
6362   if (k > 16)
6363     {
6364       /* Shift more than 16 bits: first shift up k-16 mod 16,
6365          then shift up by 16's.  */
6366       j = k - ((k >> 4) << 4);
6367       if (j == 0)
6368         j = 16;
6369       eshift (xi, j);
6370       if (WORDS_BIG_ENDIAN)
6371         *i++ = xi[M];
6372       else
6373         {
6374           i += 3;
6375           *i-- = xi[M];
6376         }
6377       k -= j;
6378       do
6379         {
6380           eshup6 (xi);
6381           if (WORDS_BIG_ENDIAN)
6382             *i++ = xi[M];
6383           else
6384             *i-- = xi[M];
6385         }
6386       while ((k -= 16) > 0);
6387     }
6388   else
6389     {
6390       /* shift not more than 16 bits */
6391       eshift (xi, k);
6392
6393 noshift:
6394
6395       if (WORDS_BIG_ENDIAN)
6396         {
6397           i += 3;
6398           *i-- = xi[M];
6399           *i-- = 0;
6400           *i-- = 0;
6401           *i = 0;
6402         }
6403       else
6404         {
6405           *i++ = xi[M];
6406           *i++ = 0;
6407           *i++ = 0;
6408           *i = 0;
6409         }
6410     }
6411 }
6412
6413
6414 /* Convert e-type to signed 64-bit int.  */
6415
6416 static void
6417 etodi (x, i)
6418      const UEMUSHORT *x;
6419      UEMUSHORT *i;
6420 {
6421   unsigned EMULONG acc;
6422   UEMUSHORT xi[NI];
6423   UEMUSHORT carry;
6424   UEMUSHORT *isave;
6425   int j, k;
6426
6427   emovi (x, xi);
6428   k = (int) xi[E] - (EXONE - 1);
6429   if (k <= 0)
6430     {
6431       for (j = 0; j < 4; j++)
6432         *i++ = 0;
6433       return;
6434     }
6435   if (k > 64)
6436     {
6437       for (j = 0; j < 4; j++)
6438         *i++ = 0xffff;
6439       if (extra_warnings)
6440         warning ("overflow on truncation to integer");
6441       return;
6442     }
6443   isave = i;
6444   if (k > 16)
6445     {
6446       /* Shift more than 16 bits: first shift up k-16 mod 16,
6447          then shift up by 16's.  */
6448       j = k - ((k >> 4) << 4);
6449       if (j == 0)
6450         j = 16;
6451       eshift (xi, j);
6452       if (WORDS_BIG_ENDIAN)
6453         *i++ = xi[M];
6454       else
6455         {
6456           i += 3;
6457           *i-- = xi[M];
6458         }
6459       k -= j;
6460       do
6461         {
6462           eshup6 (xi);
6463           if (WORDS_BIG_ENDIAN)
6464             *i++ = xi[M];
6465           else
6466             *i-- = xi[M];
6467         }
6468       while ((k -= 16) > 0);
6469     }
6470   else
6471     {
6472       /* shift not more than 16 bits */
6473       eshift (xi, k);
6474
6475       if (WORDS_BIG_ENDIAN)
6476         {
6477           i += 3;
6478           *i = xi[M];
6479           *i-- = 0;
6480           *i-- = 0;
6481           *i = 0;
6482         }
6483       else
6484         {
6485           *i++ = xi[M];
6486           *i++ = 0;
6487           *i++ = 0;
6488           *i = 0;
6489         }
6490     }
6491   /* Negate if negative */
6492   if (xi[0])
6493     {
6494       carry = 0;
6495       if (WORDS_BIG_ENDIAN)
6496         isave += 3;
6497       for (k = 0; k < 4; k++)
6498         {
6499           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6500           if (WORDS_BIG_ENDIAN)
6501             *isave-- = acc;
6502           else
6503             *isave++ = acc;
6504           carry = 0;
6505           if (acc & 0x10000)
6506             carry = 1;
6507         }
6508     }
6509 }
6510
6511
6512 /* Longhand square root routine.  */
6513
6514
6515 static int esqinited = 0;
6516 static unsigned short sqrndbit[NI];
6517
6518 static void
6519 esqrt (x, y)
6520      const UEMUSHORT *x;
6521      UEMUSHORT *y;
6522 {
6523   UEMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6524   EMULONG m, exp;
6525   int i, j, k, n, nlups;
6526
6527   if (esqinited == 0)
6528     {
6529       ecleaz (sqrndbit);
6530       sqrndbit[NI - 2] = 1;
6531       esqinited = 1;
6532     }
6533   /* Check for arg <= 0 */
6534   i = ecmp (x, ezero);
6535   if (i <= 0)
6536     {
6537       if (i == -1)
6538         {
6539           mtherr ("esqrt", DOMAIN);
6540           eclear (y);
6541         }
6542       else
6543         emov (x, y);
6544       return;
6545     }
6546
6547 #ifdef INFINITY
6548   if (eisinf (x))
6549     {
6550       eclear (y);
6551       einfin (y);
6552       return;
6553     }
6554 #endif
6555   /* Bring in the arg and renormalize if it is denormal.  */
6556   emovi (x, xx);
6557   m = (EMULONG) xx[1];          /* local long word exponent */
6558   if (m == 0)
6559     m -= enormlz (xx);
6560
6561   /* Divide exponent by 2 */
6562   m -= 0x3ffe;
6563   exp = (unsigned short) ((m / 2) + 0x3ffe);
6564
6565   /* Adjust if exponent odd */
6566   if ((m & 1) != 0)
6567     {
6568       if (m > 0)
6569         exp += 1;
6570       eshdn1 (xx);
6571     }
6572
6573   ecleaz (sq);
6574   ecleaz (num);
6575   n = 8;                        /* get 8 bits of result per inner loop */
6576   nlups = rndprc;
6577   j = 0;
6578
6579   while (nlups > 0)
6580     {
6581       /* bring in next word of arg */
6582       if (j < NE)
6583         num[NI - 1] = xx[j + 3];
6584       /* Do additional bit on last outer loop, for roundoff.  */
6585       if (nlups <= 8)
6586         n = nlups + 1;
6587       for (i = 0; i < n; i++)
6588         {
6589           /* Next 2 bits of arg */
6590           eshup1 (num);
6591           eshup1 (num);
6592           /* Shift up answer */
6593           eshup1 (sq);
6594           /* Make trial divisor */
6595           for (k = 0; k < NI; k++)
6596             temp[k] = sq[k];
6597           eshup1 (temp);
6598           eaddm (sqrndbit, temp);
6599           /* Subtract and insert answer bit if it goes in */
6600           if (ecmpm (temp, num) <= 0)
6601             {
6602               esubm (temp, num);
6603               sq[NI - 2] |= 1;
6604             }
6605         }
6606       nlups -= n;
6607       j += 1;
6608     }
6609
6610   /* Adjust for extra, roundoff loop done.  */
6611   exp += (NBITS - 1) - rndprc;
6612
6613   /* Sticky bit = 1 if the remainder is nonzero.  */
6614   k = 0;
6615   for (i = 3; i < NI; i++)
6616     k |= (int) num[i];
6617
6618   /* Renormalize and round off.  */
6619   emdnorm (sq, k, 0, exp, !ROUND_TOWARDS_ZERO);
6620   emovo (sq, y);
6621 }
6622 #endif
6623 \f
6624 /* Return the binary precision of the significand for a given
6625    floating point mode.  The mode can hold an integer value
6626    that many bits wide, without losing any bits.  */
6627
6628 unsigned int
6629 significand_size (mode)
6630      enum machine_mode mode;
6631 {
6632
6633   /* Don't test the modes, but their sizes, lest this
6634      code won't work for BITS_PER_UNIT != 8 .  */
6635
6636   switch (GET_MODE_BITSIZE (mode))
6637     {
6638      case 32:
6639
6640 #ifdef C4X
6641        return 56;
6642 #else
6643        return 24;
6644 #endif
6645
6646      case 64:
6647 #ifdef IEEE
6648        return 53;
6649 #else
6650        return 56;
6651 #endif
6652
6653      case 96:
6654        return 64;
6655
6656      case 128:
6657 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6658        return 113;
6659 #else
6660        return 64;
6661 #endif
6662
6663      default:
6664        abort ();
6665     }
6666 }