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