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