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