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