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