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