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   eclear (e);
1902   einfin (e);
1903   rndprc = rndsav;
1904
1905   PUT_REAL (e, &r);
1906   return r;
1907 }
1908
1909 /* Output an e-type NaN.
1910    This generates Intel's quiet NaN pattern for extended real.
1911    The exponent is 7fff, the leading mantissa word is c000.  */
1912
1913 #ifdef NANS
1914 static void
1915 enan (x, sign)
1916      UEMUSHORT *x;
1917      int sign;
1918 {
1919   int i;
1920
1921   for (i = 0; i < NE - 2; i++)
1922     *x++ = 0;
1923   *x++ = 0xc000;
1924   *x = (sign << 15) | 0x7fff;
1925 }
1926 #endif /* NANS */
1927
1928 /* Move in an e-type number A, converting it to exploded e-type B.  */
1929
1930 static void
1931 emovi (a, b)
1932      const UEMUSHORT *a;
1933      UEMUSHORT *b;
1934 {
1935   const UEMUSHORT *p;
1936   UEMUSHORT *q;
1937   int i;
1938
1939   q = b;
1940   p = a + (NE - 1);             /* point to last word of external number */
1941   /* get the sign bit */
1942   if (*p & 0x8000)
1943     *q++ = 0xffff;
1944   else
1945     *q++ = 0;
1946   /* get the exponent */
1947   *q = *p--;
1948   *q++ &= 0x7fff;               /* delete the sign bit */
1949 #ifdef INFINITY
1950   if ((*(q - 1) & 0x7fff) == 0x7fff)
1951     {
1952 #ifdef NANS
1953       if (eisnan (a))
1954         {
1955           *q++ = 0;
1956           for (i = 3; i < NI; i++)
1957             *q++ = *p--;
1958           return;
1959         }
1960 #endif
1961
1962       for (i = 2; i < NI; i++)
1963         *q++ = 0;
1964       return;
1965     }
1966 #endif
1967
1968   /* clear high guard word */
1969   *q++ = 0;
1970   /* move in the significand */
1971   for (i = 0; i < NE - 1; i++)
1972     *q++ = *p--;
1973   /* clear low guard word */
1974   *q = 0;
1975 }
1976
1977 /* Move out exploded e-type number A, converting it to e type B.  */
1978
1979 static void
1980 emovo (a, b)
1981      const UEMUSHORT *a;
1982      UEMUSHORT *b;
1983 {
1984   const UEMUSHORT *p;
1985   UEMUSHORT *q;
1986   UEMUSHORT i;
1987   int j;
1988
1989   p = a;
1990   q = b + (NE - 1);             /* point to output exponent */
1991   /* combine sign and exponent */
1992   i = *p++;
1993   if (i)
1994     *q-- = *p++ | 0x8000;
1995   else
1996     *q-- = *p++;
1997 #ifdef INFINITY
1998   if (*(p - 1) == 0x7fff)
1999     {
2000 #ifdef NANS
2001       if (eiisnan (a))
2002         {
2003           enan (b, eiisneg (a));
2004           return;
2005         }
2006 #endif
2007       einfin (b);
2008         return;
2009     }
2010 #endif
2011   /* skip over guard word */
2012   ++p;
2013   /* move the significand */
2014   for (j = 0; j < NE - 1; j++)
2015     *q-- = *p++;
2016 }
2017
2018 /* Clear out exploded e-type number XI.  */
2019
2020 static void
2021 ecleaz (xi)
2022      UEMUSHORT *xi;
2023 {
2024   int i;
2025
2026   for (i = 0; i < NI; i++)
2027     *xi++ = 0;
2028 }
2029
2030 /* Clear out exploded e-type XI, but don't touch the sign.  */
2031
2032 static void
2033 ecleazs (xi)
2034      UEMUSHORT *xi;
2035 {
2036   int i;
2037
2038   ++xi;
2039   for (i = 0; i < NI - 1; i++)
2040     *xi++ = 0;
2041 }
2042
2043 /* Move exploded e-type number from A to B.  */
2044
2045 static void
2046 emovz (a, b)
2047      const UEMUSHORT *a;
2048      UEMUSHORT *b;
2049 {
2050   int i;
2051
2052   for (i = 0; i < NI - 1; i++)
2053     *b++ = *a++;
2054   /* clear low guard word */
2055   *b = 0;
2056 }
2057
2058 /* Generate exploded e-type NaN.
2059    The explicit pattern for this is maximum exponent and
2060    top two significant bits set.  */
2061
2062 #ifdef NANS
2063 static void
2064 einan (x)
2065      UEMUSHORT x[];
2066 {
2067
2068   ecleaz (x);
2069   x[E] = 0x7fff;
2070   x[M + 1] = 0xc000;
2071 }
2072 #endif /* NANS */
2073
2074 /* Return nonzero if exploded e-type X is a NaN.  */
2075
2076 #ifdef NANS
2077 static int
2078 eiisnan (x)
2079      const UEMUSHORT x[];
2080 {
2081   int i;
2082
2083   if ((x[E] & 0x7fff) == 0x7fff)
2084     {
2085       for (i = M + 1; i < NI; i++)
2086         {
2087           if (x[i] != 0)
2088             return (1);
2089         }
2090     }
2091   return (0);
2092 }
2093 #endif /* NANS */
2094
2095 /* Return nonzero if sign of exploded e-type X is nonzero.  */
2096
2097 static int
2098 eiisneg (x)
2099      const UEMUSHORT x[];
2100 {
2101
2102   return x[0] != 0;
2103 }
2104
2105 #if 0
2106 /* Fill exploded e-type X with infinity pattern.
2107    This has maximum exponent and significand all zeros.  */
2108
2109 static void
2110 eiinfin (x)
2111      UEMUSHORT x[];
2112 {
2113
2114   ecleaz (x);
2115   x[E] = 0x7fff;
2116 }
2117 #endif /* 0 */
2118
2119 /* Return nonzero if exploded e-type X is infinite.  */
2120
2121 #ifdef INFINITY
2122 static int
2123 eiisinf (x)
2124      const UEMUSHORT x[];
2125 {
2126
2127 #ifdef NANS
2128   if (eiisnan (x))
2129     return (0);
2130 #endif
2131   if ((x[E] & 0x7fff) == 0x7fff)
2132     return (1);
2133   return (0);
2134 }
2135 #endif /* INFINITY */
2136
2137 /* Compare significands of numbers in internal exploded e-type format.
2138    Guard words are included in the comparison.
2139
2140    Returns      +1 if a > b
2141                  0 if a == b
2142                 -1 if a < b   */
2143
2144 static int
2145 ecmpm (a, b)
2146      const UEMUSHORT *a, *b;
2147 {
2148   int i;
2149
2150   a += M;                       /* skip up to significand area */
2151   b += M;
2152   for (i = M; i < NI; i++)
2153     {
2154       if (*a++ != *b++)
2155         goto difrnt;
2156     }
2157   return (0);
2158
2159  difrnt:
2160   if (*(--a) > *(--b))
2161     return (1);
2162   else
2163     return (-1);
2164 }
2165
2166 /* Shift significand of exploded e-type X down by 1 bit.  */
2167
2168 static void
2169 eshdn1 (x)
2170      UEMUSHORT *x;
2171 {
2172   UEMUSHORT bits;
2173   int i;
2174
2175   x += M;                       /* point to significand area */
2176
2177   bits = 0;
2178   for (i = M; i < NI; i++)
2179     {
2180       if (*x & 1)
2181         bits |= 1;
2182       *x >>= 1;
2183       if (bits & 2)
2184         *x |= 0x8000;
2185       bits <<= 1;
2186       ++x;
2187     }
2188 }
2189
2190 /* Shift significand of exploded e-type X up by 1 bit.  */
2191
2192 static void
2193 eshup1 (x)
2194      UEMUSHORT *x;
2195 {
2196   UEMUSHORT bits;
2197   int i;
2198
2199   x += NI - 1;
2200   bits = 0;
2201
2202   for (i = M; i < NI; i++)
2203     {
2204       if (*x & 0x8000)
2205         bits |= 1;
2206       *x <<= 1;
2207       if (bits & 2)
2208         *x |= 1;
2209       bits <<= 1;
2210       --x;
2211     }
2212 }
2213
2214
2215 /* Shift significand of exploded e-type X down by 8 bits.  */
2216
2217 static void
2218 eshdn8 (x)
2219      UEMUSHORT *x;
2220 {
2221   UEMUSHORT newbyt, oldbyt;
2222   int i;
2223
2224   x += M;
2225   oldbyt = 0;
2226   for (i = M; i < NI; i++)
2227     {
2228       newbyt = *x << 8;
2229       *x >>= 8;
2230       *x |= oldbyt;
2231       oldbyt = newbyt;
2232       ++x;
2233     }
2234 }
2235
2236 /* Shift significand of exploded e-type X up by 8 bits.  */
2237
2238 static void
2239 eshup8 (x)
2240      UEMUSHORT *x;
2241 {
2242   int i;
2243   UEMUSHORT newbyt, oldbyt;
2244
2245   x += NI - 1;
2246   oldbyt = 0;
2247
2248   for (i = M; i < NI; i++)
2249     {
2250       newbyt = *x >> 8;
2251       *x <<= 8;
2252       *x |= oldbyt;
2253       oldbyt = newbyt;
2254       --x;
2255     }
2256 }
2257
2258 /* Shift significand of exploded e-type X up by 16 bits.  */
2259
2260 static void
2261 eshup6 (x)
2262      UEMUSHORT *x;
2263 {
2264   int i;
2265   UEMUSHORT *p;
2266
2267   p = x + M;
2268   x += M + 1;
2269
2270   for (i = M; i < NI - 1; i++)
2271     *p++ = *x++;
2272
2273   *p = 0;
2274 }
2275
2276 /* Shift significand of exploded e-type X down by 16 bits.  */
2277
2278 static void
2279 eshdn6 (x)
2280      UEMUSHORT *x;
2281 {
2282   int i;
2283   UEMUSHORT *p;
2284
2285   x += NI - 1;
2286   p = x + 1;
2287
2288   for (i = M; i < NI - 1; i++)
2289     *(--p) = *(--x);
2290
2291   *(--p) = 0;
2292 }
2293
2294 /* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */
2295
2296 static void
2297 eaddm (x, y)
2298      const UEMUSHORT *x;
2299      UEMUSHORT *y;
2300 {
2301   unsigned EMULONG a;
2302   int i;
2303   unsigned int carry;
2304
2305   x += NI - 1;
2306   y += NI - 1;
2307   carry = 0;
2308   for (i = M; i < NI; i++)
2309     {
2310       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
2311       if (a & 0x10000)
2312         carry = 1;
2313       else
2314         carry = 0;
2315       *y = (UEMUSHORT) a;
2316       --x;
2317       --y;
2318     }
2319 }
2320
2321 /* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */
2322
2323 static void
2324 esubm (x, y)
2325      const UEMUSHORT *x;
2326      UEMUSHORT *y;
2327 {
2328   unsigned EMULONG a;
2329   int i;
2330   unsigned int carry;
2331
2332   x += NI - 1;
2333   y += NI - 1;
2334   carry = 0;
2335   for (i = M; i < NI; i++)
2336     {
2337       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
2338       if (a & 0x10000)
2339         carry = 1;
2340       else
2341         carry = 0;
2342       *y = (UEMUSHORT) a;
2343       --x;
2344       --y;
2345     }
2346 }
2347
2348
2349 static UEMUSHORT equot[NI];
2350
2351
2352 #if 0
2353 /* Radix 2 shift-and-add versions of multiply and divide  */
2354
2355
2356 /* Divide significands */
2357
2358 int
2359 edivm (den, num)
2360      UEMUSHORT den[], num[];
2361 {
2362   int i;
2363   UEMUSHORT *p, *q;
2364   UEMUSHORT j;
2365
2366   p = &equot[0];
2367   *p++ = num[0];
2368   *p++ = num[1];
2369
2370   for (i = M; i < NI; i++)
2371     {
2372       *p++ = 0;
2373     }
2374
2375   /* Use faster compare and subtraction if denominator has only 15 bits of
2376      significance.  */
2377
2378   p = &den[M + 2];
2379   if (*p++ == 0)
2380     {
2381       for (i = M + 3; i < NI; i++)
2382         {
2383           if (*p++ != 0)
2384             goto fulldiv;
2385         }
2386       if ((den[M + 1] & 1) != 0)
2387         goto fulldiv;
2388       eshdn1 (num);
2389       eshdn1 (den);
2390
2391       p = &den[M + 1];
2392       q = &num[M + 1];
2393
2394       for (i = 0; i < NBITS + 2; i++)
2395         {
2396           if (*p <= *q)
2397             {
2398               *q -= *p;
2399               j = 1;
2400             }
2401           else
2402             {
2403               j = 0;
2404             }
2405           eshup1 (equot);
2406           equot[NI - 2] |= j;
2407           eshup1 (num);
2408         }
2409       goto divdon;
2410     }
2411
2412   /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2413      bit + 1 roundoff bit.  */
2414
2415  fulldiv:
2416
2417   p = &equot[NI - 2];
2418   for (i = 0; i < NBITS + 2; i++)
2419     {
2420       if (ecmpm (den, num) <= 0)
2421         {
2422           esubm (den, num);
2423           j = 1;                /* quotient bit = 1 */
2424         }
2425       else
2426         j = 0;
2427       eshup1 (equot);
2428       *p |= j;
2429       eshup1 (num);
2430     }
2431
2432  divdon:
2433
2434   eshdn1 (equot);
2435   eshdn1 (equot);
2436
2437   /* test for nonzero remainder after roundoff bit */
2438   p = &num[M];
2439   j = 0;
2440   for (i = M; i < NI; i++)
2441     {
2442       j |= *p++;
2443     }
2444   if (j)
2445     j = 1;
2446
2447
2448   for (i = 0; i < NI; i++)
2449     num[i] = equot[i];
2450   return ((int) j);
2451 }
2452
2453
2454 /* Multiply significands */
2455
2456 int
2457 emulm (a, b)
2458      UEMUSHORT a[], b[];
2459 {
2460   UEMUSHORT *p, *q;
2461   int i, j, k;
2462
2463   equot[0] = b[0];
2464   equot[1] = b[1];
2465   for (i = M; i < NI; i++)
2466     equot[i] = 0;
2467
2468   p = &a[NI - 2];
2469   k = NBITS;
2470   while (*p == 0)               /* significand is not supposed to be zero */
2471     {
2472       eshdn6 (a);
2473       k -= 16;
2474     }
2475   if ((*p & 0xff) == 0)
2476     {
2477       eshdn8 (a);
2478       k -= 8;
2479     }
2480
2481   q = &equot[NI - 1];
2482   j = 0;
2483   for (i = 0; i < k; i++)
2484     {
2485       if (*p & 1)
2486         eaddm (b, equot);
2487       /* remember if there were any nonzero bits shifted out */
2488       if (*q & 1)
2489         j |= 1;
2490       eshdn1 (a);
2491       eshdn1 (equot);
2492     }
2493
2494   for (i = 0; i < NI; i++)
2495     b[i] = equot[i];
2496
2497   /* return flag for lost nonzero bits */
2498   return (j);
2499 }
2500
2501 #else
2502
2503 /* Radix 65536 versions of multiply and divide.  */
2504
2505 /* Multiply significand of e-type number B
2506    by 16-bit quantity A, return e-type result to C.  */
2507
2508 static void
2509 m16m (a, b, c)
2510      unsigned int a;
2511      const UEMUSHORT b[];
2512      UEMUSHORT c[];
2513 {
2514   UEMUSHORT *pp;
2515   unsigned EMULONG carry;
2516   const UEMUSHORT *ps;
2517   UEMUSHORT p[NI];
2518   unsigned EMULONG aa, m;
2519   int i;
2520
2521   aa = a;
2522   pp = &p[NI-2];
2523   *pp++ = 0;
2524   *pp = 0;
2525   ps = &b[NI-1];
2526
2527   for (i=M+1; i<NI; i++)
2528     {
2529       if (*ps == 0)
2530         {
2531           --ps;
2532           --pp;
2533           *(pp-1) = 0;
2534         }
2535       else
2536         {
2537           m = (unsigned EMULONG) aa * *ps--;
2538           carry = (m & 0xffff) + *pp;
2539           *pp-- = (UEMUSHORT) carry;
2540           carry = (carry >> 16) + (m >> 16) + *pp;
2541           *pp = (UEMUSHORT) carry;
2542           *(pp-1) = carry >> 16;
2543         }
2544     }
2545   for (i=M; i<NI; i++)
2546     c[i] = p[i];
2547 }
2548
2549 /* Divide significands of exploded e-types NUM / DEN.  Neither the
2550    numerator NUM nor the denominator DEN is permitted to have its high guard
2551    word nonzero.  */
2552
2553 static int
2554 edivm (den, num)
2555      const UEMUSHORT den[];
2556      UEMUSHORT num[];
2557 {
2558   int i;
2559   UEMUSHORT *p;
2560   unsigned EMULONG tnum;
2561   UEMUSHORT j, tdenm, tquot;
2562   UEMUSHORT tprod[NI+1];
2563
2564   p = &equot[0];
2565   *p++ = num[0];
2566   *p++ = num[1];
2567
2568   for (i=M; i<NI; i++)
2569     {
2570       *p++ = 0;
2571     }
2572   eshdn1 (num);
2573   tdenm = den[M+1];
2574   for (i=M; i<NI; i++)
2575     {
2576       /* Find trial quotient digit (the radix is 65536).  */
2577       tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2578
2579       /* Do not execute the divide instruction if it will overflow.  */
2580       if ((tdenm * (unsigned long) 0xffff) < tnum)
2581         tquot = 0xffff;
2582       else
2583         tquot = tnum / tdenm;
2584       /* Multiply denominator by trial quotient digit.  */
2585       m16m ((unsigned int) tquot, den, tprod);
2586       /* The quotient digit may have been overestimated.  */
2587       if (ecmpm (tprod, num) > 0)
2588         {
2589           tquot -= 1;
2590           esubm (den, tprod);
2591           if (ecmpm (tprod, num) > 0)
2592             {
2593               tquot -= 1;
2594               esubm (den, tprod);
2595             }
2596         }
2597       esubm (tprod, num);
2598       equot[i] = tquot;
2599       eshup6 (num);
2600     }
2601   /* test for nonzero remainder after roundoff bit */
2602   p = &num[M];
2603   j = 0;
2604   for (i=M; i<NI; i++)
2605     {
2606       j |= *p++;
2607     }
2608   if (j)
2609     j = 1;
2610
2611   for (i=0; i<NI; i++)
2612     num[i] = equot[i];
2613
2614   return ((int) j);
2615 }
2616
2617 /* Multiply significands of exploded e-type A and B, result in B.  */
2618
2619 static int
2620 emulm (a, b)
2621      const UEMUSHORT a[];
2622      UEMUSHORT b[];
2623 {
2624   const UEMUSHORT *p;
2625   UEMUSHORT *q;
2626   UEMUSHORT pprod[NI];
2627   UEMUSHORT j;
2628   int i;
2629
2630   equot[0] = b[0];
2631   equot[1] = b[1];
2632   for (i=M; i<NI; i++)
2633     equot[i] = 0;
2634
2635   j = 0;
2636   p = &a[NI-1];
2637   q = &equot[NI-1];
2638   for (i=M+1; i<NI; i++)
2639     {
2640       if (*p == 0)
2641         {
2642           --p;
2643         }
2644       else
2645         {
2646           m16m ((unsigned int) *p--, b, pprod);
2647           eaddm (pprod, equot);
2648         }
2649       j |= *q;
2650       eshdn6 (equot);
2651     }
2652
2653   for (i=0; i<NI; i++)
2654     b[i] = equot[i];
2655
2656   /* return flag for lost nonzero bits */
2657   return ((int) j);
2658 }
2659 #endif
2660
2661
2662 /* Normalize and round off.
2663
2664   The internal format number to be rounded is S.
2665   Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
2666
2667   Input SUBFLG indicates whether the number was obtained
2668   by a subtraction operation.  In that case if LOST is nonzero
2669   then the number is slightly smaller than indicated.
2670
2671   Input EXP is the biased exponent, which may be negative.
2672   the exponent field of S is ignored but is replaced by
2673   EXP as adjusted by normalization and rounding.
2674
2675   Input RCNTRL is the rounding control.  If it is nonzero, the
2676   returned value will be rounded to RNDPRC bits.
2677
2678   For future reference:  In order for emdnorm to round off denormal
2679    significands at the right point, the input exponent must be
2680    adjusted to be the actual value it would have after conversion to
2681    the final floating point type.  This adjustment has been
2682    implemented for all type conversions (etoe53, etc.) and decimal
2683    conversions, but not for the arithmetic functions (eadd, etc.).
2684    Data types having standard 15-bit exponents are not affected by
2685    this, but SFmode and DFmode are affected. For example, ediv with
2686    rndprc = 24 will not round correctly to 24-bit precision if the
2687    result is denormal.  */
2688
2689 static int rlast = -1;
2690 static int rw = 0;
2691 static UEMUSHORT rmsk = 0;
2692 static UEMUSHORT rmbit = 0;
2693 static UEMUSHORT rebit = 0;
2694 static int re = 0;
2695 static UEMUSHORT rbit[NI];
2696
2697 static void
2698 emdnorm (s, lost, subflg, exp, rcntrl)
2699      UEMUSHORT s[];
2700      int lost;
2701      int subflg ATTRIBUTE_UNUSED;
2702      EMULONG exp;
2703      int rcntrl;
2704 {
2705   int i, j;
2706   UEMUSHORT r;
2707
2708   /* Normalize */
2709   j = enormlz (s);
2710
2711   /* a blank significand could mean either zero or infinity.  */
2712 #ifndef INFINITY
2713   if (j > NBITS)
2714     {
2715       ecleazs (s);
2716       return;
2717     }
2718 #endif
2719   exp -= j;
2720 #ifndef INFINITY
2721   if (exp >= 32767L)
2722     goto overf;
2723 #else
2724   if ((j > NBITS) && (exp < 32767))
2725     {
2726       ecleazs (s);
2727       return;
2728     }
2729 #endif
2730   if (exp < 0L)
2731     {
2732       if (exp > (EMULONG) (-NBITS - 1))
2733         {
2734           j = (int) exp;
2735           i = eshift (s, j);
2736           if (i)
2737             lost = 1;
2738         }
2739       else
2740         {
2741           ecleazs (s);
2742           return;
2743         }
2744     }
2745   /* Round off, unless told not to by rcntrl.  */
2746   if (rcntrl == 0)
2747     goto mdfin;
2748   /* Set up rounding parameters if the control register changed.  */
2749   if (rndprc != rlast)
2750     {
2751       ecleaz (rbit);
2752       switch (rndprc)
2753         {
2754         default:
2755         case NBITS:
2756           rw = NI - 1;          /* low guard word */
2757           rmsk = 0xffff;
2758           rmbit = 0x8000;
2759           re = rw - 1;
2760           rebit = 1;
2761           break;
2762
2763         case 113:
2764           rw = 10;
2765           rmsk = 0x7fff;
2766           rmbit = 0x4000;
2767           rebit = 0x8000;
2768           re = rw;
2769           break;
2770
2771         case 64:
2772           rw = 7;
2773           rmsk = 0xffff;
2774           rmbit = 0x8000;
2775           re = rw - 1;
2776           rebit = 1;
2777           break;
2778
2779           /* For DEC or IBM arithmetic */
2780         case 56:
2781           rw = 6;
2782           rmsk = 0xff;
2783           rmbit = 0x80;
2784           rebit = 0x100;
2785           re = rw;
2786           break;
2787
2788         case 53:
2789           rw = 6;
2790           rmsk = 0x7ff;
2791           rmbit = 0x0400;
2792           rebit = 0x800;
2793           re = rw;
2794           break;
2795
2796           /* For C4x arithmetic */
2797         case 32:
2798           rw = 5;
2799           rmsk = 0xffff;
2800           rmbit = 0x8000;
2801           rebit = 1;
2802           re = rw - 1;
2803           break;
2804
2805         case 24:
2806           rw = 4;
2807           rmsk = 0xff;
2808           rmbit = 0x80;
2809           rebit = 0x100;
2810           re = rw;
2811           break;
2812         }
2813       rbit[re] = rebit;
2814       rlast = rndprc;
2815     }
2816
2817   /* Shift down 1 temporarily if the data structure has an implied
2818      most significant bit and the number is denormal.
2819      Intel long double denormals also lose one bit of precision.  */
2820   if ((exp <= 0) && (rndprc != NBITS)
2821       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2822     {
2823       lost |= s[NI - 1] & 1;
2824       eshdn1 (s);
2825     }
2826   /* Clear out all bits below the rounding bit,
2827      remembering in r if any were nonzero.  */
2828   r = s[rw] & rmsk;
2829   if (rndprc < NBITS)
2830     {
2831       i = rw + 1;
2832       while (i < NI)
2833         {
2834           if (s[i])
2835             r |= 1;
2836           s[i] = 0;
2837           ++i;
2838         }
2839     }
2840   s[rw] &= ~rmsk;
2841   if ((r & rmbit) != 0)
2842     {
2843 #ifndef C4X
2844       if (r == rmbit)
2845         {
2846           if (lost == 0)
2847             {                   /* round to even */
2848               if ((s[re] & rebit) == 0)
2849                 goto mddone;
2850             }
2851           else
2852             {
2853               if (subflg != 0)
2854                 goto mddone;
2855             }
2856         }
2857 #endif
2858       eaddm (rbit, s);
2859     }
2860 #ifndef C4X
2861  mddone:
2862 #endif
2863 /* Undo the temporary shift for denormal values.  */
2864   if ((exp <= 0) && (rndprc != NBITS)
2865       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2866     {
2867       eshup1 (s);
2868     }
2869   if (s[2] != 0)
2870     {                           /* overflow on roundoff */
2871       eshdn1 (s);
2872       exp += 1;
2873     }
2874  mdfin:
2875   s[NI - 1] = 0;
2876   if (exp >= 32767L)
2877     {
2878 #ifndef INFINITY
2879     overf:
2880 #endif
2881 #ifdef INFINITY
2882       s[1] = 32767;
2883       for (i = 2; i < NI - 1; i++)
2884         s[i] = 0;
2885       if (extra_warnings)
2886         warning ("floating point overflow");
2887 #else
2888       s[1] = 32766;
2889       s[2] = 0;
2890       for (i = M + 1; i < NI - 1; i++)
2891         s[i] = 0xffff;
2892       s[NI - 1] = 0;
2893       if ((rndprc < 64) || (rndprc == 113))
2894         {
2895           s[rw] &= ~rmsk;
2896           if (rndprc == 24)
2897             {
2898               s[5] = 0;
2899               s[6] = 0;
2900             }
2901         }
2902 #endif
2903       return;
2904     }
2905   if (exp < 0)
2906     s[1] = 0;
2907   else
2908     s[1] = (UEMUSHORT) exp;
2909 }
2910
2911 /*  Subtract.  C = B - A, all e type numbers.  */
2912
2913 static int subflg = 0;
2914
2915 static void
2916 esub (a, b, c)
2917      const UEMUSHORT *a, *b;
2918      UEMUSHORT *c;
2919 {
2920
2921 #ifdef NANS
2922   if (eisnan (a))
2923     {
2924       emov (a, c);
2925       return;
2926     }
2927   if (eisnan (b))
2928     {
2929       emov (b, c);
2930       return;
2931     }
2932 /* Infinity minus infinity is a NaN.
2933    Test for subtracting infinities of the same sign.  */
2934   if (eisinf (a) && eisinf (b)
2935       && ((eisneg (a) ^ eisneg (b)) == 0))
2936     {
2937       mtherr ("esub", INVALID);
2938       enan (c, 0);
2939       return;
2940     }
2941 #endif
2942   subflg = 1;
2943   eadd1 (a, b, c);
2944 }
2945
2946 /* Add.  C = A + B, all e type.  */
2947
2948 static void
2949 eadd (a, b, c)
2950      const UEMUSHORT *a, *b;
2951      UEMUSHORT *c;
2952 {
2953
2954 #ifdef NANS
2955 /* NaN plus anything is a NaN.  */
2956   if (eisnan (a))
2957     {
2958       emov (a, c);
2959       return;
2960     }
2961   if (eisnan (b))
2962     {
2963       emov (b, c);
2964       return;
2965     }
2966 /* Infinity minus infinity is a NaN.
2967    Test for adding infinities of opposite signs.  */
2968   if (eisinf (a) && eisinf (b)
2969       && ((eisneg (a) ^ eisneg (b)) != 0))
2970     {
2971       mtherr ("esub", INVALID);
2972       enan (c, 0);
2973       return;
2974     }
2975 #endif
2976   subflg = 0;
2977   eadd1 (a, b, c);
2978 }
2979
2980 /* Arithmetic common to both addition and subtraction.  */
2981
2982 static void
2983 eadd1 (a, b, c)
2984      const UEMUSHORT *a, *b;
2985      UEMUSHORT *c;
2986 {
2987   UEMUSHORT ai[NI], bi[NI], ci[NI];
2988   int i, lost, j, k;
2989   EMULONG lt, lta, ltb;
2990
2991 #ifdef INFINITY
2992   if (eisinf (a))
2993     {
2994       emov (a, c);
2995       if (subflg)
2996         eneg (c);
2997       return;
2998     }
2999   if (eisinf (b))
3000     {
3001       emov (b, c);
3002       return;
3003     }
3004 #endif
3005   emovi (a, ai);
3006   emovi (b, bi);
3007   if (subflg)
3008     ai[0] = ~ai[0];
3009
3010   /* compare exponents */
3011   lta = ai[E];
3012   ltb = bi[E];
3013   lt = lta - ltb;
3014   if (lt > 0L)
3015     {                           /* put the larger number in bi */
3016       emovz (bi, ci);
3017       emovz (ai, bi);
3018       emovz (ci, ai);
3019       ltb = bi[E];
3020       lt = -lt;
3021     }
3022   lost = 0;
3023   if (lt != 0L)
3024     {
3025       if (lt < (EMULONG) (-NBITS - 1))
3026         goto done;              /* answer same as larger addend */
3027       k = (int) lt;
3028       lost = eshift (ai, k);    /* shift the smaller number down */
3029     }
3030   else
3031     {
3032       /* exponents were the same, so must compare significands */
3033       i = ecmpm (ai, bi);
3034       if (i == 0)
3035         {                       /* the numbers are identical in magnitude */
3036           /* if different signs, result is zero */
3037           if (ai[0] != bi[0])
3038             {
3039               eclear (c);
3040               return;
3041             }
3042           /* if same sign, result is double */
3043           /* double denormalized tiny number */
3044           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
3045             {
3046               eshup1 (bi);
3047               goto done;
3048             }
3049           /* add 1 to exponent unless both are zero! */
3050           for (j = 1; j < NI - 1; j++)
3051             {
3052               if (bi[j] != 0)
3053                 {
3054                   ltb += 1;
3055                   if (ltb >= 0x7fff)
3056                     {
3057                       eclear (c);
3058                       if (ai[0] != 0)
3059                         eneg (c);
3060                       einfin (c);
3061                       return;
3062                     }
3063                   break;
3064                 }
3065             }
3066           bi[E] = (UEMUSHORT) ltb;
3067           goto done;
3068         }
3069       if (i > 0)
3070         {                       /* put the larger number in bi */
3071           emovz (bi, ci);
3072           emovz (ai, bi);
3073           emovz (ci, ai);
3074         }
3075     }
3076   if (ai[0] == bi[0])
3077     {
3078       eaddm (ai, bi);
3079       subflg = 0;
3080     }
3081   else
3082     {
3083       esubm (ai, bi);
3084       subflg = 1;
3085     }
3086   emdnorm (bi, lost, subflg, ltb, !ROUND_TOWARDS_ZERO);
3087
3088  done:
3089   emovo (bi, c);
3090 }
3091
3092 /* Divide: C = B/A, all e type.  */
3093
3094 static void
3095 ediv (a, b, c)
3096      const UEMUSHORT *a, *b;
3097      UEMUSHORT *c;
3098 {
3099   UEMUSHORT ai[NI], bi[NI];
3100   int i, sign;
3101   EMULONG lt, lta, ltb;
3102
3103 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3104    operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
3105   sign = eisneg (a) ^ eisneg (b);
3106
3107 #ifdef NANS
3108 /* Return any NaN input.  */
3109   if (eisnan (a))
3110     {
3111     emov (a, c);
3112     return;
3113     }
3114   if (eisnan (b))
3115     {
3116     emov (b, c);
3117     return;
3118     }
3119 /* Zero over zero, or infinity over infinity, is a NaN.  */
3120   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
3121       || (eisinf (a) && eisinf (b)))
3122     {
3123     mtherr ("ediv", INVALID);
3124     enan (c, sign);
3125     return;
3126     }
3127 #endif
3128 /* Infinity over anything else is infinity.  */
3129 #ifdef INFINITY
3130   if (eisinf (b))
3131     {
3132       einfin (c);
3133       goto divsign;
3134     }
3135 /* Anything else over infinity is zero.  */
3136   if (eisinf (a))
3137     {
3138       eclear (c);
3139       goto divsign;
3140     }
3141 #endif
3142   emovi (a, ai);
3143   emovi (b, bi);
3144   lta = ai[E];
3145   ltb = bi[E];
3146   if (bi[E] == 0)
3147     {                           /* See if numerator is zero.  */
3148       for (i = 1; i < NI - 1; i++)
3149         {
3150           if (bi[i] != 0)
3151             {
3152               ltb -= enormlz (bi);
3153               goto dnzro1;
3154             }
3155         }
3156       eclear (c);
3157       goto divsign;
3158     }
3159  dnzro1:
3160
3161   if (ai[E] == 0)
3162     {                           /* possible divide by zero */
3163       for (i = 1; i < NI - 1; i++)
3164         {
3165           if (ai[i] != 0)
3166             {
3167               lta -= enormlz (ai);
3168               goto dnzro2;
3169             }
3170         }
3171 /* Divide by zero is not an invalid operation.
3172    It is a divide-by-zero operation!   */
3173       einfin (c);
3174       mtherr ("ediv", SING);
3175       goto divsign;
3176     }
3177  dnzro2:
3178
3179   i = edivm (ai, bi);
3180   /* calculate exponent */
3181   lt = ltb - lta + EXONE;
3182   emdnorm (bi, i, 0, lt, !ROUND_TOWARDS_ZERO);
3183   emovo (bi, c);
3184
3185  divsign:
3186
3187   if (sign
3188 #ifndef IEEE
3189       && (ecmp (c, ezero) != 0)
3190 #endif
3191       )
3192      *(c+(NE-1)) |= 0x8000;
3193   else
3194      *(c+(NE-1)) &= ~0x8000;
3195 }
3196
3197 /* Multiply e-types A and B, return e-type product C.  */
3198
3199 static void
3200 emul (a, b, c)
3201      const UEMUSHORT *a, *b;
3202      UEMUSHORT *c;
3203 {
3204   UEMUSHORT ai[NI], bi[NI];
3205   int i, j, sign;
3206   EMULONG lt, lta, ltb;
3207
3208 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3209    operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
3210   sign = eisneg (a) ^ eisneg (b);
3211
3212 #ifdef NANS
3213 /* NaN times anything is the same NaN.  */
3214   if (eisnan (a))
3215     {
3216     emov (a, c);
3217     return;
3218     }
3219   if (eisnan (b))
3220     {
3221     emov (b, c);
3222     return;
3223     }
3224 /* Zero times infinity is a NaN.  */
3225   if ((eisinf (a) && (ecmp (b, ezero) == 0))
3226       || (eisinf (b) && (ecmp (a, ezero) == 0)))
3227     {
3228     mtherr ("emul", INVALID);
3229     enan (c, sign);
3230     return;
3231     }
3232 #endif
3233 /* Infinity times anything else is infinity.  */
3234 #ifdef INFINITY
3235   if (eisinf (a) || eisinf (b))
3236     {
3237       einfin (c);
3238       goto mulsign;
3239     }
3240 #endif
3241   emovi (a, ai);
3242   emovi (b, bi);
3243   lta = ai[E];
3244   ltb = bi[E];
3245   if (ai[E] == 0)
3246     {
3247       for (i = 1; i < NI - 1; i++)
3248         {
3249           if (ai[i] != 0)
3250             {
3251               lta -= enormlz (ai);
3252               goto mnzer1;
3253             }
3254         }
3255       eclear (c);
3256       goto mulsign;
3257     }
3258  mnzer1:
3259
3260   if (bi[E] == 0)
3261     {
3262       for (i = 1; i < NI - 1; i++)
3263         {
3264           if (bi[i] != 0)
3265             {
3266               ltb -= enormlz (bi);
3267               goto mnzer2;
3268             }
3269         }
3270       eclear (c);
3271       goto mulsign;
3272     }
3273  mnzer2:
3274
3275   /* Multiply significands */
3276   j = emulm (ai, bi);
3277   /* calculate exponent */
3278   lt = lta + ltb - (EXONE - 1);
3279   emdnorm (bi, j, 0, lt, !ROUND_TOWARDS_ZERO);
3280   emovo (bi, c);
3281
3282  mulsign:
3283
3284   if (sign
3285 #ifndef IEEE
3286       && (ecmp (c, ezero) != 0)
3287 #endif
3288       )
3289      *(c+(NE-1)) |= 0x8000;
3290   else
3291      *(c+(NE-1)) &= ~0x8000;
3292 }
3293
3294 /* Convert double precision PE to e-type Y.  */
3295
3296 static void
3297 e53toe (pe, y)
3298      const UEMUSHORT *pe;
3299      UEMUSHORT *y;
3300 {
3301 #ifdef DEC
3302
3303   dectoe (pe, y);
3304
3305 #else
3306 #ifdef IBM
3307
3308   ibmtoe (pe, y, DFmode);
3309
3310 #else
3311 #ifdef C4X
3312
3313   c4xtoe (pe, y, HFmode);
3314
3315 #else
3316
3317   ieeetoe (pe, y, &ieee_53);
3318
3319 #endif /* not C4X */
3320 #endif /* not IBM */
3321 #endif /* not DEC */
3322 }
3323
3324 /* Convert double extended precision float PE to e type Y.  */
3325
3326 static void
3327 e64toe (pe, y)
3328      const UEMUSHORT *pe;
3329      UEMUSHORT *y;
3330 {
3331   UEMUSHORT yy[NI];
3332   const UEMUSHORT *e;
3333   UEMUSHORT *p, *q;
3334   int i;
3335
3336   e = pe;
3337   p = yy;
3338   for (i = 0; i < NE - 5; i++)
3339     *p++ = 0;
3340 #ifndef C4X
3341   /* REAL_WORDS_BIG_ENDIAN is always 0 for DEC and 1 for IBM.
3342      This precision is not ordinarily supported on DEC or IBM.  */
3343   if (! REAL_WORDS_BIG_ENDIAN)
3344     {
3345       for (i = 0; i < 5; i++)
3346         *p++ = *e++;
3347
3348 #ifdef IEEE
3349       /* For denormal long double Intel format, shift significand up one
3350          -- but only if the top significand bit is zero.  A top bit of 1
3351          is "pseudodenormal" when the exponent is zero.  */
3352       if ((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3353         {
3354           UEMUSHORT temp[NI];
3355
3356           emovi (yy, temp);
3357           eshup1 (temp);
3358           emovo (temp,y);
3359           return;
3360         }
3361 #endif /* IEEE */
3362     }
3363   else
3364     {
3365       p = &yy[0] + (NE - 1);
3366       *p-- = *e++;
3367       ++e;
3368       for (i = 0; i < 4; i++)
3369         *p-- = *e++;
3370     }
3371 #endif  /* not C4X */
3372 #ifdef INFINITY
3373   /* Point to the exponent field and check max exponent cases.  */
3374   p = &yy[NE - 1];
3375   if ((*p & 0x7fff) == 0x7fff)
3376     {
3377 #ifdef NANS
3378       if (! REAL_WORDS_BIG_ENDIAN)
3379         {
3380           for (i = 0; i < 4; i++)
3381             {
3382               if ((i != 3 && pe[i] != 0)
3383                   /* Anything but 0x8000 here, including 0, is a NaN.  */
3384                   || (i == 3 && pe[i] != 0x8000))
3385                 {
3386                   enan (y, (*p & 0x8000) != 0);
3387                   return;
3388                 }
3389             }
3390         }
3391       else
3392         {
3393           /* In Motorola extended precision format, the most significant
3394              bit of an infinity mantissa could be either 1 or 0.  It is
3395              the lower order bits that tell whether the value is a NaN.  */
3396           if ((pe[2] & 0x7fff) != 0)
3397             goto bigend_nan;
3398
3399           for (i = 3; i <= 5; i++)
3400             {
3401               if (pe[i] != 0)
3402                 {
3403 bigend_nan:
3404                   enan (y, (*p & 0x8000) != 0);
3405                   return;
3406                 }
3407             }
3408         }
3409 #endif /* NANS */
3410       eclear (y);
3411       einfin (y);
3412       if (*p & 0x8000)
3413         eneg (y);
3414       return;
3415     }
3416 #endif  /* INFINITY */
3417   p = yy;
3418   q = y;
3419   for (i = 0; i < NE; i++)
3420     *q++ = *p++;
3421 }
3422
3423 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3424 /* Convert 128-bit long double precision float PE to e type Y.  */
3425
3426 static void
3427 e113toe (pe, y)
3428      const UEMUSHORT *pe;
3429      UEMUSHORT *y;
3430 {
3431   ieeetoe (pe, y, &ieee_113);
3432 }
3433 #endif  /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3434
3435 /* Convert single precision float PE to e type Y.  */
3436
3437 static void
3438 e24toe (pe, y)
3439      const UEMUSHORT *pe;
3440      UEMUSHORT *y;
3441 {
3442 #ifdef IBM
3443
3444   ibmtoe (pe, y, SFmode);
3445
3446 #else
3447
3448 #ifdef C4X
3449
3450   c4xtoe (pe, y, QFmode);
3451
3452 #else
3453 #ifdef DEC
3454
3455   ieeetoe (pe, y, &dec_f);
3456
3457 #else
3458
3459   ieeetoe (pe, y, &ieee_24);
3460
3461 #endif /* not DEC */
3462 #endif /* not C4X */
3463 #endif /* not IBM */
3464 }
3465
3466 /* Convert machine format float of specified format PE to e type Y.  */
3467
3468 static void
3469 ieeetoe (pe, y, fmt)
3470      const UEMUSHORT *pe;
3471      UEMUSHORT *y;
3472      const struct ieee_format *fmt;
3473 {
3474   UEMUSHORT r;
3475   const UEMUSHORT *e;
3476   UEMUSHORT *p;
3477   UEMUSHORT yy[NI];
3478   int denorm, i, k;
3479   int shortsm1 = fmt->bits / 16 - 1;
3480 #ifdef INFINITY
3481   int expmask = (1 << fmt->expbits) - 1;
3482 #endif
3483 int expshift = (fmt->precision - 1) & 0x0f;
3484   int highbit = 1 << expshift;
3485
3486   e = pe;
3487   denorm = 0;
3488   ecleaz (yy);
3489   if (! REAL_WORDS_BIG_ENDIAN)
3490     e += shortsm1;
3491   r = *e;
3492   yy[0] = 0;
3493   if (r & 0x8000)
3494     yy[0] = 0xffff;
3495   yy[M] = (r & (highbit - 1)) | highbit;
3496   r = (r & 0x7fff) >> expshift;
3497 #ifdef INFINITY
3498   if (!LARGEST_EXPONENT_IS_NORMAL (fmt->precision) && r == expmask)
3499     {
3500 #ifdef NANS
3501       /* First check the word where high order mantissa and exponent live */
3502       if ((*e & (highbit - 1)) != 0)
3503         {
3504           enan (y, yy[0] != 0);
3505           return;
3506         }
3507       if (! REAL_WORDS_BIG_ENDIAN)
3508         {
3509           for (i = 0; i < shortsm1; i++)
3510             {
3511               if (pe[i] != 0)
3512                 {
3513                   enan (y, yy[0] != 0);
3514                   return;
3515                 }
3516             }
3517         }
3518       else
3519         {
3520           for (i = 1; i < shortsm1 + 1; i++)
3521             {
3522               if (pe[i] != 0)
3523                 {
3524                   enan (y, yy[0] != 0);
3525                   return;
3526                 }
3527             }
3528         }
3529 #endif /* NANS */
3530       eclear (y);
3531       einfin (y);
3532       if (yy[0])
3533         eneg (y);
3534       return;
3535     }
3536 #endif  /* INFINITY */
3537   /* If zero exponent, then the significand is denormalized.
3538      So take back the understood high significand bit.  */
3539   if (r == 0)
3540     {
3541       denorm = 1;
3542       yy[M] &= ~highbit;
3543     }
3544   r += fmt->adjustment;
3545   yy[E] = r;
3546   p = &yy[M + 1];
3547   if (! REAL_WORDS_BIG_ENDIAN)
3548     {
3549       for (i = 0; i < shortsm1; i++)
3550         *p++ = *(--e);
3551     }
3552   else
3553     {
3554       ++e;
3555       for (i = 0; i < shortsm1; i++)
3556         *p++ = *e++;
3557     }
3558   if (fmt->precision == 113)
3559     {
3560       /* denorm is left alone in 113 bit format */
3561       if (!denorm)
3562         eshift (yy, -1);
3563     }
3564   else
3565     {
3566       eshift (yy, -(expshift + 1));
3567       if (denorm)
3568         { /* if zero exponent, then normalize the significand */
3569           if ((k = enormlz (yy)) > NBITS)
3570             ecleazs (yy);
3571           else
3572             yy[E] -= (UEMUSHORT) (k - 1);
3573         }
3574     }
3575   emovo (yy, y);
3576 }
3577
3578 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3579 /* Convert e-type X to IEEE 128-bit long double format E.  */
3580
3581 static void
3582 etoe113 (x, e)
3583      const UEMUSHORT *x;
3584      UEMUSHORT *e;
3585 {
3586   etoieee (x, e, &ieee_113);
3587 }
3588
3589 /* Convert exploded e-type X, that has already been rounded to
3590    113-bit precision, to IEEE 128-bit long double format Y.  */
3591
3592 static void
3593 toe113 (x, y)
3594      UEMUSHORT *x, *y;
3595 {
3596   toieee (x, y, &ieee_113);
3597 }
3598
3599 #endif  /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3600
3601 /* Convert e-type X to IEEE double extended format E.  */
3602
3603 static void
3604 etoe64 (x, e)
3605      const UEMUSHORT *x;
3606      UEMUSHORT *e;
3607 {
3608   etoieee (x, e, &ieee_64);
3609 }
3610
3611 /* Convert exploded e-type X, that has already been rounded to
3612    64-bit precision, to IEEE double extended format Y.  */
3613
3614 static void
3615 toe64 (x, y)
3616      UEMUSHORT *x, *y;
3617 {
3618   toieee (x, y, &ieee_64);
3619 }
3620
3621 /* e type to double precision.  */
3622
3623 #ifdef DEC
3624 /* Convert e-type X to DEC-format double E.  */
3625
3626 static void
3627 etoe53 (x, e)
3628      const UEMUSHORT *x;
3629      UEMUSHORT *e;
3630 {
3631   etodec (x, e);
3632 }
3633
3634 /* Convert exploded e-type X, that has already been rounded to
3635    56-bit double precision, to DEC double Y.  */
3636
3637 static void
3638 toe53 (x, y)
3639      UEMUSHORT *x, *y;
3640 {
3641   todec (x, y);
3642 }
3643
3644 #else
3645 #ifdef IBM
3646 /* Convert e-type X to IBM 370-format double E.  */
3647
3648 static void
3649 etoe53 (x, e)
3650      const UEMUSHORT *x;
3651      UEMUSHORT *e;
3652 {
3653   etoibm (x, e, DFmode);
3654 }
3655
3656 /* Convert exploded e-type X, that has already been rounded to
3657    56-bit precision, to IBM 370 double Y.  */
3658
3659 static void
3660 toe53 (x, y)
3661      UEMUSHORT *x, *y;
3662 {
3663   toibm (x, y, DFmode);
3664 }
3665
3666 #else /* it's neither DEC nor IBM */
3667 #ifdef C4X
3668 /* Convert e-type X to C4X-format long double E.  */
3669
3670 static void
3671 etoe53 (x, e)
3672      const UEMUSHORT *x;
3673      UEMUSHORT *e;
3674 {
3675   etoc4x (x, e, HFmode);
3676 }
3677
3678 /* Convert exploded e-type X, that has already been rounded to
3679    56-bit precision, to IBM 370 double Y.  */
3680
3681 static void
3682 toe53 (x, y)
3683      UEMUSHORT *x, *y;
3684 {
3685   toc4x (x, y, HFmode);
3686 }
3687
3688 #else  /* it's neither DEC nor IBM nor C4X */
3689
3690 /* Convert e-type X to IEEE double E.  */
3691
3692 static void
3693 etoe53 (x, e)
3694       const UEMUSHORT *x;
3695       UEMUSHORT *e;
3696 {
3697   etoieee (x, e, &ieee_53);
3698 }
3699
3700 /* Convert exploded e-type X, that has already been rounded to
3701    53-bit precision, to IEEE double Y.  */
3702
3703 static void
3704 toe53 (x, y)
3705      UEMUSHORT *x, *y;
3706 {
3707   toieee (x, y, &ieee_53);
3708 }
3709
3710 #endif /* not C4X */
3711 #endif /* not IBM */
3712 #endif /* not DEC */
3713
3714
3715
3716 /* e type to single precision.  */
3717
3718 #ifdef IBM
3719 /* Convert e-type X to IBM 370 float E.  */
3720
3721 static void
3722 etoe24 (x, e)
3723      const UEMUSHORT *x;
3724      UEMUSHORT *e;
3725 {
3726   etoibm (x, e, SFmode);
3727 }
3728
3729 /* Convert exploded e-type X, that has already been rounded to
3730    float precision, to IBM 370 float Y.  */
3731
3732 static void
3733 toe24 (x, y)
3734      UEMUSHORT *x, *y;
3735 {
3736   toibm (x, y, SFmode);
3737 }
3738
3739 #else /* it's not IBM */
3740
3741 #ifdef C4X
3742 /* Convert e-type X to C4X float E.  */
3743
3744 static void
3745 etoe24 (x, e)
3746      const UEMUSHORT *x;
3747      UEMUSHORT *e;
3748 {
3749   etoc4x (x, e, QFmode);
3750 }
3751
3752 /* Convert exploded e-type X, that has already been rounded to
3753    float precision, to IBM 370 float Y.  */
3754
3755 static void
3756 toe24 (x, y)
3757      UEMUSHORT *x, *y;
3758 {
3759   toc4x (x, y, QFmode);
3760 }
3761
3762 #else /* it's neither IBM nor C4X */
3763
3764 #ifdef DEC
3765
3766 /* Convert e-type X to DEC F-float E.  */
3767
3768 static void
3769 etoe24 (x, e)
3770       const UEMUSHORT *x;
3771       UEMUSHORT *e;
3772 {
3773   etoieee (x, e, &dec_f);
3774 }
3775
3776 /* Convert exploded e-type X, that has already been rounded to
3777    float precision, to DEC F-float Y.  */
3778
3779 static void
3780 toe24 (x, y)
3781       UEMUSHORT *x, *y;
3782 {
3783   toieee (x, y, &dec_f);
3784 }
3785
3786 #else
3787
3788 /* Convert e-type X to IEEE float E.  */
3789
3790 static void
3791 etoe24 (x, e)
3792       const UEMUSHORT *x;
3793       UEMUSHORT *e;
3794 {
3795   etoieee (x, e, &ieee_24);
3796 }
3797
3798 /* Convert exploded e-type X, that has already been rounded to
3799    float precision, to IEEE float Y.  */
3800
3801 static void
3802 toe24 (x, y)
3803       UEMUSHORT *x, *y;
3804 {
3805   toieee (x, y, &ieee_24);
3806 }
3807
3808 #endif  /* not DEC */
3809 #endif  /* not C4X */
3810 #endif  /* not IBM */
3811
3812
3813 /* Convert e-type X to the IEEE format described by FMT.  */
3814
3815 static void
3816 etoieee (x, e, fmt)
3817       const UEMUSHORT *x;
3818       UEMUSHORT *e;
3819       const struct ieee_format *fmt;
3820 {
3821   UEMUSHORT xi[NI];
3822   EMULONG exp;
3823   int rndsav;
3824
3825 #ifdef NANS
3826   if (eisnan (x))
3827     {
3828       make_nan (e, eisneg (x), fmt->mode);
3829       return;
3830     }
3831 #endif
3832
3833   emovi (x, xi);
3834
3835 #ifdef INFINITY
3836   if (eisinf (x))
3837     goto nonorm;
3838 #endif
3839   /* Adjust exponent for offset.  */
3840   exp = (EMULONG) xi[E] - fmt->adjustment;
3841
3842   /* Round off to nearest or even.  */
3843   rndsav = rndprc;
3844   rndprc = fmt->precision;
3845   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3846   rndprc = rndsav;
3847 #ifdef INFINITY
3848  nonorm:
3849 #endif
3850   toieee (xi, e, fmt);
3851 }
3852
3853 /* Convert exploded e-type X, that has already been rounded to
3854    the necessary precision, to the IEEE format described by FMT.  */
3855
3856 static void
3857 toieee (x, y, fmt)
3858       UEMUSHORT *x, *y;
3859       const struct ieee_format *fmt;
3860 {
3861   UEMUSHORT maxexp;
3862   UEMUSHORT *q;
3863   int words;
3864   int i;
3865
3866   maxexp = (1 << fmt->expbits) - 1;
3867   words = (fmt->bits - fmt->expbits) / EMUSHORT_SIZE;
3868
3869 #ifdef NANS
3870   if (eiisnan (x))
3871     {
3872       make_nan (y, eiisneg (x), fmt->mode);
3873       return;
3874     }
3875 #endif
3876
3877   if (fmt->expbits < 15
3878       && LARGEST_EXPONENT_IS_NORMAL (fmt->bits)
3879       && x[E] > maxexp)
3880     {
3881       saturate (y, eiisneg (x), fmt->bits, 1);
3882       return;
3883     }
3884
3885   /* Point to the exponent.  */
3886   if (REAL_WORDS_BIG_ENDIAN)
3887     q = y;
3888   else
3889     q = y + words;
3890
3891   /* Copy the sign.  */
3892   if (x[0])
3893     *q = 0x8000;
3894   else
3895     *q = 0;
3896
3897   if (fmt->expbits < 15
3898       && !LARGEST_EXPONENT_IS_NORMAL (fmt->bits)
3899       && x[E] >= maxexp)
3900     {
3901       /* Saturate at largest number less that infinity.  */
3902       UEMUSHORT fill;
3903 #ifdef INFINITY
3904       *q |= maxexp << (15 - fmt->expbits);
3905       fill = 0;
3906 #else
3907       *q |= (maxexp << (15 - fmt->expbits)) - 1;
3908       fill = 0xffff;
3909 #endif
3910
3911       if (!REAL_WORDS_BIG_ENDIAN)
3912         {
3913           for (i = 0; i < words; i++)
3914             *(--q) = fill;
3915         }
3916       else
3917         {
3918           for (i = 0; i < words; i++)
3919             *(++q) = fill;
3920         }
3921 #if defined(INFINITY) && defined(ERANGE)
3922       errno = ERANGE;
3923 #endif
3924       return;
3925     }
3926
3927   /* If denormal and DEC float, return zero (DEC has no denormals) */
3928 #ifdef DEC
3929   if (x[E] == 0)
3930     {
3931       for (i = 0; i < fmt->bits / EMUSHORT_SIZE ; i++)
3932         q[i] = 0;
3933       return;
3934     }
3935 #endif /* DEC */
3936
3937   /* Delete the implied bit unless denormal, except for
3938      64-bit precision.  */
3939   if (fmt->precision != 64 && x[E] != 0)
3940     {
3941       eshup1 (x);
3942     }
3943
3944   /* Shift denormal double extended Intel format significand down
3945      one bit.  */
3946   if (fmt->precision == 64 && x[E] == 0 && ! REAL_WORDS_BIG_ENDIAN)
3947     eshdn1 (x);
3948
3949   if (fmt->expbits < 15)
3950     {
3951       /* Shift the significand.  */
3952       eshift (x, 15 - fmt->expbits);
3953
3954       /* Combine the exponent and upper bits of the significand.  */
3955       *q |= x[E] << (15 - fmt->expbits);
3956       *q |= x[M] & (UEMUSHORT) ~((maxexp << (15 - fmt->expbits)) | 0x8000);
3957     }
3958   else
3959     {
3960       /* Copy the exponent.  */
3961       *q |= x[E];
3962     }
3963
3964   /* Add padding after the exponent. At the moment, this is only necessary for
3965      64-bit precision; in this case, the padding is 16 bits.  */
3966   if (fmt->precision == 64)
3967     {
3968       *(q + 1) = 0;
3969
3970       /* Skip padding.  */
3971       if (REAL_WORDS_BIG_ENDIAN)
3972         ++q;
3973     }
3974
3975   /* Copy the significand.  */
3976   if (REAL_WORDS_BIG_ENDIAN)
3977     {
3978       for (i = 0; i < words; i++)
3979         *(++q) = x[i + M + 1];
3980     }
3981 #ifdef INFINITY
3982   else if (fmt->precision == 64 && eiisinf (x))
3983     {
3984       /* Intel double extended infinity significand.  */
3985       *(--q) = 0x8000;
3986       *(--q) = 0;
3987       *(--q) = 0;
3988       *(--q) = 0;
3989     }
3990 #endif
3991   else
3992     {
3993       for (i = 0; i < words; i++)
3994         *(--q) = x[i + M + 1];
3995     }
3996 }
3997
3998
3999 /* Compare two e type numbers.
4000    Return +1 if a > b
4001            0 if a == b
4002           -1 if a < b
4003           -2 if either a or b is a NaN.  */
4004
4005 static int
4006 ecmp (a, b)
4007      const UEMUSHORT *a, *b;
4008 {
4009   UEMUSHORT ai[NI], bi[NI];
4010   UEMUSHORT *p, *q;
4011   int i;
4012   int msign;
4013
4014 #ifdef NANS
4015   if (eisnan (a)  || eisnan (b))
4016       return (-2);
4017 #endif
4018   emovi (a, ai);
4019   p = ai;
4020   emovi (b, bi);
4021   q = bi;
4022
4023   if (*p != *q)
4024     {                           /* the signs are different */
4025       /* -0 equals + 0 */
4026       for (i = 1; i < NI - 1; i++)
4027         {
4028           if (ai[i] != 0)
4029             goto nzro;
4030           if (bi[i] != 0)
4031             goto nzro;
4032         }
4033       return (0);
4034     nzro:
4035       if (*p == 0)
4036         return (1);
4037       else
4038         return (-1);
4039     }
4040   /* both are the same sign */
4041   if (*p == 0)
4042     msign = 1;
4043   else
4044     msign = -1;
4045   i = NI - 1;
4046   do
4047     {
4048       if (*p++ != *q++)
4049         {
4050           goto diff;
4051         }
4052     }
4053   while (--i > 0);
4054
4055   return (0);                   /* equality */
4056
4057  diff:
4058
4059   if (*(--p) > *(--q))
4060     return (msign);             /* p is bigger */
4061   else
4062     return (-msign);            /* p is littler */
4063 }
4064
4065 #if 0
4066 /* Find e-type nearest integer to X, as floor (X + 0.5).  */
4067
4068 static void
4069 eround (x, y)
4070      const UEMUSHORT *x;
4071      UEMUSHORT *y;
4072 {
4073   eadd (ehalf, x, y);
4074   efloor (y, y);
4075 }
4076 #endif /* 0 */
4077
4078 /* Convert HOST_WIDE_INT LP to e type Y.  */
4079
4080 static void
4081 ltoe (lp, y)
4082      const HOST_WIDE_INT *lp;
4083      UEMUSHORT *y;
4084 {
4085   UEMUSHORT yi[NI];
4086   unsigned HOST_WIDE_INT ll;
4087   int k;
4088
4089   ecleaz (yi);
4090   if (*lp < 0)
4091     {
4092       /* make it positive */
4093       ll = (unsigned HOST_WIDE_INT) (-(*lp));
4094       yi[0] = 0xffff;           /* put correct sign in the e type number */
4095     }
4096   else
4097     {
4098       ll = (unsigned HOST_WIDE_INT) (*lp);
4099     }
4100   /* move the long integer to yi significand area */
4101 #if HOST_BITS_PER_WIDE_INT == 64
4102   yi[M] = (UEMUSHORT) (ll >> 48);
4103   yi[M + 1] = (UEMUSHORT) (ll >> 32);
4104   yi[M + 2] = (UEMUSHORT) (ll >> 16);
4105   yi[M + 3] = (UEMUSHORT) ll;
4106   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
4107 #else
4108   yi[M] = (UEMUSHORT) (ll >> 16);
4109   yi[M + 1] = (UEMUSHORT) ll;
4110   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4111 #endif
4112
4113   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4114     ecleaz (yi);                /* it was zero */
4115   else
4116     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
4117   emovo (yi, y);                /* output the answer */
4118 }
4119
4120 /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
4121
4122 static void
4123 ultoe (lp, y)
4124      const unsigned HOST_WIDE_INT *lp;
4125      UEMUSHORT *y;
4126 {
4127   UEMUSHORT yi[NI];
4128   unsigned HOST_WIDE_INT ll;
4129   int k;
4130
4131   ecleaz (yi);
4132   ll = *lp;
4133
4134   /* move the long integer to ayi significand area */
4135 #if HOST_BITS_PER_WIDE_INT == 64
4136   yi[M] = (UEMUSHORT) (ll >> 48);
4137   yi[M + 1] = (UEMUSHORT) (ll >> 32);
4138   yi[M + 2] = (UEMUSHORT) (ll >> 16);
4139   yi[M + 3] = (UEMUSHORT) ll;
4140   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
4141 #else
4142   yi[M] = (UEMUSHORT) (ll >> 16);
4143   yi[M + 1] = (UEMUSHORT) ll;
4144   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4145 #endif
4146
4147   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4148     ecleaz (yi);                /* it was zero */
4149   else
4150     yi[E] -= (UEMUSHORT) k;  /* subtract shift count from exponent */
4151   emovo (yi, y);                /* output the answer */
4152 }
4153
4154
4155 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4156    part FRAC of e-type (packed internal format) floating point input X.
4157    The integer output I has the sign of the input, except that
4158    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4159    The output e-type fraction FRAC is the positive fractional
4160    part of abs (X).  */
4161
4162 static void
4163 eifrac (x, i, frac)
4164      const UEMUSHORT *x;
4165      HOST_WIDE_INT *i;
4166      UEMUSHORT *frac;
4167 {
4168   UEMUSHORT xi[NI];
4169   int j, k;
4170   unsigned HOST_WIDE_INT ll;
4171
4172   emovi (x, xi);
4173   k = (int) xi[E] - (EXONE - 1);
4174   if (k <= 0)
4175     {
4176       /* if exponent <= 0, integer = 0 and real output is fraction */
4177       *i = 0L;
4178       emovo (xi, frac);
4179       return;
4180     }
4181   if (k > (HOST_BITS_PER_WIDE_INT - 1))
4182     {
4183       /* long integer overflow: output large integer
4184          and correct fraction  */
4185       if (xi[0])
4186         *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4187       else
4188         {
4189 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4190           /* In this case, let it overflow and convert as if unsigned.  */
4191           euifrac (x, &ll, frac);
4192           *i = (HOST_WIDE_INT) ll;
4193           return;
4194 #else
4195           /* In other cases, return the largest positive integer.  */
4196           *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4197 #endif
4198         }
4199       eshift (xi, k);
4200       if (extra_warnings)
4201         warning ("overflow on truncation to integer");
4202     }
4203   else if (k > 16)
4204     {
4205       /* Shift more than 16 bits: first shift up k-16 mod 16,
4206          then shift up by 16's.  */
4207       j = k - ((k >> 4) << 4);
4208       eshift (xi, j);
4209       ll = xi[M];
4210       k -= j;
4211       do
4212         {
4213           eshup6 (xi);
4214           ll = (ll << 16) | xi[M];
4215         }
4216       while ((k -= 16) > 0);
4217       *i = ll;
4218       if (xi[0])
4219         *i = -(*i);
4220     }
4221   else
4222     {
4223       /* shift not more than 16 bits */
4224       eshift (xi, k);
4225       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4226       if (xi[0])
4227         *i = -(*i);
4228     }
4229   xi[0] = 0;
4230   xi[E] = EXONE - 1;
4231   xi[M] = 0;
4232   if ((k = enormlz (xi)) > NBITS)
4233     ecleaz (xi);
4234   else
4235     xi[E] -= (UEMUSHORT) k;
4236
4237   emovo (xi, frac);
4238 }
4239
4240
4241 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4242    FRAC of e-type X.  A negative input yields integer output = 0 but
4243    correct fraction.  */
4244
4245 static void
4246 euifrac (x, i, frac)
4247      const UEMUSHORT *x;
4248      unsigned HOST_WIDE_INT *i;
4249      UEMUSHORT *frac;
4250 {
4251   unsigned HOST_WIDE_INT ll;
4252   UEMUSHORT xi[NI];
4253   int j, k;
4254
4255   emovi (x, xi);
4256   k = (int) xi[E] - (EXONE - 1);
4257   if (k <= 0)
4258     {
4259       /* if exponent <= 0, integer = 0 and argument is fraction */
4260       *i = 0L;
4261       emovo (xi, frac);
4262       return;
4263     }
4264   if (k > HOST_BITS_PER_WIDE_INT)
4265     {
4266       /* Long integer overflow: output large integer
4267          and correct fraction.
4268          Note, the BSD MicroVAX compiler says that ~(0UL)
4269          is a syntax error.  */
4270       *i = ~(0L);
4271       eshift (xi, k);
4272       if (extra_warnings)
4273         warning ("overflow on truncation to unsigned integer");
4274     }
4275   else if (k > 16)
4276     {
4277       /* Shift more than 16 bits: first shift up k-16 mod 16,
4278          then shift up by 16's.  */
4279       j = k - ((k >> 4) << 4);
4280       eshift (xi, j);
4281       ll = xi[M];
4282       k -= j;
4283       do
4284         {
4285           eshup6 (xi);
4286           ll = (ll << 16) | xi[M];
4287         }
4288       while ((k -= 16) > 0);
4289       *i = ll;
4290     }
4291   else
4292     {
4293       /* shift not more than 16 bits */
4294       eshift (xi, k);
4295       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4296     }
4297
4298   if (xi[0])  /* A negative value yields unsigned integer 0.  */
4299     *i = 0L;
4300
4301   xi[0] = 0;
4302   xi[E] = EXONE - 1;
4303   xi[M] = 0;
4304   if ((k = enormlz (xi)) > NBITS)
4305     ecleaz (xi);
4306   else
4307     xi[E] -= (UEMUSHORT) k;
4308
4309   emovo (xi, frac);
4310 }
4311
4312 /* Shift the significand of exploded e-type X up or down by SC bits.  */
4313
4314 static int
4315 eshift (x, sc)
4316      UEMUSHORT *x;
4317      int sc;
4318 {
4319   UEMUSHORT lost;
4320   UEMUSHORT *p;
4321
4322   if (sc == 0)
4323     return (0);
4324
4325   lost = 0;
4326   p = x + NI - 1;
4327
4328   if (sc < 0)
4329     {
4330       sc = -sc;
4331       while (sc >= 16)
4332         {
4333           lost |= *p;           /* remember lost bits */
4334           eshdn6 (x);
4335           sc -= 16;
4336         }
4337
4338       while (sc >= 8)
4339         {
4340           lost |= *p & 0xff;
4341           eshdn8 (x);
4342           sc -= 8;
4343         }
4344
4345       while (sc > 0)
4346         {
4347           lost |= *p & 1;
4348           eshdn1 (x);
4349           sc -= 1;
4350         }
4351     }
4352   else
4353     {
4354       while (sc >= 16)
4355         {
4356           eshup6 (x);
4357           sc -= 16;
4358         }
4359
4360       while (sc >= 8)
4361         {
4362           eshup8 (x);
4363           sc -= 8;
4364         }
4365
4366       while (sc > 0)
4367         {
4368           eshup1 (x);
4369           sc -= 1;
4370         }
4371     }
4372   if (lost)
4373     lost = 1;
4374   return ((int) lost);
4375 }
4376
4377 /* Shift normalize the significand area of exploded e-type X.
4378    Return the shift count (up = positive).  */
4379
4380 static int
4381 enormlz (x)
4382      UEMUSHORT x[];
4383 {
4384   UEMUSHORT *p;
4385   int sc;
4386
4387   sc = 0;
4388   p = &x[M];
4389   if (*p != 0)
4390     goto normdn;
4391   ++p;
4392   if (*p & 0x8000)
4393     return (0);                 /* already normalized */
4394   while (*p == 0)
4395     {
4396       eshup6 (x);
4397       sc += 16;
4398
4399       /* With guard word, there are NBITS+16 bits available.
4400        Return true if all are zero.  */
4401       if (sc > NBITS)
4402         return (sc);
4403     }
4404   /* see if high byte is zero */
4405   while ((*p & 0xff00) == 0)
4406     {
4407       eshup8 (x);
4408       sc += 8;
4409     }
4410   /* now shift 1 bit at a time */
4411   while ((*p & 0x8000) == 0)
4412     {
4413       eshup1 (x);
4414       sc += 1;
4415       if (sc > NBITS)
4416         {
4417           mtherr ("enormlz", UNDERFLOW);
4418           return (sc);
4419         }
4420     }
4421   return (sc);
4422
4423   /* Normalize by shifting down out of the high guard word
4424      of the significand */
4425  normdn:
4426
4427   if (*p & 0xff00)
4428     {
4429       eshdn8 (x);
4430       sc -= 8;
4431     }
4432   while (*p != 0)
4433     {
4434       eshdn1 (x);
4435       sc -= 1;
4436
4437       if (sc < -NBITS)
4438         {
4439           mtherr ("enormlz", OVERFLOW);
4440           return (sc);
4441         }
4442     }
4443   return (sc);
4444 }
4445
4446 /* Powers of ten used in decimal <-> binary conversions.  */
4447
4448 #define NTEN 12
4449 #define MAXP 4096
4450
4451 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4452 static const UEMUSHORT etens[NTEN + 1][NE] =
4453 {
4454   {0x6576, 0x4a92, 0x804a, 0x153f,
4455    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4456   {0x6a32, 0xce52, 0x329a, 0x28ce,
4457    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4458   {0x526c, 0x50ce, 0xf18b, 0x3d28,
4459    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4460   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4461    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4462   {0x851e, 0xeab7, 0x98fe, 0x901b,
4463    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4464   {0x0235, 0x0137, 0x36b1, 0x336c,
4465    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4466   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4467    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4468   {0x0000, 0x0000, 0x0000, 0x0000,
4469    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4470   {0x0000, 0x0000, 0x0000, 0x0000,
4471    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4472   {0x0000, 0x0000, 0x0000, 0x0000,
4473    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4474   {0x0000, 0x0000, 0x0000, 0x0000,
4475    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4476   {0x0000, 0x0000, 0x0000, 0x0000,
4477    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4478   {0x0000, 0x0000, 0x0000, 0x0000,
4479    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4480 };
4481
4482 static const UEMUSHORT emtens[NTEN + 1][NE] =
4483 {
4484   {0x2030, 0xcffc, 0xa1c3, 0x8123,
4485    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4486   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4487    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4488   {0xf53f, 0xf698, 0x6bd3, 0x0158,
4489    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4490   {0xe731, 0x04d4, 0xe3f2, 0xd332,
4491    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4492   {0xa23e, 0x5308, 0xfefb, 0x1155,
4493    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4494   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4495    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4496   {0x2a20, 0x6224, 0x47b3, 0x98d7,
4497    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4498   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4499    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4500   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4501    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4502   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4503    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4504   {0xc155, 0xa4a8, 0x404e, 0x6113,
4505    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4506   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4507    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4508   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4509    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4510 };
4511 #else
4512 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4513 static const UEMUSHORT etens[NTEN + 1][NE] =
4514 {
4515   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4516   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4517   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4518   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4519   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4520   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4521   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4522   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4523   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4524   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4525   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4526   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4527   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4528 };
4529
4530 static const UEMUSHORT emtens[NTEN + 1][NE] =
4531 {
4532   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4533   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4534   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4535   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4536   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4537   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4538   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4539   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4540   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4541   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4542   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4543   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4544   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4545 };
4546 #endif
4547
4548 #if 0
4549 /* Convert float value X to ASCII string STRING with NDIG digits after
4550    the decimal point.  */
4551
4552 static void
4553 e24toasc (x, string, ndigs)
4554      const UEMUSHORT x[];
4555      char *string;
4556      int ndigs;
4557 {
4558   UEMUSHORT w[NI];
4559
4560   e24toe (x, w);
4561   etoasc (w, string, ndigs);
4562 }
4563
4564 /* Convert double value X to ASCII string STRING with NDIG digits after
4565    the decimal point.  */
4566
4567 static void
4568 e53toasc (x, string, ndigs)
4569      const UEMUSHORT x[];
4570      char *string;
4571      int ndigs;
4572 {
4573   UEMUSHORT w[NI];
4574
4575   e53toe (x, w);
4576   etoasc (w, string, ndigs);
4577 }
4578
4579 /* Convert double extended value X to ASCII string STRING with NDIG digits
4580    after the decimal point.  */
4581
4582 static void
4583 e64toasc (x, string, ndigs)
4584      const UEMUSHORT x[];
4585      char *string;
4586      int ndigs;
4587 {
4588   UEMUSHORT w[NI];
4589
4590   e64toe (x, w);
4591   etoasc (w, string, ndigs);
4592 }
4593
4594 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4595    after the decimal point.  */
4596
4597 static void
4598 e113toasc (x, string, ndigs)
4599      const UEMUSHORT x[];
4600      char *string;
4601      int ndigs;
4602 {
4603   UEMUSHORT w[NI];
4604
4605   e113toe (x, w);
4606   etoasc (w, string, ndigs);
4607 }
4608 #endif /* 0 */
4609
4610 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4611    the decimal point.  */
4612
4613 static char wstring[80];        /* working storage for ASCII output */
4614
4615 static void
4616 etoasc (x, string, ndigs)
4617      const UEMUSHORT x[];
4618      char *string;
4619      int ndigs;
4620 {
4621   EMUSHORT digit;
4622   UEMUSHORT y[NI], t[NI], u[NI], w[NI];
4623   const UEMUSHORT *p, *r, *ten;
4624   UEMUSHORT sign;
4625   int i, j, k, expon, rndsav;
4626   char *s, *ss;
4627   UEMUSHORT m;
4628
4629
4630   rndsav = rndprc;
4631   ss = string;
4632   s = wstring;
4633   *ss = '\0';
4634   *s = '\0';
4635 #ifdef NANS
4636   if (eisnan (x))
4637     {
4638       sprintf (wstring, " NaN ");
4639       goto bxit;
4640     }
4641 #endif
4642   rndprc = NBITS;               /* set to full precision */
4643   emov (x, y);                  /* retain external format */
4644   if (y[NE - 1] & 0x8000)
4645     {
4646       sign = 0xffff;
4647       y[NE - 1] &= 0x7fff;
4648     }
4649   else
4650     {
4651       sign = 0;
4652     }
4653   expon = 0;
4654   ten = &etens[NTEN][0];
4655   emov (eone, t);
4656   /* Test for zero exponent */
4657   if (y[NE - 1] == 0)
4658     {
4659       for (k = 0; k < NE - 1; k++)
4660         {
4661           if (y[k] != 0)
4662             goto tnzro;         /* denormalized number */
4663         }
4664       goto isone;               /* valid all zeros */
4665     }
4666  tnzro:
4667
4668   /* Test for infinity.  */
4669   if (y[NE - 1] == 0x7fff)
4670     {
4671       if (sign)
4672         sprintf (wstring, " -Infinity ");
4673       else
4674         sprintf (wstring, " Infinity ");
4675       goto bxit;
4676     }
4677
4678   /* Test for exponent nonzero but significand denormalized.
4679    * This is an error condition.
4680    */
4681   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4682     {
4683       mtherr ("etoasc", DOMAIN);
4684       sprintf (wstring, "NaN");
4685       goto bxit;
4686     }
4687
4688   /* Compare to 1.0 */
4689   i = ecmp (eone, y);
4690   if (i == 0)
4691     goto isone;
4692
4693   if (i == -2)
4694     abort ();
4695
4696   if (i < 0)
4697     {                           /* Number is greater than 1 */
4698       /* Convert significand to an integer and strip trailing decimal zeros.  */
4699       emov (y, u);
4700       u[NE - 1] = EXONE + NBITS - 1;
4701
4702       p = &etens[NTEN - 4][0];
4703       m = 16;
4704       do
4705         {
4706           ediv (p, u, t);
4707           efloor (t, w);
4708           for (j = 0; j < NE - 1; j++)
4709             {
4710               if (t[j] != w[j])
4711                 goto noint;
4712             }
4713           emov (t, u);
4714           expon += (int) m;
4715         noint:
4716           p += NE;
4717           m >>= 1;
4718         }
4719       while (m != 0);
4720
4721       /* Rescale from integer significand */
4722       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4723       emov (u, y);
4724       /* Find power of 10 */
4725       emov (eone, t);
4726       m = MAXP;
4727       p = &etens[0][0];
4728       /* An unordered compare result shouldn't happen here.  */
4729       while (ecmp (ten, u) <= 0)
4730         {
4731           if (ecmp (p, u) <= 0)
4732             {
4733               ediv (p, u, u);
4734               emul (p, t, t);
4735               expon += (int) m;
4736             }
4737           m >>= 1;
4738           if (m == 0)
4739             break;
4740           p += NE;
4741         }
4742     }
4743   else
4744     {                           /* Number is less than 1.0 */
4745       /* Pad significand with trailing decimal zeros.  */
4746       if (y[NE - 1] == 0)
4747         {
4748           while ((y[NE - 2] & 0x8000) == 0)
4749             {
4750               emul (ten, y, y);
4751               expon -= 1;
4752             }
4753         }
4754       else
4755         {
4756           emovi (y, w);
4757           for (i = 0; i < NDEC + 1; i++)
4758             {
4759               if ((w[NI - 1] & 0x7) != 0)
4760                 break;
4761               /* multiply by 10 */
4762               emovz (w, u);
4763               eshdn1 (u);
4764               eshdn1 (u);
4765               eaddm (w, u);
4766               u[1] += 3;
4767               while (u[2] != 0)
4768                 {
4769                   eshdn1 (u);
4770                   u[1] += 1;
4771                 }
4772               if (u[NI - 1] != 0)
4773                 break;
4774               if (eone[NE - 1] <= u[1])
4775                 break;
4776               emovz (u, w);
4777               expon -= 1;
4778             }
4779           emovo (w, y);
4780         }
4781       k = -MAXP;
4782       p = &emtens[0][0];
4783       r = &etens[0][0];
4784       emov (y, w);
4785       emov (eone, t);
4786       while (ecmp (eone, w) > 0)
4787         {
4788           if (ecmp (p, w) >= 0)
4789             {
4790               emul (r, w, w);
4791               emul (r, t, t);
4792               expon += k;
4793             }
4794           k /= 2;
4795           if (k == 0)
4796             break;
4797           p += NE;
4798           r += NE;
4799         }
4800       ediv (t, eone, t);
4801     }
4802  isone:
4803   /* Find the first (leading) digit.  */
4804   emovi (t, w);
4805   emovz (w, t);
4806   emovi (y, w);
4807   emovz (w, y);
4808   eiremain (t, y);
4809   digit = equot[NI - 1];
4810   while ((digit == 0) && (ecmp (y, ezero) != 0))
4811     {
4812       eshup1 (y);
4813       emovz (y, u);
4814       eshup1 (u);
4815       eshup1 (u);
4816       eaddm (u, y);
4817       eiremain (t, y);
4818       digit = equot[NI - 1];
4819       expon -= 1;
4820     }
4821   s = wstring;
4822   if (sign)
4823     *s++ = '-';
4824   else
4825     *s++ = ' ';
4826   /* Examine number of digits requested by caller.  */
4827   if (ndigs < 0)
4828     ndigs = 0;
4829   if (ndigs > NDEC)
4830     ndigs = NDEC;
4831   if (digit == 10)
4832     {
4833       *s++ = '1';
4834       *s++ = '.';
4835       if (ndigs > 0)
4836         {
4837           *s++ = '0';
4838           ndigs -= 1;
4839         }
4840       expon += 1;
4841     }
4842   else
4843     {
4844       *s++ = (char) digit + '0';
4845       *s++ = '.';
4846     }
4847   /* Generate digits after the decimal point.  */
4848   for (k = 0; k <= ndigs; k++)
4849     {
4850       /* multiply current number by 10, without normalizing */
4851       eshup1 (y);
4852       emovz (y, u);
4853       eshup1 (u);
4854       eshup1 (u);
4855       eaddm (u, y);
4856       eiremain (t, y);
4857       *s++ = (char) equot[NI - 1] + '0';
4858     }
4859   digit = equot[NI - 1];
4860   --s;
4861   ss = s;
4862   /* round off the ASCII string */
4863   if (digit > 4)
4864     {
4865       /* Test for critical rounding case in ASCII output.  */
4866       if (digit == 5)
4867         {
4868           emovo (y, t);
4869           if (ecmp (t, ezero) != 0)
4870             goto roun;          /* round to nearest */
4871 #ifndef C4X
4872           if ((*(s - 1) & 1) == 0)
4873             goto doexp;         /* round to even */
4874 #endif
4875         }
4876       /* Round up and propagate carry-outs */
4877     roun:
4878       --s;
4879       k = *s & CHARMASK;
4880       /* Carry out to most significant digit? */
4881       if (k == '.')
4882         {
4883           --s;
4884           k = *s;
4885           k += 1;
4886           *s = (char) k;
4887           /* Most significant digit carries to 10? */
4888           if (k > '9')
4889             {
4890               expon += 1;
4891               *s = '1';
4892             }
4893           goto doexp;
4894         }
4895       /* Round up and carry out from less significant digits */
4896       k += 1;
4897       *s = (char) k;
4898       if (k > '9')
4899         {
4900           *s = '0';
4901           goto roun;
4902         }
4903     }
4904  doexp:
4905   /* Strip trailing zeros, but leave at least one.  */
4906   while (ss[-1] == '0' && ss[-2] != '.')
4907     --ss;
4908   sprintf (ss, "e%d", expon);
4909  bxit:
4910   rndprc = rndsav;
4911   /* copy out the working string */
4912   s = string;
4913   ss = wstring;
4914   while (*ss == ' ')            /* strip possible leading space */
4915     ++ss;
4916   while ((*s++ = *ss++) != '\0')
4917     ;
4918 }
4919
4920
4921 /* Convert ASCII string to floating point.
4922
4923    Numeric input is a free format decimal number of any length, with
4924    or without decimal point.  Entering E after the number followed by an
4925    integer number causes the second number to be interpreted as a power of
4926    10 to be multiplied by the first number (i.e., "scientific" notation).  */
4927
4928 /* Convert ASCII string S to single precision float value Y.  */
4929
4930 static void
4931 asctoe24 (s, y)
4932      const char *s;
4933      UEMUSHORT *y;
4934 {
4935   asctoeg (s, y, 24);
4936 }
4937
4938
4939 /* Convert ASCII string S to double precision value Y.  */
4940
4941 static void
4942 asctoe53 (s, y)
4943      const char *s;
4944      UEMUSHORT *y;
4945 {
4946 #if defined(DEC) || defined(IBM)
4947   asctoeg (s, y, 56);
4948 #else
4949 #if defined(C4X)
4950   asctoeg (s, y, 32);
4951 #else
4952   asctoeg (s, y, 53);
4953 #endif
4954 #endif
4955 }
4956
4957
4958 /* Convert ASCII string S to double extended value Y.  */
4959
4960 static void
4961 asctoe64 (s, y)
4962      const char *s;
4963      UEMUSHORT *y;
4964 {
4965   asctoeg (s, y, 64);
4966 }
4967
4968 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
4969 /* Convert ASCII string S to 128-bit long double Y.  */
4970
4971 static void
4972 asctoe113 (s, y)
4973      const char *s;
4974      UEMUSHORT *y;
4975 {
4976   asctoeg (s, y, 113);
4977 }
4978 #endif
4979
4980 /* Convert ASCII string S to e type Y.  */
4981
4982 static void
4983 asctoe (s, y)
4984      const char *s;
4985      UEMUSHORT *y;
4986 {
4987   asctoeg (s, y, NBITS);
4988 }
4989
4990 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4991    of OPREC bits.  BASE is 16 for C99 hexadecimal floating constants.  */
4992
4993 static void
4994 asctoeg (ss, y, oprec)
4995      const char *ss;
4996      UEMUSHORT *y;
4997      int oprec;
4998 {
4999   UEMUSHORT yy[NI], xt[NI], tt[NI];
5000   int esign, decflg, sgnflg, nexp, exp, prec, lost;
5001   int i, k, trail, c, rndsav;
5002   EMULONG lexp;
5003   UEMUSHORT nsign;
5004   char *sp, *s, *lstr;
5005   int base = 10;
5006
5007   /* Copy the input string.  */
5008   lstr = (char *) alloca (strlen (ss) + 1);
5009
5010   while (*ss == ' ')            /* skip leading spaces */
5011     ++ss;
5012
5013   sp = lstr;
5014   while ((*sp++ = *ss++) != '\0')
5015     ;
5016   s = lstr;
5017
5018   if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
5019     {
5020       base = 16;
5021       s += 2;
5022     }
5023
5024   rndsav = rndprc;
5025   rndprc = NBITS;               /* Set to full precision */
5026   lost = 0;
5027   nsign = 0;
5028   decflg = 0;
5029   sgnflg = 0;
5030   nexp = 0;
5031   exp = 0;
5032   prec = 0;
5033   ecleaz (yy);
5034   trail = 0;
5035
5036  nxtcom:
5037   k = hex_value (*s);
5038   if ((k >= 0) && (k < base))
5039     {
5040       /* Ignore leading zeros */
5041       if ((prec == 0) && (decflg == 0) && (k == 0))
5042         goto donchr;
5043       /* Identify and strip trailing zeros after the decimal point.  */
5044       if ((trail == 0) && (decflg != 0))
5045         {
5046           sp = s;
5047           while (ISDIGIT (*sp) || (base == 16 && ISXDIGIT (*sp)))
5048             ++sp;
5049           /* Check for syntax error */
5050           c = *sp & CHARMASK;
5051           if ((base != 10 || ((c != 'e') && (c != 'E')))
5052               && (base != 16 || ((c != 'p') && (c != 'P')))
5053               && (c != '\0')
5054               && (c != '\n') && (c != '\r') && (c != ' ')
5055               && (c != ','))
5056             goto unexpected_char_error;
5057           --sp;
5058           while (*sp == '0')
5059             *sp-- = 'z';
5060           trail = 1;
5061           if (*s == 'z')
5062             goto donchr;
5063         }
5064
5065       /* If enough digits were given to more than fill up the yy register,
5066          continuing until overflow into the high guard word yy[2]
5067          guarantees that there will be a roundoff bit at the top
5068          of the low guard word after normalization.  */
5069
5070       if (yy[2] == 0)
5071         {
5072           if (base == 16)
5073             {
5074               if (decflg)
5075                 nexp += 4;      /* count digits after decimal point */
5076
5077               eshup1 (yy);      /* multiply current number by 16 */
5078               eshup1 (yy);
5079               eshup1 (yy);
5080               eshup1 (yy);
5081             }
5082           else
5083             {
5084               if (decflg)
5085                 nexp += 1;              /* count digits after decimal point */
5086
5087               eshup1 (yy);              /* multiply current number by 10 */
5088               emovz (yy, xt);
5089               eshup1 (xt);
5090               eshup1 (xt);
5091               eaddm (xt, yy);
5092             }
5093           /* Insert the current digit.  */
5094           ecleaz (xt);
5095           xt[NI - 2] = (UEMUSHORT) k;
5096           eaddm (xt, yy);
5097         }
5098       else
5099         {
5100           /* Mark any lost non-zero digit.  */
5101           lost |= k;
5102           /* Count lost digits before the decimal point.  */
5103           if (decflg == 0)
5104             {
5105               if (base == 10)
5106                 nexp -= 1;
5107               else
5108                 nexp -= 4;
5109             }
5110         }
5111       prec += 1;
5112       goto donchr;
5113     }
5114
5115   switch (*s)
5116     {
5117     case 'z':
5118       break;
5119     case 'E':
5120     case 'e':
5121     case 'P':
5122     case 'p':
5123       goto expnt;
5124     case '.':                   /* decimal point */
5125       if (decflg)
5126         goto unexpected_char_error;
5127       ++decflg;
5128       break;
5129     case '-':
5130       nsign = 0xffff;
5131       if (sgnflg)
5132         goto unexpected_char_error;
5133       ++sgnflg;
5134       break;
5135     case '+':
5136       if (sgnflg)
5137         goto unexpected_char_error;
5138       ++sgnflg;
5139       break;
5140     case ',':
5141     case ' ':
5142     case '\0':
5143     case '\n':
5144     case '\r':
5145       goto daldone;
5146     case 'i':
5147     case 'I':
5148       goto infinite;
5149     default:
5150     unexpected_char_error:
5151 #ifdef NANS
5152       einan (yy);
5153 #else
5154       mtherr ("asctoe", DOMAIN);
5155       eclear (yy);
5156 #endif
5157       goto aexit;
5158     }
5159  donchr:
5160   ++s;
5161   goto nxtcom;
5162
5163   /* Exponent interpretation */
5164  expnt:
5165   /* 0.0eXXX is zero, regardless of XXX.  Check for the 0.0.  */
5166   for (k = 0; k < NI; k++)
5167     {
5168       if (yy[k] != 0)
5169         goto read_expnt;
5170     }
5171   goto aexit;
5172
5173 read_expnt:
5174   esign = 1;
5175   exp = 0;
5176   ++s;
5177   /* check for + or - */
5178   if (*s == '-')
5179     {
5180       esign = -1;
5181       ++s;
5182     }
5183   if (*s == '+')
5184     ++s;
5185   while (ISDIGIT (*s))
5186     {
5187       exp *= 10;
5188       exp += *s++ - '0';
5189       if (exp > 999999)
5190         break;
5191     }
5192   if (esign < 0)
5193     exp = -exp;
5194   if ((exp > MAXDECEXP) && (base == 10))
5195     {
5196  infinite:
5197       ecleaz (yy);
5198       yy[E] = 0x7fff;           /* infinity */
5199       goto aexit;
5200     }
5201   if ((exp < MINDECEXP) && (base == 10))
5202     {
5203  zero:
5204       ecleaz (yy);
5205       goto aexit;
5206     }
5207
5208  daldone:
5209   if (base == 16)
5210     {
5211       /* Base 16 hexadecimal floating constant.  */
5212       if ((k = enormlz (yy)) > NBITS)
5213         {
5214           ecleaz (yy);
5215           goto aexit;
5216         }
5217       /* Adjust the exponent.  NEXP is the number of hex digits,
5218          EXP is a power of 2.  */
5219       lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5220       if (lexp > 0x7fff)
5221         goto infinite;
5222       if (lexp < 0)
5223         goto zero;
5224       yy[E] = lexp;
5225       goto expdon;
5226     }
5227
5228   nexp = exp - nexp;
5229   /* Pad trailing zeros to minimize power of 10, per IEEE spec.  */
5230   while ((nexp > 0) && (yy[2] == 0))
5231     {
5232       emovz (yy, xt);
5233       eshup1 (xt);
5234       eshup1 (xt);
5235       eaddm (yy, xt);
5236       eshup1 (xt);
5237       if (xt[2] != 0)
5238         break;
5239       nexp -= 1;
5240       emovz (xt, yy);
5241     }
5242   if ((k = enormlz (yy)) > NBITS)
5243     {
5244       ecleaz (yy);
5245       goto aexit;
5246     }
5247   lexp = (EXONE - 1 + NBITS) - k;
5248   emdnorm (yy, lost, 0, lexp, 64);
5249   lost = 0;
5250
5251   /* Convert to external format:
5252
5253      Multiply by 10**nexp.  If precision is 64 bits,
5254      the maximum relative error incurred in forming 10**n
5255      for 0 <= n <= 324 is 8.2e-20, at 10**180.
5256      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5257      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
5258
5259   lexp = yy[E];
5260   if (nexp == 0)
5261     {
5262       k = 0;
5263       goto expdon;
5264     }
5265   esign = 1;
5266   if (nexp < 0)
5267     {
5268       nexp = -nexp;
5269       esign = -1;
5270       if (nexp > 4096)
5271         {
5272           /* Punt.  Can't handle this without 2 divides.  */
5273           emovi (etens[0], tt);
5274           lexp -= tt[E];
5275           k = edivm (tt, yy);
5276           lexp += EXONE;
5277           nexp -= 4096;
5278         }
5279     }
5280   emov (eone, xt);
5281   exp = 1;
5282   i = NTEN;
5283   do
5284     {
5285       if (exp & nexp)
5286         emul (etens[i], xt, xt);
5287       i--;
5288       exp = exp + exp;
5289     }
5290   while (exp <= MAXP);
5291
5292   emovi (xt, tt);
5293   if (esign < 0)
5294     {
5295       lexp -= tt[E];
5296       k = edivm (tt, yy);
5297       lexp += EXONE;
5298     }
5299   else
5300     {
5301       lexp += tt[E];
5302       k = emulm (tt, yy);
5303       lexp -= EXONE - 1;
5304     }
5305   lost = k;
5306
5307  expdon:
5308
5309   /* Round and convert directly to the destination type */
5310   if (oprec == 53)
5311     lexp -= EXONE - 0x3ff;
5312 #ifdef C4X
5313   else if (oprec == 24 || oprec == 32)
5314     lexp -= (EXONE - 0x7f);
5315 #else
5316 #ifdef IBM
5317   else if (oprec == 24 || oprec == 56)
5318     lexp -= EXONE - (0x41 << 2);
5319 #else
5320 #ifdef DEC
5321   else if (oprec == 24)
5322     lexp -= dec_f.adjustment;
5323   else if (oprec == 56)
5324     {
5325       if (TARGET_G_FLOAT)
5326         lexp -= dec_g.adjustment;
5327       else
5328         lexp -= dec_d.adjustment;
5329     }
5330 #else
5331 else if (oprec == 24)
5332     lexp -= EXONE - 0177;
5333 #endif /* DEC */
5334 #endif /* IBM */
5335 #endif /* C4X */
5336   rndprc = oprec;
5337   emdnorm (yy, lost, 0, lexp, 64);
5338
5339  aexit:
5340
5341   rndprc = rndsav;
5342   yy[0] = nsign;
5343   switch (oprec)
5344     {
5345 #ifdef DEC
5346     case 56:
5347       todec (yy, y);
5348       break;
5349 #endif
5350 #ifdef IBM
5351     case 56:
5352       toibm (yy, y, DFmode);
5353       break;
5354 #endif
5355 #ifdef C4X
5356     case 32:
5357       toc4x (yy, y, HFmode);
5358       break;
5359 #endif
5360
5361     case 53:
5362       toe53 (yy, y);
5363       break;
5364     case 24:
5365       toe24 (yy, y);
5366       break;
5367     case 64:
5368       toe64 (yy, y);
5369       break;
5370 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5371     case 113:
5372       toe113 (yy, y);
5373       break;
5374 #endif
5375     case NBITS:
5376       emovo (yy, y);
5377       break;
5378     }
5379 }
5380
5381
5382
5383 /* Return Y = largest integer not greater than X (truncated toward minus
5384    infinity).  */
5385
5386 static const UEMUSHORT bmask[] =
5387 {
5388   0xffff,
5389   0xfffe,
5390   0xfffc,
5391   0xfff8,
5392   0xfff0,
5393   0xffe0,
5394   0xffc0,
5395   0xff80,
5396   0xff00,
5397   0xfe00,
5398   0xfc00,
5399   0xf800,
5400   0xf000,
5401   0xe000,
5402   0xc000,
5403   0x8000,
5404   0x0000,
5405 };
5406
5407 static void
5408 efloor (x, y)
5409      const UEMUSHORT x[];
5410      UEMUSHORT y[];
5411 {
5412   UEMUSHORT *p;
5413   int e, expon, i;
5414   UEMUSHORT f[NE];
5415
5416   emov (x, f);                  /* leave in external format */
5417   expon = (int) f[NE - 1];
5418   e = (expon & 0x7fff) - (EXONE - 1);
5419   if (e <= 0)
5420     {
5421       eclear (y);
5422       goto isitneg;
5423     }
5424   /* number of bits to clear out */
5425   e = NBITS - e;
5426   emov (f, y);
5427   if (e <= 0)
5428     return;
5429
5430   p = &y[0];
5431   while (e >= 16)
5432     {
5433       *p++ = 0;
5434       e -= 16;
5435     }
5436   /* clear the remaining bits */
5437   *p &= bmask[e];
5438   /* truncate negatives toward minus infinity */
5439  isitneg:
5440
5441   if ((UEMUSHORT) expon & (UEMUSHORT) 0x8000)
5442     {
5443       for (i = 0; i < NE - 1; i++)
5444         {
5445           if (f[i] != y[i])
5446             {
5447               esub (eone, y, y);
5448               break;
5449             }
5450         }
5451     }
5452 }
5453
5454
5455 #if 0
5456 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5457    For example, 1.1 = 0.55 * 2^1.  */
5458
5459 static void
5460 efrexp (x, exp, s)
5461      const UEMUSHORT x[];
5462      int *exp;
5463      UEMUSHORT s[];
5464 {
5465   UEMUSHORT xi[NI];
5466   EMULONG li;
5467
5468   emovi (x, xi);
5469   /*  Handle denormalized numbers properly using long integer exponent.  */
5470   li = (EMULONG) ((EMUSHORT) xi[1]);
5471
5472   if (li == 0)
5473     {
5474       li -= enormlz (xi);
5475     }
5476   xi[1] = 0x3ffe;
5477   emovo (xi, s);
5478   *exp = (int) (li - 0x3ffe);
5479 }
5480 #endif
5481
5482 /* Return e type Y = X * 2^PWR2.  */
5483
5484 static void
5485 eldexp (x, pwr2, y)
5486      const UEMUSHORT x[];
5487      int pwr2;
5488      UEMUSHORT y[];
5489 {
5490   UEMUSHORT xi[NI];
5491   EMULONG li;
5492   int i;
5493
5494   emovi (x, xi);
5495   li = xi[1];
5496   li += pwr2;
5497   i = 0;
5498   emdnorm (xi, i, i, li, !ROUND_TOWARDS_ZERO);
5499   emovo (xi, y);
5500 }
5501
5502
5503 #if 0
5504 /* C = remainder after dividing B by A, all e type values.
5505    Least significant integer quotient bits left in EQUOT.  */
5506
5507 static void
5508 eremain (a, b, c)
5509      const UEMUSHORT a[], b[];
5510      UEMUSHORT c[];
5511 {
5512   UEMUSHORT den[NI], num[NI];
5513
5514 #ifdef NANS
5515   if (eisinf (b)
5516       || (ecmp (a, ezero) == 0)
5517       || eisnan (a)
5518       || eisnan (b))
5519     {
5520       enan (c, 0);
5521       return;
5522     }
5523 #endif
5524   if (ecmp (a, ezero) == 0)
5525     {
5526       mtherr ("eremain", SING);
5527       eclear (c);
5528       return;
5529     }
5530   emovi (a, den);
5531   emovi (b, num);
5532   eiremain (den, num);
5533   /* Sign of remainder = sign of quotient */
5534   if (a[0] == b[0])
5535     num[0] = 0;
5536   else
5537     num[0] = 0xffff;
5538   emovo (num, c);
5539 }
5540 #endif
5541
5542 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5543     remainder in NUM.  */
5544
5545 static void
5546 eiremain (den, num)
5547      UEMUSHORT den[], num[];
5548 {
5549   EMULONG ld, ln;
5550   UEMUSHORT j;
5551
5552   ld = den[E];
5553   ld -= enormlz (den);
5554   ln = num[E];
5555   ln -= enormlz (num);
5556   ecleaz (equot);
5557   while (ln >= ld)
5558     {
5559       if (ecmpm (den, num) <= 0)
5560         {
5561           esubm (den, num);
5562           j = 1;
5563         }
5564       else
5565           j = 0;
5566       eshup1 (equot);
5567       equot[NI - 1] |= j;
5568       eshup1 (num);
5569       ln -= 1;
5570     }
5571   emdnorm (num, 0, 0, ln, 0);
5572 }
5573
5574 /* Report an error condition CODE encountered in function NAME.
5575
5576     Mnemonic        Value          Significance
5577
5578      DOMAIN            1       argument domain error
5579      SING              2       function singularity
5580      OVERFLOW          3       overflow range error
5581      UNDERFLOW         4       underflow range error
5582      TLOSS             5       total loss of precision
5583      PLOSS             6       partial loss of precision
5584      INVALID           7       NaN - producing operation
5585      EDOM             33       Unix domain error code
5586      ERANGE           34       Unix range error code
5587
5588    The order of appearance of the following messages is bound to the
5589    error codes defined above.  */
5590
5591 int merror = 0;
5592 extern int merror;
5593
5594 static void
5595 mtherr (name, code)
5596      const char *name;
5597      int code;
5598 {
5599   /* The string passed by the calling program is supposed to be the
5600      name of the function in which the error occurred.
5601      The code argument selects which error message string will be printed.  */
5602
5603   if (strcmp (name, "esub") == 0)
5604     name = "subtraction";
5605   else if (strcmp (name, "ediv") == 0)
5606     name = "division";
5607   else if (strcmp (name, "emul") == 0)
5608     name = "multiplication";
5609   else if (strcmp (name, "enormlz") == 0)
5610     name = "normalization";
5611   else if (strcmp (name, "etoasc") == 0)
5612     name = "conversion to text";
5613   else if (strcmp (name, "asctoe") == 0)
5614     name = "parsing";
5615   else if (strcmp (name, "eremain") == 0)
5616     name = "modulus";
5617   else if (strcmp (name, "esqrt") == 0)
5618     name = "square root";
5619   if (extra_warnings)
5620     {
5621       switch (code)
5622         {
5623         case DOMAIN:    warning ("%s: argument domain error"    , name); break;
5624         case SING:      warning ("%s: function singularity"     , name); break;
5625         case OVERFLOW:  warning ("%s: overflow range error"     , name); break;
5626         case UNDERFLOW: warning ("%s: underflow range error"    , name); break;
5627         case TLOSS:     warning ("%s: total loss of precision"  , name); break;
5628         case PLOSS:     warning ("%s: partial loss of precision", name); break;
5629         case INVALID:   warning ("%s: NaN - producing operation", name); break;
5630         default:        abort ();
5631         }
5632     }
5633
5634   /* Set global error message word */
5635   merror = code + 1;
5636 }
5637
5638 #ifdef DEC
5639 /* Convert DEC double precision D to e type E.  */
5640
5641 static void
5642 dectoe (d, e)
5643      const UEMUSHORT *d;
5644      UEMUSHORT *e;
5645 {
5646   if (TARGET_G_FLOAT)
5647     ieeetoe (d, e, &dec_g);
5648   else
5649     ieeetoe (d, e, &dec_d);
5650 }
5651
5652 /* Convert e type X to DEC double precision D.  */
5653
5654 static void
5655 etodec (x, d)
5656      const UEMUSHORT *x;
5657      UEMUSHORT *d;
5658 {
5659   UEMUSHORT xi[NI];
5660   EMULONG exp;
5661   int rndsav;
5662   const struct ieee_format *fmt;
5663
5664   if (TARGET_G_FLOAT)
5665     fmt = &dec_g;
5666   else
5667     fmt = &dec_d;
5668
5669   emovi (x, xi);
5670   /* Adjust exponent for offsets.  */
5671   exp = (EMULONG) xi[E] - fmt->adjustment;
5672   /* Round off to nearest or even.  */
5673   rndsav = rndprc;
5674   rndprc = fmt->precision;
5675   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5676   rndprc = rndsav;
5677   todec (xi, d);
5678 }
5679
5680 /* Convert exploded e-type X, that has already been rounded to
5681    56-bit precision, to DEC format double Y.  */
5682
5683 static void
5684 todec (x, y)
5685      UEMUSHORT *x, *y;
5686 {
5687   if (TARGET_G_FLOAT)
5688     toieee (x, y, &dec_g);
5689   else
5690     toieee (x, y, &dec_d);
5691 }
5692 #endif /* DEC */
5693
5694 #ifdef IBM
5695 /* Convert IBM single/double precision to e type.  */
5696
5697 static void
5698 ibmtoe (d, e, mode)
5699      const UEMUSHORT *d;
5700      UEMUSHORT *e;
5701      enum machine_mode mode;
5702 {
5703   UEMUSHORT y[NI];
5704   UEMUSHORT r, *p;
5705
5706   ecleaz (y);                   /* start with a zero */
5707   p = y;                        /* point to our number */
5708   r = *d;                       /* get IBM exponent word */
5709   if (*d & (unsigned int) 0x8000)
5710     *p = 0xffff;                /* fill in our sign */
5711   ++p;                          /* bump pointer to our exponent word */
5712   r &= 0x7f00;                  /* strip the sign bit */
5713   r >>= 6;                      /* shift exponent word down 6 bits */
5714                                 /* in fact shift by 8 right and 2 left */
5715   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5716                                 /* add our e type exponent offset */
5717   *p++ = r;                     /* to form our exponent */
5718
5719   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5720                                 /* strip off the IBM exponent and sign bits */
5721   if (mode != SFmode)           /* there are only 2 words in SFmode */
5722     {
5723       *p++ = *d++;              /* fill in the rest of our mantissa */
5724       *p++ = *d++;
5725     }
5726   *p = *d;
5727
5728   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5729     y[0] = y[E] = 0;
5730   else
5731     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5732                               /* handle change in RADIX */
5733   emovo (y, e);
5734 }
5735
5736
5737
5738 /* Convert e type to IBM single/double precision.  */
5739
5740 static void
5741 etoibm (x, d, mode)
5742      const UEMUSHORT *x;
5743      UEMUSHORT *d;
5744      enum machine_mode mode;
5745 {
5746   UEMUSHORT xi[NI];
5747   EMULONG exp;
5748   int rndsav;
5749
5750   emovi (x, xi);
5751   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5752                                                         /* round off to nearest or even */
5753   rndsav = rndprc;
5754   rndprc = 56;
5755   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5756   rndprc = rndsav;
5757   toibm (xi, d, mode);
5758 }
5759
5760 static void
5761 toibm (x, y, mode)
5762      UEMUSHORT *x, *y;
5763      enum machine_mode mode;
5764 {
5765   UEMUSHORT i;
5766   UEMUSHORT *p;
5767   int r;
5768
5769   p = x;
5770   *y = 0;
5771   if (*p++)
5772     *y = 0x8000;
5773   i = *p++;
5774   if (i == 0)
5775     {
5776       *y++ = 0;
5777       *y++ = 0;
5778       if (mode != SFmode)
5779         {
5780           *y++ = 0;
5781           *y++ = 0;
5782         }
5783       return;
5784     }
5785   r = i & 0x3;
5786   i >>= 2;
5787   if (i > 0x7f)
5788     {
5789       *y++ |= 0x7fff;
5790       *y++ = 0xffff;
5791       if (mode != SFmode)
5792         {
5793           *y++ = 0xffff;
5794           *y++ = 0xffff;
5795         }
5796 #ifdef ERANGE
5797       errno = ERANGE;
5798 #endif
5799       return;
5800     }
5801   i &= 0x7f;
5802   *y |= (i << 8);
5803   eshift (x, r + 5);
5804   *y++ |= x[M];
5805   *y++ = x[M + 1];
5806   if (mode != SFmode)
5807     {
5808       *y++ = x[M + 2];
5809       *y++ = x[M + 3];
5810     }
5811 }
5812 #endif /* IBM */
5813
5814
5815 #ifdef C4X
5816 /* Convert C4X single/double precision to e type.  */
5817
5818 static void
5819 c4xtoe (d, e, mode)
5820      const UEMUSHORT *d;
5821      UEMUSHORT *e;
5822      enum machine_mode mode;
5823 {
5824   UEMUSHORT y[NI];
5825   UEMUSHORT dn[4];
5826   int r;
5827   int isnegative;
5828   int size;
5829   int i;
5830   int carry;
5831
5832   dn[0] = d[0];
5833   dn[1] = d[1];
5834   if (mode != QFmode)
5835     {
5836       dn[2] = d[3] << 8;
5837       dn[3] = 0;
5838     }
5839
5840   /* Short-circuit the zero case.  */
5841   if ((dn[0] == 0x8000)
5842       && (dn[1] == 0x0000)
5843       && ((mode == QFmode) || ((dn[2] == 0x0000) && (dn[3] == 0x0000))))
5844     {
5845       e[0] = 0;
5846       e[1] = 0;
5847       e[2] = 0;
5848       e[3] = 0;
5849       e[4] = 0;
5850       e[5] = 0;
5851       return;
5852     }
5853
5854   ecleaz (y);                   /* start with a zero */
5855   r = dn[0];                    /* get sign/exponent part */
5856   if (r & (unsigned int) 0x0080)
5857     {
5858       y[0] = 0xffff;            /* fill in our sign */
5859       isnegative = TRUE;
5860     }
5861   else
5862     isnegative = FALSE;
5863
5864   r >>= 8;                      /* Shift exponent word down 8 bits.  */
5865   if (r & 0x80)                 /* Make the exponent negative if it is.  */
5866     r = r | (~0 & ~0xff);
5867
5868   if (isnegative)
5869     {
5870       /* Now do the high order mantissa.  We don't "or" on the high bit
5871          because it is 2 (not 1) and is handled a little differently
5872          below.  */
5873       y[M] = dn[0] & 0x7f;
5874
5875       y[M+1] = dn[1];
5876       if (mode != QFmode)       /* There are only 2 words in QFmode.  */
5877         {
5878           y[M+2] = dn[2];       /* Fill in the rest of our mantissa.  */
5879           y[M+3] = dn[3];
5880           size = 4;
5881         }
5882       else
5883         size = 2;
5884       eshift (y, -8);
5885
5886       /* Now do the two's complement on the data.  */
5887
5888       carry = 1;        /* Initially add 1 for the two's complement.  */
5889       for (i=size + M; i > M; i--)
5890         {
5891           if (carry && (y[i] == 0x0000))
5892             /* We overflowed into the next word, carry is the same.  */
5893             y[i] = carry ? 0x0000 : 0xffff;
5894           else
5895             {
5896               /* No overflow, just invert and add carry.  */
5897               y[i] = ((~y[i]) + carry) & 0xffff;
5898               carry = 0;
5899             }
5900         }
5901
5902       if (carry)
5903         {
5904           eshift (y, -1);
5905           y[M+1] |= 0x8000;
5906           r++;
5907         }
5908        y[1] = r + EXONE;
5909     }
5910   else
5911     {
5912       /* Add our e type exponent offset to form our exponent.  */
5913       r += EXONE;
5914       y[1] = r;
5915
5916      /* Now do the high order mantissa strip off the exponent and sign
5917         bits and add the high 1 bit.  */
5918      y[M] = (dn[0] & 0x7f) | 0x80;
5919
5920      y[M+1] = dn[1];
5921      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
5922        {
5923          y[M+2] = dn[2];        /* Fill in the rest of our mantissa.  */
5924          y[M+3] = dn[3];
5925        }
5926      eshift (y, -8);
5927     }
5928
5929   emovo (y, e);
5930 }
5931
5932
5933 /* Convert e type to C4X single/double precision.  */
5934
5935 static void
5936 etoc4x (x, d, mode)
5937      const UEMUSHORT *x;
5938      UEMUSHORT *d;
5939      enum machine_mode mode;
5940 {
5941   UEMUSHORT xi[NI];
5942   EMULONG exp;
5943   int rndsav;
5944
5945   emovi (x, xi);
5946
5947   /* Adjust exponent for offsets.  */
5948   exp = (EMULONG) xi[E] - (EXONE - 0x7f);
5949
5950   /* Round off to nearest or even.  */
5951   rndsav = rndprc;
5952   rndprc = mode == QFmode ? 24 : 32;
5953   emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5954   rndprc = rndsav;
5955   toc4x (xi, d, mode);
5956 }
5957
5958 static void
5959 toc4x (x, y, mode)
5960      UEMUSHORT *x, *y;
5961      enum machine_mode mode;
5962 {
5963   int i;
5964   int v;
5965   int carry;
5966
5967   /* Short-circuit the zero case */
5968   if ((x[0] == 0)       /* Zero exponent and sign */
5969       && (x[1] == 0)
5970       && (x[M] == 0)    /* The rest is for zero mantissa */
5971       && (x[M+1] == 0)
5972       /* Only check for double if necessary */
5973       && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
5974     {
5975       /* We have a zero.  Put it into the output and return.  */
5976       *y++ = 0x8000;
5977       *y++ = 0x0000;
5978       if (mode != QFmode)
5979         {
5980           *y++ = 0x0000;
5981           *y++ = 0x0000;
5982         }
5983       return;
5984     }
5985
5986   *y = 0;
5987
5988   /* Negative number require a two's complement conversion of the
5989      mantissa.  */
5990   if (x[0])
5991     {
5992       *y = 0x0080;
5993
5994       i = ((int) x[1]) - 0x7f;
5995
5996       /* Now add 1 to the inverted data to do the two's complement.  */
5997       if (mode != QFmode)
5998         v = 4 + M;
5999       else
6000         v = 2 + M;
6001       carry = 1;
6002       while (v > M)
6003         {
6004           if (x[v] == 0x0000)
6005             x[v] = carry ? 0x0000 : 0xffff;
6006           else
6007             {
6008               x[v] = ((~x[v]) + carry) & 0xffff;
6009               carry = 0;
6010             }
6011           v--;
6012         }
6013
6014       /* The following is a special case.  The C4X negative float requires
6015          a zero in the high bit (because the format is (2 - x) x 2^m), so
6016          if a one is in that bit, we have to shift left one to get rid
6017          of it.  This only occurs if the number is -1 x 2^m.  */
6018       if (x[M+1] & 0x8000)
6019         {
6020           /* This is the case of -1 x 2^m, we have to rid ourselves of the
6021              high sign bit and shift the exponent.  */
6022           eshift (x, 1);
6023           i--;
6024         }
6025     }
6026   else
6027     i = ((int) x[1]) - 0x7f;
6028
6029   if ((i < -128) || (i > 127))
6030     {
6031       y[0] |= 0xff7f;
6032       y[1] = 0xffff;
6033       if (mode != QFmode)
6034         {
6035           y[2] = 0xffff;
6036           y[3] = 0xffff;
6037           y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6038           y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6039         }
6040 #ifdef ERANGE
6041       errno = ERANGE;
6042 #endif
6043       return;
6044     }
6045
6046   y[0] |= ((i & 0xff) << 8);
6047
6048   eshift (x, 8);
6049
6050   y[0] |= x[M] & 0x7f;
6051   y[1] = x[M + 1];
6052   if (mode != QFmode)
6053     {
6054       y[2] = x[M + 2];
6055       y[3] = x[M + 3];
6056       y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6057       y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6058     }
6059 }
6060 #endif /* C4X */
6061
6062 /* Output a binary NaN bit pattern in the target machine's format.  */
6063
6064 /* If special NaN bit patterns are required, define them in tm.h
6065    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
6066    patterns here.  */
6067 #ifdef TFMODE_NAN
6068 TFMODE_NAN;
6069 #else
6070 #if defined (IEEE) && (INTEL_EXTENDED_IEEE_FORMAT == 0)
6071 static const UEMUSHORT TFbignan[8] =
6072  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6073 static const UEMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6074 #endif
6075 #endif
6076
6077 #ifdef XFMODE_NAN
6078 XFMODE_NAN;
6079 #else
6080 #ifdef IEEE
6081 static const UEMUSHORT XFbignan[6] =
6082  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6083 static const UEMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6084 #endif
6085 #endif
6086
6087 #ifdef DFMODE_NAN
6088 DFMODE_NAN;
6089 #else
6090 #ifdef IEEE
6091 static const UEMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6092 static const UEMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
6093 #endif
6094 #endif
6095
6096 #ifdef SFMODE_NAN
6097 SFMODE_NAN;
6098 #else
6099 #ifdef IEEE
6100 static const UEMUSHORT SFbignan[2] = {0x7fff, 0xffff};
6101 static const UEMUSHORT SFlittlenan[2] = {0, 0xffc0};
6102 #endif
6103 #endif
6104
6105
6106 #ifdef NANS
6107 static void
6108 make_nan (nan, sign, mode)
6109      UEMUSHORT *nan;
6110      int sign;
6111      enum machine_mode mode;
6112 {
6113   int n;
6114   const UEMUSHORT *p;
6115   int size;
6116
6117   size = GET_MODE_BITSIZE (mode);
6118   if (LARGEST_EXPONENT_IS_NORMAL (size))
6119     {
6120       warning ("%d-bit floats cannot hold NaNs", size);
6121       saturate (nan, sign, size, 0);
6122       return;
6123     }
6124   switch (mode)
6125     {
6126 /* Possibly the `reserved operand' patterns on a VAX can be
6127    used like NaN's, but probably not in the same way as IEEE.  */
6128 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6129     case TFmode:
6130 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6131       n = 8;
6132       if (REAL_WORDS_BIG_ENDIAN)
6133         p = TFbignan;
6134       else
6135         p = TFlittlenan;
6136       break;
6137 #endif
6138       /* FALLTHRU */
6139
6140     case XFmode:
6141       n = 6;
6142       if (REAL_WORDS_BIG_ENDIAN)
6143         p = XFbignan;
6144       else
6145         p = XFlittlenan;
6146       break;
6147
6148     case DFmode:
6149       n = 4;
6150       if (REAL_WORDS_BIG_ENDIAN)
6151         p = DFbignan;
6152       else
6153         p = DFlittlenan;
6154       break;
6155
6156     case SFmode:
6157     case HFmode:
6158       n = 2;
6159       if (REAL_WORDS_BIG_ENDIAN)
6160         p = SFbignan;
6161       else
6162         p = SFlittlenan;
6163       break;
6164 #endif
6165
6166     default:
6167       abort ();
6168     }
6169   if (REAL_WORDS_BIG_ENDIAN)
6170     *nan++ = (sign << 15) | (*p++ & 0x7fff);
6171   while (--n != 0)
6172     *nan++ = *p++;
6173   if (! REAL_WORDS_BIG_ENDIAN)
6174     *nan = (sign << 15) | (*p & 0x7fff);
6175 }
6176 #endif /* NANS */
6177
6178
6179 /* Create a saturation value for a SIZE-bit float, assuming that
6180    LARGEST_EXPONENT_IS_NORMAL (SIZE).
6181
6182    If SIGN is true, fill X with the most negative value, otherwise fill
6183    it with the most positive value.  WARN is true if the function should
6184    warn about overflow.  */
6185
6186 static void
6187 saturate (x, sign, size, warn)
6188      UEMUSHORT *x;
6189      int sign, size, warn;
6190 {
6191   int i;
6192
6193   if (warn && extra_warnings)
6194     warning ("value exceeds the range of a %d-bit float", size);
6195
6196   /* Create the most negative value.  */
6197   for (i = 0; i < size / EMUSHORT_SIZE; i++)
6198     x[i] = 0xffff;
6199
6200   /* Make it positive, if necessary.  */
6201   if (!sign)
6202     x[REAL_WORDS_BIG_ENDIAN? 0 : i - 1] = 0x7fff;
6203 }
6204
6205
6206 /* This is the inverse of the function `etarsingle' invoked by
6207    REAL_VALUE_TO_TARGET_SINGLE.  */
6208
6209 REAL_VALUE_TYPE
6210 ereal_unto_float (f)
6211      long f;
6212 {
6213   REAL_VALUE_TYPE r;
6214   UEMUSHORT s[2];
6215   UEMUSHORT e[NE];
6216
6217   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6218    This is the inverse operation to what the function `endian' does.  */
6219   if (REAL_WORDS_BIG_ENDIAN)
6220     {
6221       s[0] = (UEMUSHORT) (f >> 16);
6222       s[1] = (UEMUSHORT) f;
6223     }
6224   else
6225     {
6226       s[0] = (UEMUSHORT) f;
6227       s[1] = (UEMUSHORT) (f >> 16);
6228     }
6229   /* Convert and promote the target float to E-type.  */
6230   e24toe (s, e);
6231   /* Output E-type to REAL_VALUE_TYPE.  */
6232   PUT_REAL (e, &r);
6233   return r;
6234 }
6235
6236
6237 /* This is the inverse of the function `etardouble' invoked by
6238    REAL_VALUE_TO_TARGET_DOUBLE.  */
6239
6240 REAL_VALUE_TYPE
6241 ereal_unto_double (d)
6242      long d[];
6243 {
6244   REAL_VALUE_TYPE r;
6245   UEMUSHORT s[4];
6246   UEMUSHORT e[NE];
6247
6248   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6249   if (REAL_WORDS_BIG_ENDIAN)
6250     {
6251       s[0] = (UEMUSHORT) (d[0] >> 16);
6252       s[1] = (UEMUSHORT) d[0];
6253       s[2] = (UEMUSHORT) (d[1] >> 16);
6254       s[3] = (UEMUSHORT) d[1];
6255     }
6256   else
6257     {
6258       /* Target float words are little-endian.  */
6259       s[0] = (UEMUSHORT) d[0];
6260       s[1] = (UEMUSHORT) (d[0] >> 16);
6261       s[2] = (UEMUSHORT) d[1];
6262       s[3] = (UEMUSHORT) (d[1] >> 16);
6263     }
6264   /* Convert target double to E-type.  */
6265   e53toe (s, e);
6266   /* Output E-type to REAL_VALUE_TYPE.  */
6267   PUT_REAL (e, &r);
6268   return r;
6269 }
6270
6271
6272 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6273    This is somewhat like ereal_unto_float, but the input types
6274    for these are different.  */
6275
6276 REAL_VALUE_TYPE
6277 ereal_from_float (f)
6278      HOST_WIDE_INT f;
6279 {
6280   REAL_VALUE_TYPE r;
6281   UEMUSHORT s[2];
6282   UEMUSHORT e[NE];
6283
6284   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6285    This is the inverse operation to what the function `endian' does.  */
6286   if (REAL_WORDS_BIG_ENDIAN)
6287     {
6288       s[0] = (UEMUSHORT) (f >> 16);
6289       s[1] = (UEMUSHORT) f;
6290     }
6291   else
6292     {
6293       s[0] = (UEMUSHORT) f;
6294       s[1] = (UEMUSHORT) (f >> 16);
6295     }
6296   /* Convert and promote the target float to E-type.  */
6297   e24toe (s, e);
6298   /* Output E-type to REAL_VALUE_TYPE.  */
6299   PUT_REAL (e, &r);
6300   return r;
6301 }
6302
6303
6304 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6305    This is somewhat like ereal_unto_double, but the input types
6306    for these are different.
6307
6308    The DFmode is stored as an array of HOST_WIDE_INT in the target's
6309    data format, with no holes in the bit packing.  The first element
6310    of the input array holds the bits that would come first in the
6311    target computer's memory.  */
6312
6313 REAL_VALUE_TYPE
6314 ereal_from_double (d)
6315      HOST_WIDE_INT d[];
6316 {
6317   REAL_VALUE_TYPE r;
6318   UEMUSHORT s[4];
6319   UEMUSHORT e[NE];
6320
6321   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6322   if (REAL_WORDS_BIG_ENDIAN)
6323     {
6324 #if HOST_BITS_PER_WIDE_INT == 32
6325       s[0] = (UEMUSHORT) (d[0] >> 16);
6326       s[1] = (UEMUSHORT) d[0];
6327       s[2] = (UEMUSHORT) (d[1] >> 16);
6328       s[3] = (UEMUSHORT) d[1];
6329 #else
6330       /* In this case the entire target double is contained in the
6331          first array element.  The second element of the input is
6332          ignored.  */
6333       s[0] = (UEMUSHORT) (d[0] >> 48);
6334       s[1] = (UEMUSHORT) (d[0] >> 32);
6335       s[2] = (UEMUSHORT) (d[0] >> 16);
6336       s[3] = (UEMUSHORT) d[0];
6337 #endif
6338     }
6339   else
6340     {
6341       /* Target float words are little-endian.  */
6342       s[0] = (UEMUSHORT) d[0];
6343       s[1] = (UEMUSHORT) (d[0] >> 16);
6344 #if HOST_BITS_PER_WIDE_INT == 32
6345       s[2] = (UEMUSHORT) d[1];
6346       s[3] = (UEMUSHORT) (d[1] >> 16);
6347 #else
6348       s[2] = (UEMUSHORT) (d[0] >> 32);
6349       s[3] = (UEMUSHORT) (d[0] >> 48);
6350 #endif
6351     }
6352   /* Convert target double to E-type.  */
6353   e53toe (s, e);
6354   /* Output E-type to REAL_VALUE_TYPE.  */
6355   PUT_REAL (e, &r);
6356   return r;
6357 }
6358
6359
6360 #if 0
6361 /* Convert target computer unsigned 64-bit integer to e-type.
6362    The endian-ness of DImode follows the convention for integers,
6363    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
6364
6365 static void
6366 uditoe (di, e)
6367      const UEMUSHORT *di;  /* Address of the 64-bit int.  */
6368      UEMUSHORT *e;
6369 {
6370   UEMUSHORT yi[NI];
6371   int k;
6372
6373   ecleaz (yi);
6374   if (WORDS_BIG_ENDIAN)
6375     {
6376       for (k = M; k < M + 4; k++)
6377         yi[k] = *di++;
6378     }
6379   else
6380     {
6381       for (k = M + 3; k >= M; k--)
6382         yi[k] = *di++;
6383     }
6384   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6385   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6386     ecleaz (yi);                /* it was zero */
6387   else
6388     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6389   emovo (yi, e);
6390 }
6391
6392 /* Convert target computer signed 64-bit integer to e-type.  */
6393
6394 static void
6395 ditoe (di, e)
6396      const UEMUSHORT *di;  /* Address of the 64-bit int.  */
6397      UEMUSHORT *e;
6398 {
6399   unsigned EMULONG acc;
6400   UEMUSHORT yi[NI];
6401   UEMUSHORT carry;
6402   int k, sign;
6403
6404   ecleaz (yi);
6405   if (WORDS_BIG_ENDIAN)
6406     {
6407       for (k = M; k < M + 4; k++)
6408         yi[k] = *di++;
6409     }
6410   else
6411     {
6412       for (k = M + 3; k >= M; k--)
6413         yi[k] = *di++;
6414     }
6415   /* Take absolute value */
6416   sign = 0;
6417   if (yi[M] & 0x8000)
6418     {
6419       sign = 1;
6420       carry = 0;
6421       for (k = M + 3; k >= M; k--)
6422         {
6423           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6424           yi[k] = acc;
6425           carry = 0;
6426           if (acc & 0x10000)
6427             carry = 1;
6428         }
6429     }
6430   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6431   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6432     ecleaz (yi);                /* it was zero */
6433   else
6434     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6435   emovo (yi, e);
6436   if (sign)
6437         eneg (e);
6438 }
6439
6440
6441 /* Convert e-type to unsigned 64-bit int.  */
6442
6443 static void
6444 etoudi (x, i)
6445      const UEMUSHORT *x;
6446      UEMUSHORT *i;
6447 {
6448   UEMUSHORT xi[NI];
6449   int j, k;
6450
6451   emovi (x, xi);
6452   if (xi[0])
6453     {
6454       xi[M] = 0;
6455       goto noshift;
6456     }
6457   k = (int) xi[E] - (EXONE - 1);
6458   if (k <= 0)
6459     {
6460       for (j = 0; j < 4; j++)
6461         *i++ = 0;
6462       return;
6463     }
6464   if (k > 64)
6465     {
6466       for (j = 0; j < 4; j++)
6467         *i++ = 0xffff;
6468       if (extra_warnings)
6469         warning ("overflow on truncation to integer");
6470       return;
6471     }
6472   if (k > 16)
6473     {
6474       /* Shift more than 16 bits: first shift up k-16 mod 16,
6475          then shift up by 16's.  */
6476       j = k - ((k >> 4) << 4);
6477       if (j == 0)
6478         j = 16;
6479       eshift (xi, j);
6480       if (WORDS_BIG_ENDIAN)
6481         *i++ = xi[M];
6482       else
6483         {
6484           i += 3;
6485           *i-- = xi[M];
6486         }
6487       k -= j;
6488       do
6489         {
6490           eshup6 (xi);
6491           if (WORDS_BIG_ENDIAN)
6492             *i++ = xi[M];
6493           else
6494             *i-- = xi[M];
6495         }
6496       while ((k -= 16) > 0);
6497     }
6498   else
6499     {
6500       /* shift not more than 16 bits */
6501       eshift (xi, k);
6502
6503 noshift:
6504
6505       if (WORDS_BIG_ENDIAN)
6506         {
6507           i += 3;
6508           *i-- = xi[M];
6509           *i-- = 0;
6510           *i-- = 0;
6511           *i = 0;
6512         }
6513       else
6514         {
6515           *i++ = xi[M];
6516           *i++ = 0;
6517           *i++ = 0;
6518           *i = 0;
6519         }
6520     }
6521 }
6522
6523
6524 /* Convert e-type to signed 64-bit int.  */
6525
6526 static void
6527 etodi (x, i)
6528      const UEMUSHORT *x;
6529      UEMUSHORT *i;
6530 {
6531   unsigned EMULONG acc;
6532   UEMUSHORT xi[NI];
6533   UEMUSHORT carry;
6534   UEMUSHORT *isave;
6535   int j, k;
6536
6537   emovi (x, xi);
6538   k = (int) xi[E] - (EXONE - 1);
6539   if (k <= 0)
6540     {
6541       for (j = 0; j < 4; j++)
6542         *i++ = 0;
6543       return;
6544     }
6545   if (k > 64)
6546     {
6547       for (j = 0; j < 4; j++)
6548         *i++ = 0xffff;
6549       if (extra_warnings)
6550         warning ("overflow on truncation to integer");
6551       return;
6552     }
6553   isave = i;
6554   if (k > 16)
6555     {
6556       /* Shift more than 16 bits: first shift up k-16 mod 16,
6557          then shift up by 16's.  */
6558       j = k - ((k >> 4) << 4);
6559       if (j == 0)
6560         j = 16;
6561       eshift (xi, j);
6562       if (WORDS_BIG_ENDIAN)
6563         *i++ = xi[M];
6564       else
6565         {
6566           i += 3;
6567           *i-- = xi[M];
6568         }
6569       k -= j;
6570       do
6571         {
6572           eshup6 (xi);
6573           if (WORDS_BIG_ENDIAN)
6574             *i++ = xi[M];
6575           else
6576             *i-- = xi[M];
6577         }
6578       while ((k -= 16) > 0);
6579     }
6580   else
6581     {
6582       /* shift not more than 16 bits */
6583       eshift (xi, k);
6584
6585       if (WORDS_BIG_ENDIAN)
6586         {
6587           i += 3;
6588           *i = xi[M];
6589           *i-- = 0;
6590           *i-- = 0;
6591           *i = 0;
6592         }
6593       else
6594         {
6595           *i++ = xi[M];
6596           *i++ = 0;
6597           *i++ = 0;
6598           *i = 0;
6599         }
6600     }
6601   /* Negate if negative */
6602   if (xi[0])
6603     {
6604       carry = 0;
6605       if (WORDS_BIG_ENDIAN)
6606         isave += 3;
6607       for (k = 0; k < 4; k++)
6608         {
6609           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6610           if (WORDS_BIG_ENDIAN)
6611             *isave-- = acc;
6612           else
6613             *isave++ = acc;
6614           carry = 0;
6615           if (acc & 0x10000)
6616             carry = 1;
6617         }
6618     }
6619 }
6620
6621
6622 /* Longhand square root routine.  */
6623
6624
6625 static int esqinited = 0;
6626 static unsigned short sqrndbit[NI];
6627
6628 static void
6629 esqrt (x, y)
6630      const UEMUSHORT *x;
6631      UEMUSHORT *y;
6632 {
6633   UEMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6634   EMULONG m, exp;
6635   int i, j, k, n, nlups;
6636
6637   if (esqinited == 0)
6638     {
6639       ecleaz (sqrndbit);
6640       sqrndbit[NI - 2] = 1;
6641       esqinited = 1;
6642     }
6643   /* Check for arg <= 0 */
6644   i = ecmp (x, ezero);
6645   if (i <= 0)
6646     {
6647       if (i == -1)
6648         {
6649           mtherr ("esqrt", DOMAIN);
6650           eclear (y);
6651         }
6652       else
6653         emov (x, y);
6654       return;
6655     }
6656
6657 #ifdef INFINITY
6658   if (eisinf (x))
6659     {
6660       eclear (y);
6661       einfin (y);
6662       return;
6663     }
6664 #endif
6665   /* Bring in the arg and renormalize if it is denormal.  */
6666   emovi (x, xx);
6667   m = (EMULONG) xx[1];          /* local long word exponent */
6668   if (m == 0)
6669     m -= enormlz (xx);
6670
6671   /* Divide exponent by 2 */
6672   m -= 0x3ffe;
6673   exp = (unsigned short) ((m / 2) + 0x3ffe);
6674
6675   /* Adjust if exponent odd */
6676   if ((m & 1) != 0)
6677     {
6678       if (m > 0)
6679         exp += 1;
6680       eshdn1 (xx);
6681     }
6682
6683   ecleaz (sq);
6684   ecleaz (num);
6685   n = 8;                        /* get 8 bits of result per inner loop */
6686   nlups = rndprc;
6687   j = 0;
6688
6689   while (nlups > 0)
6690     {
6691       /* bring in next word of arg */
6692       if (j < NE)
6693         num[NI - 1] = xx[j + 3];
6694       /* Do additional bit on last outer loop, for roundoff.  */
6695       if (nlups <= 8)
6696         n = nlups + 1;
6697       for (i = 0; i < n; i++)
6698         {
6699           /* Next 2 bits of arg */
6700           eshup1 (num);
6701           eshup1 (num);
6702           /* Shift up answer */
6703           eshup1 (sq);
6704           /* Make trial divisor */
6705           for (k = 0; k < NI; k++)
6706             temp[k] = sq[k];
6707           eshup1 (temp);
6708           eaddm (sqrndbit, temp);
6709           /* Subtract and insert answer bit if it goes in */
6710           if (ecmpm (temp, num) <= 0)
6711             {
6712               esubm (temp, num);
6713               sq[NI - 2] |= 1;
6714             }
6715         }
6716       nlups -= n;
6717       j += 1;
6718     }
6719
6720   /* Adjust for extra, roundoff loop done.  */
6721   exp += (NBITS - 1) - rndprc;
6722
6723   /* Sticky bit = 1 if the remainder is nonzero.  */
6724   k = 0;
6725   for (i = 3; i < NI; i++)
6726     k |= (int) num[i];
6727
6728   /* Renormalize and round off.  */
6729   emdnorm (sq, k, 0, exp, !ROUND_TOWARDS_ZERO);
6730   emovo (sq, y);
6731 }
6732 #endif
6733 \f
6734 /* Return the binary precision of the significand for a given
6735    floating point mode.  The mode can hold an integer value
6736    that many bits wide, without losing any bits.  */
6737
6738 unsigned int
6739 significand_size (mode)
6740      enum machine_mode mode;
6741 {
6742
6743   /* Don't test the modes, but their sizes, lest this
6744      code won't work for BITS_PER_UNIT != 8 .  */
6745
6746   switch (GET_MODE_BITSIZE (mode))
6747     {
6748      case 32:
6749
6750 #ifdef C4X
6751        return 56;
6752 #else
6753        return 24;
6754 #endif
6755
6756      case 64:
6757 #ifdef IEEE
6758        return 53;
6759 #else
6760        return 56;
6761 #endif
6762
6763      case 96:
6764        return 64;
6765
6766      case 128:
6767 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6768        return 113;
6769 #else
6770        return 64;
6771 #endif
6772
6773      default:
6774        abort ();
6775     }
6776 }