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   k = hex_value(*s);
5194   if ((k >= 0) && (k < base))
5195     {
5196       /* Ignore leading zeros */
5197       if ((prec == 0) && (decflg == 0) && (k == 0))
5198         goto donchr;
5199       /* Identify and strip trailing zeros after the decimal point.  */
5200       if ((trail == 0) && (decflg != 0))
5201         {
5202           sp = s;
5203           while (ISDIGIT (*sp) || (base == 16 && ISXDIGIT (*sp)))
5204             ++sp;
5205           /* Check for syntax error */
5206           c = *sp & CHARMASK;
5207           if ((base != 10 || ((c != 'e') && (c != 'E')))
5208               && (base != 16 || ((c != 'p') && (c != 'P')))
5209               && (c != '\0')
5210               && (c != '\n') && (c != '\r') && (c != ' ')
5211               && (c != ','))
5212             goto unexpected_char_error;
5213           --sp;
5214           while (*sp == '0')
5215             *sp-- = 'z';
5216           trail = 1;
5217           if (*s == 'z')
5218             goto donchr;
5219         }
5220
5221       /* If enough digits were given to more than fill up the yy register,
5222          continuing until overflow into the high guard word yy[2]
5223          guarantees that there will be a roundoff bit at the top
5224          of the low guard word after normalization.  */
5225
5226       if (yy[2] == 0)
5227         {
5228           if (base == 16)
5229             {
5230               if (decflg)
5231                 nexp += 4;      /* count digits after decimal point */
5232
5233               eshup1 (yy);      /* multiply current number by 16 */
5234               eshup1 (yy);
5235               eshup1 (yy);
5236               eshup1 (yy);
5237             }
5238           else
5239             {
5240               if (decflg)
5241                 nexp += 1;              /* count digits after decimal point */
5242
5243               eshup1 (yy);              /* multiply current number by 10 */
5244               emovz (yy, xt);
5245               eshup1 (xt);
5246               eshup1 (xt);
5247               eaddm (xt, yy);
5248             }
5249           /* Insert the current digit.  */
5250           ecleaz (xt);
5251           xt[NI - 2] = (UEMUSHORT) k;
5252           eaddm (xt, yy);
5253         }
5254       else
5255         {
5256           /* Mark any lost non-zero digit.  */
5257           lost |= k;
5258           /* Count lost digits before the decimal point.  */
5259           if (decflg == 0)
5260             {
5261               if (base == 10)
5262                 nexp -= 1;
5263               else
5264                 nexp -= 4;
5265             }
5266         }
5267       prec += 1;
5268       goto donchr;
5269     }
5270
5271   switch (*s)
5272     {
5273     case 'z':
5274       break;
5275     case 'E':
5276     case 'e':
5277     case 'P':
5278     case 'p':
5279       goto expnt;
5280     case '.':                   /* decimal point */
5281       if (decflg)
5282         goto unexpected_char_error;
5283       ++decflg;
5284       break;
5285     case '-':
5286       nsign = 0xffff;
5287       if (sgnflg)
5288         goto unexpected_char_error;
5289       ++sgnflg;
5290       break;
5291     case '+':
5292       if (sgnflg)
5293         goto unexpected_char_error;
5294       ++sgnflg;
5295       break;
5296     case ',':
5297     case ' ':
5298     case '\0':
5299     case '\n':
5300     case '\r':
5301       goto daldone;
5302     case 'i':
5303     case 'I':
5304       goto infinite;
5305     default:
5306     unexpected_char_error:
5307 #ifdef NANS
5308       einan (yy);
5309 #else
5310       mtherr ("asctoe", DOMAIN);
5311       eclear (yy);
5312 #endif
5313       goto aexit;
5314     }
5315  donchr:
5316   ++s;
5317   goto nxtcom;
5318
5319   /* Exponent interpretation */
5320  expnt:
5321   /* 0.0eXXX is zero, regardless of XXX.  Check for the 0.0.  */
5322   for (k = 0; k < NI; k++)
5323     {
5324       if (yy[k] != 0)
5325         goto read_expnt;
5326     }
5327   goto aexit;
5328
5329 read_expnt:
5330   esign = 1;
5331   exp = 0;
5332   ++s;
5333   /* check for + or - */
5334   if (*s == '-')
5335     {
5336       esign = -1;
5337       ++s;
5338     }
5339   if (*s == '+')
5340     ++s;
5341   while (ISDIGIT (*s))
5342     {
5343       exp *= 10;
5344       exp += *s++ - '0';
5345       if (exp > 999999)
5346         break;
5347     }
5348   if (esign < 0)
5349     exp = -exp;
5350   if ((exp > MAXDECEXP) && (base == 10))
5351     {
5352  infinite:
5353       ecleaz (yy);
5354       yy[E] = 0x7fff;           /* infinity */
5355       goto aexit;
5356     }
5357   if ((exp < MINDECEXP) && (base == 10))
5358     {
5359  zero:
5360       ecleaz (yy);
5361       goto aexit;
5362     }
5363
5364  daldone:
5365   if (base == 16)
5366     {
5367       /* Base 16 hexadecimal floating constant.  */
5368       if ((k = enormlz (yy)) > NBITS)
5369         {
5370           ecleaz (yy);
5371           goto aexit;
5372         }
5373       /* Adjust the exponent.  NEXP is the number of hex digits,
5374          EXP is a power of 2.  */
5375       lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5376       if (lexp > 0x7fff)
5377         goto infinite;
5378       if (lexp < 0)
5379         goto zero;
5380       yy[E] = lexp;
5381       goto expdon;
5382     }
5383
5384   nexp = exp - nexp;
5385   /* Pad trailing zeros to minimize power of 10, per IEEE spec.  */
5386   while ((nexp > 0) && (yy[2] == 0))
5387     {
5388       emovz (yy, xt);
5389       eshup1 (xt);
5390       eshup1 (xt);
5391       eaddm (yy, xt);
5392       eshup1 (xt);
5393       if (xt[2] != 0)
5394         break;
5395       nexp -= 1;
5396       emovz (xt, yy);
5397     }
5398   if ((k = enormlz (yy)) > NBITS)
5399     {
5400       ecleaz (yy);
5401       goto aexit;
5402     }
5403   lexp = (EXONE - 1 + NBITS) - k;
5404   emdnorm (yy, lost, 0, lexp, 64);
5405   lost = 0;
5406
5407   /* Convert to external format:
5408
5409      Multiply by 10**nexp.  If precision is 64 bits,
5410      the maximum relative error incurred in forming 10**n
5411      for 0 <= n <= 324 is 8.2e-20, at 10**180.
5412      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5413      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
5414
5415   lexp = yy[E];
5416   if (nexp == 0)
5417     {
5418       k = 0;
5419       goto expdon;
5420     }
5421   esign = 1;
5422   if (nexp < 0)
5423     {
5424       nexp = -nexp;
5425       esign = -1;
5426       if (nexp > 4096)
5427         {
5428           /* Punt.  Can't handle this without 2 divides.  */
5429           emovi (etens[0], tt);
5430           lexp -= tt[E];
5431           k = edivm (tt, yy);
5432           lexp += EXONE;
5433           nexp -= 4096;
5434         }
5435     }
5436   emov (eone, xt);
5437   exp = 1;
5438   i = NTEN;
5439   do
5440     {
5441       if (exp & nexp)
5442         emul (etens[i], xt, xt);
5443       i--;
5444       exp = exp + exp;
5445     }
5446   while (exp <= MAXP);
5447
5448   emovi (xt, tt);
5449   if (esign < 0)
5450     {
5451       lexp -= tt[E];
5452       k = edivm (tt, yy);
5453       lexp += EXONE;
5454     }
5455   else
5456     {
5457       lexp += tt[E];
5458       k = emulm (tt, yy);
5459       lexp -= EXONE - 1;
5460     }
5461   lost = k;
5462
5463  expdon:
5464
5465   /* Round and convert directly to the destination type */
5466   if (oprec == 53)
5467     lexp -= EXONE - 0x3ff;
5468 #ifdef C4X
5469   else if (oprec == 24 || oprec == 32)
5470     lexp -= (EXONE - 0x7f);
5471 #else
5472 #ifdef IBM
5473   else if (oprec == 24 || oprec == 56)
5474     lexp -= EXONE - (0x41 << 2);
5475 #else
5476   else if (oprec == 24)
5477     lexp -= EXONE - 0177;
5478 #endif /* IBM */
5479 #endif /* C4X */
5480 #ifdef DEC
5481   else if (oprec == 56)
5482     lexp -= EXONE - 0201;
5483 #endif
5484   rndprc = oprec;
5485   emdnorm (yy, lost, 0, lexp, 64);
5486
5487  aexit:
5488
5489   rndprc = rndsav;
5490   yy[0] = nsign;
5491   switch (oprec)
5492     {
5493 #ifdef DEC
5494     case 56:
5495       todec (yy, y);            /* see etodec.c */
5496       break;
5497 #endif
5498 #ifdef IBM
5499     case 56:
5500       toibm (yy, y, DFmode);
5501       break;
5502 #endif
5503 #ifdef C4X
5504     case 32:
5505       toc4x (yy, y, HFmode);
5506       break;
5507 #endif
5508
5509     case 53:
5510       toe53 (yy, y);
5511       break;
5512     case 24:
5513       toe24 (yy, y);
5514       break;
5515     case 64:
5516       toe64 (yy, y);
5517       break;
5518     case 113:
5519       toe113 (yy, y);
5520       break;
5521     case NBITS:
5522       emovo (yy, y);
5523       break;
5524     }
5525 }
5526
5527
5528
5529 /* Return Y = largest integer not greater than X (truncated toward minus
5530    infinity).  */
5531
5532 static const UEMUSHORT bmask[] =
5533 {
5534   0xffff,
5535   0xfffe,
5536   0xfffc,
5537   0xfff8,
5538   0xfff0,
5539   0xffe0,
5540   0xffc0,
5541   0xff80,
5542   0xff00,
5543   0xfe00,
5544   0xfc00,
5545   0xf800,
5546   0xf000,
5547   0xe000,
5548   0xc000,
5549   0x8000,
5550   0x0000,
5551 };
5552
5553 static void
5554 efloor (x, y)
5555      UEMUSHORT x[], y[];
5556 {
5557   UEMUSHORT *p;
5558   int e, expon, i;
5559   UEMUSHORT f[NE];
5560
5561   emov (x, f);                  /* leave in external format */
5562   expon = (int) f[NE - 1];
5563   e = (expon & 0x7fff) - (EXONE - 1);
5564   if (e <= 0)
5565     {
5566       eclear (y);
5567       goto isitneg;
5568     }
5569   /* number of bits to clear out */
5570   e = NBITS - e;
5571   emov (f, y);
5572   if (e <= 0)
5573     return;
5574
5575   p = &y[0];
5576   while (e >= 16)
5577     {
5578       *p++ = 0;
5579       e -= 16;
5580     }
5581   /* clear the remaining bits */
5582   *p &= bmask[e];
5583   /* truncate negatives toward minus infinity */
5584  isitneg:
5585
5586   if ((UEMUSHORT) expon & (UEMUSHORT) 0x8000)
5587     {
5588       for (i = 0; i < NE - 1; i++)
5589         {
5590           if (f[i] != y[i])
5591             {
5592               esub (eone, y, y);
5593               break;
5594             }
5595         }
5596     }
5597 }
5598
5599
5600 #if 0
5601 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5602    For example, 1.1 = 0.55 * 2^1.  */
5603
5604 static void
5605 efrexp (x, exp, s)
5606      UEMUSHORT x[];
5607      int *exp;
5608      UEMUSHORT s[];
5609 {
5610   UEMUSHORT xi[NI];
5611   EMULONG li;
5612
5613   emovi (x, xi);
5614   /*  Handle denormalized numbers properly using long integer exponent.  */
5615   li = (EMULONG) ((EMUSHORT) xi[1]);
5616
5617   if (li == 0)
5618     {
5619       li -= enormlz (xi);
5620     }
5621   xi[1] = 0x3ffe;
5622   emovo (xi, s);
5623   *exp = (int) (li - 0x3ffe);
5624 }
5625 #endif
5626
5627 /* Return e type Y = X * 2^PWR2.  */
5628
5629 static void
5630 eldexp (x, pwr2, y)
5631      UEMUSHORT x[];
5632      int pwr2;
5633      UEMUSHORT y[];
5634 {
5635   UEMUSHORT xi[NI];
5636   EMULONG li;
5637   int i;
5638
5639   emovi (x, xi);
5640   li = xi[1];
5641   li += pwr2;
5642   i = 0;
5643   emdnorm (xi, i, i, li, 64);
5644   emovo (xi, y);
5645 }
5646
5647
5648 #if 0
5649 /* C = remainder after dividing B by A, all e type values.
5650    Least significant integer quotient bits left in EQUOT.  */
5651
5652 static void
5653 eremain (a, b, c)
5654      UEMUSHORT a[], b[], c[];
5655 {
5656   UEMUSHORT den[NI], num[NI];
5657
5658 #ifdef NANS
5659   if (eisinf (b)
5660       || (ecmp (a, ezero) == 0)
5661       || eisnan (a)
5662       || eisnan (b))
5663     {
5664       enan (c, 0);
5665       return;
5666     }
5667 #endif
5668   if (ecmp (a, ezero) == 0)
5669     {
5670       mtherr ("eremain", SING);
5671       eclear (c);
5672       return;
5673     }
5674   emovi (a, den);
5675   emovi (b, num);
5676   eiremain (den, num);
5677   /* Sign of remainder = sign of quotient */
5678   if (a[0] == b[0])
5679     num[0] = 0;
5680   else
5681     num[0] = 0xffff;
5682   emovo (num, c);
5683 }
5684 #endif
5685
5686 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5687     remainder in NUM.  */
5688
5689 static void
5690 eiremain (den, num)
5691      UEMUSHORT den[], num[];
5692 {
5693   EMULONG ld, ln;
5694   UEMUSHORT j;
5695
5696   ld = den[E];
5697   ld -= enormlz (den);
5698   ln = num[E];
5699   ln -= enormlz (num);
5700   ecleaz (equot);
5701   while (ln >= ld)
5702     {
5703       if (ecmpm (den, num) <= 0)
5704         {
5705           esubm (den, num);
5706           j = 1;
5707         }
5708       else
5709           j = 0;
5710       eshup1 (equot);
5711       equot[NI - 1] |= j;
5712       eshup1 (num);
5713       ln -= 1;
5714     }
5715   emdnorm (num, 0, 0, ln, 0);
5716 }
5717
5718 /* Report an error condition CODE encountered in function NAME.
5719
5720     Mnemonic        Value          Significance
5721
5722      DOMAIN            1       argument domain error
5723      SING              2       function singularity
5724      OVERFLOW          3       overflow range error
5725      UNDERFLOW         4       underflow range error
5726      TLOSS             5       total loss of precision
5727      PLOSS             6       partial loss of precision
5728      INVALID           7       NaN - producing operation
5729      EDOM             33       Unix domain error code
5730      ERANGE           34       Unix range error code
5731
5732    The order of appearance of the following messages is bound to the
5733    error codes defined above.  */
5734
5735 int merror = 0;
5736 extern int merror;
5737
5738 static void
5739 mtherr (name, code)
5740      const char *name;
5741      int code;
5742 {
5743   /* The string passed by the calling program is supposed to be the
5744      name of the function in which the error occurred.
5745      The code argument selects which error message string will be printed.  */
5746
5747   if (strcmp (name, "esub") == 0)
5748     name = "subtraction";
5749   else if (strcmp (name, "ediv") == 0)
5750     name = "division";
5751   else if (strcmp (name, "emul") == 0)
5752     name = "multiplication";
5753   else if (strcmp (name, "enormlz") == 0)
5754     name = "normalization";
5755   else if (strcmp (name, "etoasc") == 0)
5756     name = "conversion to text";
5757   else if (strcmp (name, "asctoe") == 0)
5758     name = "parsing";
5759   else if (strcmp (name, "eremain") == 0)
5760     name = "modulus";
5761   else if (strcmp (name, "esqrt") == 0)
5762     name = "square root";
5763   if (extra_warnings)
5764     {
5765       switch (code)
5766         {
5767         case DOMAIN:    warning ("%s: argument domain error"    , name); break;
5768         case SING:      warning ("%s: function singularity"     , name); break;
5769         case OVERFLOW:  warning ("%s: overflow range error"     , name); break;
5770         case UNDERFLOW: warning ("%s: underflow range error"    , name); break;
5771         case TLOSS:     warning ("%s: total loss of precision"  , name); break;
5772         case PLOSS:     warning ("%s: partial loss of precision", name); break;
5773         case INVALID:   warning ("%s: NaN - producing operation", name); break;
5774         default:        abort ();
5775         }
5776     }
5777
5778   /* Set global error message word */
5779   merror = code + 1;
5780 }
5781
5782 #ifdef DEC
5783 /* Convert DEC double precision D to e type E.  */
5784
5785 static void
5786 dectoe (d, e)
5787      UEMUSHORT *d;
5788      UEMUSHORT *e;
5789 {
5790   UEMUSHORT y[NI];
5791   UEMUSHORT r, *p;
5792
5793   ecleaz (y);                   /* start with a zero */
5794   p = y;                        /* point to our number */
5795   r = *d;                       /* get DEC exponent word */
5796   if (*d & (unsigned int) 0x8000)
5797     *p = 0xffff;                /* fill in our sign */
5798   ++p;                          /* bump pointer to our exponent word */
5799   r &= 0x7fff;                  /* strip the sign bit */
5800   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
5801     goto done;
5802
5803
5804   r >>= 7;                      /* shift exponent word down 7 bits */
5805   r += EXONE - 0201;            /* subtract DEC exponent offset */
5806   /* add our e type exponent offset */
5807   *p++ = r;                     /* to form our exponent */
5808
5809   r = *d++;                     /* now do the high order mantissa */
5810   r &= 0177;                    /* strip off the DEC exponent and sign bits */
5811   r |= 0200;                    /* the DEC understood high order mantissa bit */
5812   *p++ = r;                     /* put result in our high guard word */
5813
5814   *p++ = *d++;                  /* fill in the rest of our mantissa */
5815   *p++ = *d++;
5816   *p = *d;
5817
5818   eshdn8 (y);                   /* shift our mantissa down 8 bits */
5819  done:
5820   emovo (y, e);
5821 }
5822
5823 /* Convert e type X to DEC double precision D.  */
5824
5825 static void
5826 etodec (x, d)
5827      UEMUSHORT *x, *d;
5828 {
5829   UEMUSHORT xi[NI];
5830   EMULONG exp;
5831   int rndsav;
5832
5833   emovi (x, xi);
5834   /* Adjust exponent for offsets.  */
5835   exp = (EMULONG) xi[E] - (EXONE - 0201);
5836   /* Round off to nearest or even.  */
5837   rndsav = rndprc;
5838   rndprc = 56;
5839   emdnorm (xi, 0, 0, exp, 64);
5840   rndprc = rndsav;
5841   todec (xi, d);
5842 }
5843
5844 /* Convert exploded e-type X, that has already been rounded to
5845    56-bit precision, to DEC format double Y.  */
5846
5847 static void
5848 todec (x, y)
5849      UEMUSHORT *x, *y;
5850 {
5851   UEMUSHORT i;
5852   UEMUSHORT *p;
5853
5854   p = x;
5855   *y = 0;
5856   if (*p++)
5857     *y = 0100000;
5858   i = *p++;
5859   if (i == 0)
5860     {
5861       *y++ = 0;
5862       *y++ = 0;
5863       *y++ = 0;
5864       *y++ = 0;
5865       return;
5866     }
5867   if (i > 0377)
5868     {
5869       *y++ |= 077777;
5870       *y++ = 0xffff;
5871       *y++ = 0xffff;
5872       *y++ = 0xffff;
5873 #ifdef ERANGE
5874       errno = ERANGE;
5875 #endif
5876       return;
5877     }
5878   i &= 0377;
5879   i <<= 7;
5880   eshup8 (x);
5881   x[M] &= 0177;
5882   i |= x[M];
5883   *y++ |= i;
5884   *y++ = x[M + 1];
5885   *y++ = x[M + 2];
5886   *y++ = x[M + 3];
5887 }
5888 #endif /* DEC */
5889
5890 #ifdef IBM
5891 /* Convert IBM single/double precision to e type.  */
5892
5893 static void
5894 ibmtoe (d, e, mode)
5895      UEMUSHORT *d;
5896      UEMUSHORT *e;
5897      enum machine_mode mode;
5898 {
5899   UEMUSHORT y[NI];
5900   UEMUSHORT r, *p;
5901
5902   ecleaz (y);                   /* start with a zero */
5903   p = y;                        /* point to our number */
5904   r = *d;                       /* get IBM exponent word */
5905   if (*d & (unsigned int) 0x8000)
5906     *p = 0xffff;                /* fill in our sign */
5907   ++p;                          /* bump pointer to our exponent word */
5908   r &= 0x7f00;                  /* strip the sign bit */
5909   r >>= 6;                      /* shift exponent word down 6 bits */
5910                                 /* in fact shift by 8 right and 2 left */
5911   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5912                                 /* add our e type exponent offset */
5913   *p++ = r;                     /* to form our exponent */
5914
5915   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5916                                 /* strip off the IBM exponent and sign bits */
5917   if (mode != SFmode)           /* there are only 2 words in SFmode */
5918     {
5919       *p++ = *d++;              /* fill in the rest of our mantissa */
5920       *p++ = *d++;
5921     }
5922   *p = *d;
5923
5924   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5925     y[0] = y[E] = 0;
5926   else
5927     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5928                               /* handle change in RADIX */
5929   emovo (y, e);
5930 }
5931
5932
5933
5934 /* Convert e type to IBM single/double precision.  */
5935
5936 static void
5937 etoibm (x, d, mode)
5938      UEMUSHORT *x, *d;
5939      enum machine_mode mode;
5940 {
5941   UEMUSHORT xi[NI];
5942   EMULONG exp;
5943   int rndsav;
5944
5945   emovi (x, xi);
5946   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5947                                                         /* round off to nearest or even */
5948   rndsav = rndprc;
5949   rndprc = 56;
5950   emdnorm (xi, 0, 0, exp, 64);
5951   rndprc = rndsav;
5952   toibm (xi, d, mode);
5953 }
5954
5955 static void
5956 toibm (x, y, mode)
5957      UEMUSHORT *x, *y;
5958      enum machine_mode mode;
5959 {
5960   UEMUSHORT i;
5961   UEMUSHORT *p;
5962   int r;
5963
5964   p = x;
5965   *y = 0;
5966   if (*p++)
5967     *y = 0x8000;
5968   i = *p++;
5969   if (i == 0)
5970     {
5971       *y++ = 0;
5972       *y++ = 0;
5973       if (mode != SFmode)
5974         {
5975           *y++ = 0;
5976           *y++ = 0;
5977         }
5978       return;
5979     }
5980   r = i & 0x3;
5981   i >>= 2;
5982   if (i > 0x7f)
5983     {
5984       *y++ |= 0x7fff;
5985       *y++ = 0xffff;
5986       if (mode != SFmode)
5987         {
5988           *y++ = 0xffff;
5989           *y++ = 0xffff;
5990         }
5991 #ifdef ERANGE
5992       errno = ERANGE;
5993 #endif
5994       return;
5995     }
5996   i &= 0x7f;
5997   *y |= (i << 8);
5998   eshift (x, r + 5);
5999   *y++ |= x[M];
6000   *y++ = x[M + 1];
6001   if (mode != SFmode)
6002     {
6003       *y++ = x[M + 2];
6004       *y++ = x[M + 3];
6005     }
6006 }
6007 #endif /* IBM */
6008
6009
6010 #ifdef C4X
6011 /* Convert C4X single/double precision to e type.  */
6012
6013 static void
6014 c4xtoe (d, e, mode)
6015      UEMUSHORT *d;
6016      UEMUSHORT *e;
6017      enum machine_mode mode;
6018 {
6019   UEMUSHORT y[NI];
6020   int r;
6021   int isnegative;
6022   int size;
6023   int i;
6024   int carry;
6025
6026   /* Short-circuit the zero case.  */
6027   if ((d[0] == 0x8000)
6028       && (d[1] == 0x0000)
6029       && ((mode == QFmode) || ((d[2] == 0x0000) && (d[3] == 0x0000))))
6030     {
6031       e[0] = 0;
6032       e[1] = 0;
6033       e[2] = 0;
6034       e[3] = 0;
6035       e[4] = 0;
6036       e[5] = 0;
6037       return;
6038     }
6039
6040   ecleaz (y);                   /* start with a zero */
6041   r = d[0];                     /* get sign/exponent part */
6042   if (r & (unsigned int) 0x0080)
6043   {
6044      y[0] = 0xffff;             /* fill in our sign */
6045      isnegative = TRUE;
6046   }
6047   else
6048   {
6049      isnegative = FALSE;
6050   }
6051
6052   r >>= 8;                      /* Shift exponent word down 8 bits.  */
6053   if (r & 0x80)                 /* Make the exponent negative if it is.  */
6054   {
6055      r = r | (~0 & ~0xff);
6056   }
6057
6058   if (isnegative)
6059   {
6060      /* Now do the high order mantissa.  We don't "or" on the high bit
6061         because it is 2 (not 1) and is handled a little differently
6062         below.  */
6063      y[M] = d[0] & 0x7f;
6064
6065      y[M+1] = d[1];
6066      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
6067      {
6068         y[M+2] = d[2];          /* Fill in the rest of our mantissa.  */
6069         y[M+3] = d[3];
6070         size = 4;
6071      }
6072      else
6073      {
6074         size = 2;
6075      }
6076      eshift(y, -8);
6077
6078      /* Now do the two's complement on the data.  */
6079
6080      carry = 1; /* Initially add 1 for the two's complement.  */
6081      for (i=size + M; i > M; i--)
6082      {
6083         if (carry && (y[i] == 0x0000))
6084         {
6085            /* We overflowed into the next word, carry is the same.  */
6086            y[i] = carry ? 0x0000 : 0xffff;
6087         }
6088         else
6089         {
6090            /* No overflow, just invert and add carry.  */
6091            y[i] = ((~y[i]) + carry) & 0xffff;
6092            carry = 0;
6093         }
6094      }
6095
6096      if (carry)
6097      {
6098         eshift(y, -1);
6099         y[M+1] |= 0x8000;
6100         r++;
6101      }
6102      y[1] = r + EXONE;
6103   }
6104   else
6105   {
6106     /* Add our e type exponent offset to form our exponent.  */
6107      r += EXONE;
6108      y[1] = r;
6109
6110      /* Now do the high order mantissa strip off the exponent and sign
6111         bits and add the high 1 bit.  */
6112      y[M] = (d[0] & 0x7f) | 0x80;
6113
6114      y[M+1] = d[1];
6115      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
6116      {
6117         y[M+2] = d[2];          /* Fill in the rest of our mantissa.  */
6118         y[M+3] = d[3];
6119      }
6120      eshift(y, -8);
6121   }
6122
6123   emovo (y, e);
6124 }
6125
6126
6127 /* Convert e type to C4X single/double precision.  */
6128
6129 static void
6130 etoc4x (x, d, mode)
6131      UEMUSHORT *x, *d;
6132      enum machine_mode mode;
6133 {
6134   UEMUSHORT xi[NI];
6135   EMULONG exp;
6136   int rndsav;
6137
6138   emovi (x, xi);
6139
6140   /* Adjust exponent for offsets.  */
6141   exp = (EMULONG) xi[E] - (EXONE - 0x7f);
6142
6143   /* Round off to nearest or even.  */
6144   rndsav = rndprc;
6145   rndprc = mode == QFmode ? 24 : 32;
6146   emdnorm (xi, 0, 0, exp, 64);
6147   rndprc = rndsav;
6148   toc4x (xi, d, mode);
6149 }
6150
6151 static void
6152 toc4x (x, y, mode)
6153      UEMUSHORT *x, *y;
6154      enum machine_mode mode;
6155 {
6156   int i;
6157   int v;
6158   int carry;
6159
6160   /* Short-circuit the zero case */
6161   if ((x[0] == 0)       /* Zero exponent and sign */
6162       && (x[1] == 0)
6163       && (x[M] == 0)    /* The rest is for zero mantissa */
6164       && (x[M+1] == 0)
6165       /* Only check for double if necessary */
6166       && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
6167     {
6168       /* We have a zero.  Put it into the output and return.  */
6169       *y++ = 0x8000;
6170       *y++ = 0x0000;
6171       if (mode != QFmode)
6172         {
6173           *y++ = 0x0000;
6174           *y++ = 0x0000;
6175         }
6176       return;
6177     }
6178
6179   *y = 0;
6180
6181   /* Negative number require a two's complement conversion of the
6182      mantissa.  */
6183   if (x[0])
6184     {
6185       *y = 0x0080;
6186
6187       i = ((int) x[1]) - 0x7f;
6188
6189       /* Now add 1 to the inverted data to do the two's complement.  */
6190       if (mode != QFmode)
6191         v = 4 + M;
6192       else
6193         v = 2 + M;
6194       carry = 1;
6195       while (v > M)
6196         {
6197           if (x[v] == 0x0000)
6198             {
6199               x[v] = carry ? 0x0000 : 0xffff;
6200             }
6201           else
6202             {
6203               x[v] = ((~x[v]) + carry) & 0xffff;
6204               carry = 0;
6205             }
6206           v--;
6207         }
6208
6209       /* The following is a special case.  The C4X negative float requires
6210          a zero in the high bit (because the format is (2 - x) x 2^m), so
6211          if a one is in that bit, we have to shift left one to get rid
6212          of it.  This only occurs if the number is -1 x 2^m.  */
6213       if (x[M+1] & 0x8000)
6214         {
6215           /* This is the case of -1 x 2^m, we have to rid ourselves of the
6216              high sign bit and shift the exponent.  */
6217           eshift(x, 1);
6218           i--;
6219         }
6220     }
6221   else
6222     {
6223       i = ((int) x[1]) - 0x7f;
6224     }
6225
6226   if ((i < -128) || (i > 127))
6227     {
6228       y[0] |= 0xff7f;
6229       y[1] = 0xffff;
6230       if (mode != QFmode)
6231         {
6232           y[2] = 0xffff;
6233           y[3] = 0xffff;
6234         }
6235 #ifdef ERANGE
6236       errno = ERANGE;
6237 #endif
6238       return;
6239     }
6240
6241   y[0] |= ((i & 0xff) << 8);
6242
6243   eshift (x, 8);
6244
6245   y[0] |= x[M] & 0x7f;
6246   y[1] = x[M + 1];
6247   if (mode != QFmode)
6248     {
6249       y[2] = x[M + 2];
6250       y[3] = x[M + 3];
6251     }
6252 }
6253 #endif /* C4X */
6254
6255 /* Output a binary NaN bit pattern in the target machine's format.  */
6256
6257 /* If special NaN bit patterns are required, define them in tm.h
6258    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
6259    patterns here.  */
6260 #ifdef TFMODE_NAN
6261 TFMODE_NAN;
6262 #else
6263 #ifdef IEEE
6264 UEMUSHORT TFbignan[8] =
6265  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6266 UEMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6267 #endif
6268 #endif
6269
6270 #ifdef XFMODE_NAN
6271 XFMODE_NAN;
6272 #else
6273 #ifdef IEEE
6274 UEMUSHORT XFbignan[6] =
6275  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6276 UEMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6277 #endif
6278 #endif
6279
6280 #ifdef DFMODE_NAN
6281 DFMODE_NAN;
6282 #else
6283 #ifdef IEEE
6284 UEMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6285 UEMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
6286 #endif
6287 #endif
6288
6289 #ifdef SFMODE_NAN
6290 SFMODE_NAN;
6291 #else
6292 #ifdef IEEE
6293 UEMUSHORT SFbignan[2] = {0x7fff, 0xffff};
6294 UEMUSHORT SFlittlenan[2] = {0, 0xffc0};
6295 #endif
6296 #endif
6297
6298
6299 #ifdef NANS
6300 static void
6301 make_nan (nan, sign, mode)
6302      UEMUSHORT *nan;
6303      int sign;
6304      enum machine_mode mode;
6305 {
6306   int n;
6307   UEMUSHORT *p;
6308
6309   switch (mode)
6310     {
6311 /* Possibly the `reserved operand' patterns on a VAX can be
6312    used like NaN's, but probably not in the same way as IEEE.  */
6313 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6314     case TFmode:
6315 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6316       n = 8;
6317       if (REAL_WORDS_BIG_ENDIAN)
6318         p = TFbignan;
6319       else
6320         p = TFlittlenan;
6321       break;
6322 #endif
6323       /* FALLTHRU */
6324
6325     case XFmode:
6326       n = 6;
6327       if (REAL_WORDS_BIG_ENDIAN)
6328         p = XFbignan;
6329       else
6330         p = XFlittlenan;
6331       break;
6332
6333     case DFmode:
6334       n = 4;
6335       if (REAL_WORDS_BIG_ENDIAN)
6336         p = DFbignan;
6337       else
6338         p = DFlittlenan;
6339       break;
6340
6341     case SFmode:
6342     case HFmode:
6343       n = 2;
6344       if (REAL_WORDS_BIG_ENDIAN)
6345         p = SFbignan;
6346       else
6347         p = SFlittlenan;
6348       break;
6349 #endif
6350
6351     default:
6352       abort ();
6353     }
6354   if (REAL_WORDS_BIG_ENDIAN)
6355     *nan++ = (sign << 15) | (*p++ & 0x7fff);
6356   while (--n != 0)
6357     *nan++ = *p++;
6358   if (! REAL_WORDS_BIG_ENDIAN)
6359     *nan = (sign << 15) | (*p & 0x7fff);
6360 }
6361 #endif /* NANS */
6362
6363 /* This is the inverse of the function `etarsingle' invoked by
6364    REAL_VALUE_TO_TARGET_SINGLE.  */
6365
6366 REAL_VALUE_TYPE
6367 ereal_unto_float (f)
6368      long f;
6369 {
6370   REAL_VALUE_TYPE r;
6371   UEMUSHORT s[2];
6372   UEMUSHORT e[NE];
6373
6374   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6375    This is the inverse operation to what the function `endian' does.  */
6376   if (REAL_WORDS_BIG_ENDIAN)
6377     {
6378       s[0] = (UEMUSHORT) (f >> 16);
6379       s[1] = (UEMUSHORT) f;
6380     }
6381   else
6382     {
6383       s[0] = (UEMUSHORT) f;
6384       s[1] = (UEMUSHORT) (f >> 16);
6385     }
6386   /* Convert and promote the target float to E-type.  */
6387   e24toe (s, e);
6388   /* Output E-type to REAL_VALUE_TYPE.  */
6389   PUT_REAL (e, &r);
6390   return r;
6391 }
6392
6393
6394 /* This is the inverse of the function `etardouble' invoked by
6395    REAL_VALUE_TO_TARGET_DOUBLE.  */
6396
6397 REAL_VALUE_TYPE
6398 ereal_unto_double (d)
6399      long d[];
6400 {
6401   REAL_VALUE_TYPE r;
6402   UEMUSHORT s[4];
6403   UEMUSHORT e[NE];
6404
6405   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6406   if (REAL_WORDS_BIG_ENDIAN)
6407     {
6408       s[0] = (UEMUSHORT) (d[0] >> 16);
6409       s[1] = (UEMUSHORT) d[0];
6410       s[2] = (UEMUSHORT) (d[1] >> 16);
6411       s[3] = (UEMUSHORT) d[1];
6412     }
6413   else
6414     {
6415       /* Target float words are little-endian.  */
6416       s[0] = (UEMUSHORT) d[0];
6417       s[1] = (UEMUSHORT) (d[0] >> 16);
6418       s[2] = (UEMUSHORT) d[1];
6419       s[3] = (UEMUSHORT) (d[1] >> 16);
6420     }
6421   /* Convert target double to E-type.  */
6422   e53toe (s, e);
6423   /* Output E-type to REAL_VALUE_TYPE.  */
6424   PUT_REAL (e, &r);
6425   return r;
6426 }
6427
6428
6429 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6430    This is somewhat like ereal_unto_float, but the input types
6431    for these are different.  */
6432
6433 REAL_VALUE_TYPE
6434 ereal_from_float (f)
6435      HOST_WIDE_INT f;
6436 {
6437   REAL_VALUE_TYPE r;
6438   UEMUSHORT s[2];
6439   UEMUSHORT e[NE];
6440
6441   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6442    This is the inverse operation to what the function `endian' does.  */
6443   if (REAL_WORDS_BIG_ENDIAN)
6444     {
6445       s[0] = (UEMUSHORT) (f >> 16);
6446       s[1] = (UEMUSHORT) f;
6447     }
6448   else
6449     {
6450       s[0] = (UEMUSHORT) f;
6451       s[1] = (UEMUSHORT) (f >> 16);
6452     }
6453   /* Convert and promote the target float to E-type.  */
6454   e24toe (s, e);
6455   /* Output E-type to REAL_VALUE_TYPE.  */
6456   PUT_REAL (e, &r);
6457   return r;
6458 }
6459
6460
6461 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6462    This is somewhat like ereal_unto_double, but the input types
6463    for these are different.
6464
6465    The DFmode is stored as an array of HOST_WIDE_INT in the target's
6466    data format, with no holes in the bit packing.  The first element
6467    of the input array holds the bits that would come first in the
6468    target computer's memory.  */
6469
6470 REAL_VALUE_TYPE
6471 ereal_from_double (d)
6472      HOST_WIDE_INT d[];
6473 {
6474   REAL_VALUE_TYPE r;
6475   UEMUSHORT s[4];
6476   UEMUSHORT e[NE];
6477
6478   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6479   if (REAL_WORDS_BIG_ENDIAN)
6480     {
6481 #if HOST_BITS_PER_WIDE_INT == 32
6482       s[0] = (UEMUSHORT) (d[0] >> 16);
6483       s[1] = (UEMUSHORT) d[0];
6484       s[2] = (UEMUSHORT) (d[1] >> 16);
6485       s[3] = (UEMUSHORT) d[1];
6486 #else
6487       /* In this case the entire target double is contained in the
6488          first array element.  The second element of the input is
6489          ignored.  */
6490       s[0] = (UEMUSHORT) (d[0] >> 48);
6491       s[1] = (UEMUSHORT) (d[0] >> 32);
6492       s[2] = (UEMUSHORT) (d[0] >> 16);
6493       s[3] = (UEMUSHORT) d[0];
6494 #endif
6495     }
6496   else
6497     {
6498       /* Target float words are little-endian.  */
6499       s[0] = (UEMUSHORT) d[0];
6500       s[1] = (UEMUSHORT) (d[0] >> 16);
6501 #if HOST_BITS_PER_WIDE_INT == 32
6502       s[2] = (UEMUSHORT) d[1];
6503       s[3] = (UEMUSHORT) (d[1] >> 16);
6504 #else
6505       s[2] = (UEMUSHORT) (d[0] >> 32);
6506       s[3] = (UEMUSHORT) (d[0] >> 48);
6507 #endif
6508     }
6509   /* Convert target double to E-type.  */
6510   e53toe (s, e);
6511   /* Output E-type to REAL_VALUE_TYPE.  */
6512   PUT_REAL (e, &r);
6513   return r;
6514 }
6515
6516
6517 #if 0
6518 /* Convert target computer unsigned 64-bit integer to e-type.
6519    The endian-ness of DImode follows the convention for integers,
6520    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
6521
6522 static void
6523 uditoe (di, e)
6524      UEMUSHORT *di;  /* Address of the 64-bit int.  */
6525      UEMUSHORT *e;
6526 {
6527   UEMUSHORT yi[NI];
6528   int k;
6529
6530   ecleaz (yi);
6531   if (WORDS_BIG_ENDIAN)
6532     {
6533       for (k = M; k < M + 4; k++)
6534         yi[k] = *di++;
6535     }
6536   else
6537     {
6538       for (k = M + 3; k >= M; k--)
6539         yi[k] = *di++;
6540     }
6541   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6542   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6543     ecleaz (yi);                /* it was zero */
6544   else
6545     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6546   emovo (yi, e);
6547 }
6548
6549 /* Convert target computer signed 64-bit integer to e-type.  */
6550
6551 static void
6552 ditoe (di, e)
6553      UEMUSHORT *di;  /* Address of the 64-bit int.  */
6554      UEMUSHORT *e;
6555 {
6556   unsigned EMULONG acc;
6557   UEMUSHORT yi[NI];
6558   UEMUSHORT carry;
6559   int k, sign;
6560
6561   ecleaz (yi);
6562   if (WORDS_BIG_ENDIAN)
6563     {
6564       for (k = M; k < M + 4; k++)
6565         yi[k] = *di++;
6566     }
6567   else
6568     {
6569       for (k = M + 3; k >= M; k--)
6570         yi[k] = *di++;
6571     }
6572   /* Take absolute value */
6573   sign = 0;
6574   if (yi[M] & 0x8000)
6575     {
6576       sign = 1;
6577       carry = 0;
6578       for (k = M + 3; k >= M; k--)
6579         {
6580           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6581           yi[k] = acc;
6582           carry = 0;
6583           if (acc & 0x10000)
6584             carry = 1;
6585         }
6586     }
6587   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6588   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6589     ecleaz (yi);                /* it was zero */
6590   else
6591     yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6592   emovo (yi, e);
6593   if (sign)
6594         eneg (e);
6595 }
6596
6597
6598 /* Convert e-type to unsigned 64-bit int.  */
6599
6600 static void
6601 etoudi (x, i)
6602      UEMUSHORT *x;
6603      UEMUSHORT *i;
6604 {
6605   UEMUSHORT xi[NI];
6606   int j, k;
6607
6608   emovi (x, xi);
6609   if (xi[0])
6610     {
6611       xi[M] = 0;
6612       goto noshift;
6613     }
6614   k = (int) xi[E] - (EXONE - 1);
6615   if (k <= 0)
6616     {
6617       for (j = 0; j < 4; j++)
6618         *i++ = 0;
6619       return;
6620     }
6621   if (k > 64)
6622     {
6623       for (j = 0; j < 4; j++)
6624         *i++ = 0xffff;
6625       if (extra_warnings)
6626         warning ("overflow on truncation to integer");
6627       return;
6628     }
6629   if (k > 16)
6630     {
6631       /* Shift more than 16 bits: first shift up k-16 mod 16,
6632          then shift up by 16's.  */
6633       j = k - ((k >> 4) << 4);
6634       if (j == 0)
6635         j = 16;
6636       eshift (xi, j);
6637       if (WORDS_BIG_ENDIAN)
6638         *i++ = xi[M];
6639       else
6640         {
6641           i += 3;
6642           *i-- = xi[M];
6643         }
6644       k -= j;
6645       do
6646         {
6647           eshup6 (xi);
6648           if (WORDS_BIG_ENDIAN)
6649             *i++ = xi[M];
6650           else
6651             *i-- = xi[M];
6652         }
6653       while ((k -= 16) > 0);
6654     }
6655   else
6656     {
6657         /* shift not more than 16 bits */
6658       eshift (xi, k);
6659
6660 noshift:
6661
6662       if (WORDS_BIG_ENDIAN)
6663         {
6664           i += 3;
6665           *i-- = xi[M];
6666           *i-- = 0;
6667           *i-- = 0;
6668           *i = 0;
6669         }
6670       else
6671         {
6672           *i++ = xi[M];
6673           *i++ = 0;
6674           *i++ = 0;
6675           *i = 0;
6676         }
6677     }
6678 }
6679
6680
6681 /* Convert e-type to signed 64-bit int.  */
6682
6683 static void
6684 etodi (x, i)
6685      UEMUSHORT *x;
6686      UEMUSHORT *i;
6687 {
6688   unsigned EMULONG acc;
6689   UEMUSHORT xi[NI];
6690   UEMUSHORT carry;
6691   UEMUSHORT *isave;
6692   int j, k;
6693
6694   emovi (x, xi);
6695   k = (int) xi[E] - (EXONE - 1);
6696   if (k <= 0)
6697     {
6698       for (j = 0; j < 4; j++)
6699         *i++ = 0;
6700       return;
6701     }
6702   if (k > 64)
6703     {
6704       for (j = 0; j < 4; j++)
6705         *i++ = 0xffff;
6706       if (extra_warnings)
6707         warning ("overflow on truncation to integer");
6708       return;
6709     }
6710   isave = i;
6711   if (k > 16)
6712     {
6713       /* Shift more than 16 bits: first shift up k-16 mod 16,
6714          then shift up by 16's.  */
6715       j = k - ((k >> 4) << 4);
6716       if (j == 0)
6717         j = 16;
6718       eshift (xi, j);
6719       if (WORDS_BIG_ENDIAN)
6720         *i++ = xi[M];
6721       else
6722         {
6723           i += 3;
6724           *i-- = xi[M];
6725         }
6726       k -= j;
6727       do
6728         {
6729           eshup6 (xi);
6730           if (WORDS_BIG_ENDIAN)
6731             *i++ = xi[M];
6732           else
6733             *i-- = xi[M];
6734         }
6735       while ((k -= 16) > 0);
6736     }
6737   else
6738     {
6739         /* shift not more than 16 bits */
6740       eshift (xi, k);
6741
6742       if (WORDS_BIG_ENDIAN)
6743         {
6744           i += 3;
6745           *i = xi[M];
6746           *i-- = 0;
6747           *i-- = 0;
6748           *i = 0;
6749         }
6750       else
6751         {
6752           *i++ = xi[M];
6753           *i++ = 0;
6754           *i++ = 0;
6755           *i = 0;
6756         }
6757     }
6758   /* Negate if negative */
6759   if (xi[0])
6760     {
6761       carry = 0;
6762       if (WORDS_BIG_ENDIAN)
6763         isave += 3;
6764       for (k = 0; k < 4; k++)
6765         {
6766           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6767           if (WORDS_BIG_ENDIAN)
6768             *isave-- = acc;
6769           else
6770             *isave++ = acc;
6771           carry = 0;
6772           if (acc & 0x10000)
6773             carry = 1;
6774         }
6775     }
6776 }
6777
6778
6779 /* Longhand square root routine.  */
6780
6781
6782 static int esqinited = 0;
6783 static unsigned short sqrndbit[NI];
6784
6785 static void
6786 esqrt (x, y)
6787      UEMUSHORT *x, *y;
6788 {
6789   UEMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6790   EMULONG m, exp;
6791   int i, j, k, n, nlups;
6792
6793   if (esqinited == 0)
6794     {
6795       ecleaz (sqrndbit);
6796       sqrndbit[NI - 2] = 1;
6797       esqinited = 1;
6798     }
6799   /* Check for arg <= 0 */
6800   i = ecmp (x, ezero);
6801   if (i <= 0)
6802     {
6803       if (i == -1)
6804         {
6805           mtherr ("esqrt", DOMAIN);
6806           eclear (y);
6807         }
6808       else
6809         emov (x, y);
6810       return;
6811     }
6812
6813 #ifdef INFINITY
6814   if (eisinf (x))
6815     {
6816       eclear (y);
6817       einfin (y);
6818       return;
6819     }
6820 #endif
6821   /* Bring in the arg and renormalize if it is denormal.  */
6822   emovi (x, xx);
6823   m = (EMULONG) xx[1];          /* local long word exponent */
6824   if (m == 0)
6825     m -= enormlz (xx);
6826
6827   /* Divide exponent by 2 */
6828   m -= 0x3ffe;
6829   exp = (unsigned short) ((m / 2) + 0x3ffe);
6830
6831   /* Adjust if exponent odd */
6832   if ((m & 1) != 0)
6833     {
6834       if (m > 0)
6835         exp += 1;
6836       eshdn1 (xx);
6837     }
6838
6839   ecleaz (sq);
6840   ecleaz (num);
6841   n = 8;                        /* get 8 bits of result per inner loop */
6842   nlups = rndprc;
6843   j = 0;
6844
6845   while (nlups > 0)
6846     {
6847       /* bring in next word of arg */
6848       if (j < NE)
6849         num[NI - 1] = xx[j + 3];
6850       /* Do additional bit on last outer loop, for roundoff.  */
6851       if (nlups <= 8)
6852         n = nlups + 1;
6853       for (i = 0; i < n; i++)
6854         {
6855           /* Next 2 bits of arg */
6856           eshup1 (num);
6857           eshup1 (num);
6858           /* Shift up answer */
6859           eshup1 (sq);
6860           /* Make trial divisor */
6861           for (k = 0; k < NI; k++)
6862             temp[k] = sq[k];
6863           eshup1 (temp);
6864           eaddm (sqrndbit, temp);
6865           /* Subtract and insert answer bit if it goes in */
6866           if (ecmpm (temp, num) <= 0)
6867             {
6868               esubm (temp, num);
6869               sq[NI - 2] |= 1;
6870             }
6871         }
6872       nlups -= n;
6873       j += 1;
6874     }
6875
6876   /* Adjust for extra, roundoff loop done.  */
6877   exp += (NBITS - 1) - rndprc;
6878
6879   /* Sticky bit = 1 if the remainder is nonzero.  */
6880   k = 0;
6881   for (i = 3; i < NI; i++)
6882     k |= (int) num[i];
6883
6884   /* Renormalize and round off.  */
6885   emdnorm (sq, k, 0, exp, 64);
6886   emovo (sq, y);
6887 }
6888 #endif
6889 #endif /* EMU_NON_COMPILE not defined */
6890 \f
6891 /* Return the binary precision of the significand for a given
6892    floating point mode.  The mode can hold an integer value
6893    that many bits wide, without losing any bits.  */
6894
6895 unsigned int
6896 significand_size (mode)
6897      enum machine_mode mode;
6898 {
6899
6900 /* Don't test the modes, but their sizes, lest this
6901    code won't work for BITS_PER_UNIT != 8 .  */
6902
6903 switch (GET_MODE_BITSIZE (mode))
6904   {
6905   case 32:
6906
6907 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6908     return 56;
6909 #endif
6910
6911     return 24;
6912
6913   case 64:
6914 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6915     return 53;
6916 #else
6917 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6918     return 56;
6919 #else
6920 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6921     return 56;
6922 #else
6923 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6924     return 56;
6925 #else
6926     abort ();
6927 #endif
6928 #endif
6929 #endif
6930 #endif
6931
6932   case 96:
6933     return 64;
6934
6935   case 128:
6936 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6937     return 113;
6938 #else
6939     return 64;
6940 #endif
6941
6942   default:
6943     abort ();
6944   }
6945 }