gcc/real.c

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