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