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