gcc/real.c

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