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       /* The purpose of this conditional is to avoid scribbling beyond
3650          the end of a long double, in case the type is only 80 bits wide.  */
3651       if (LONG_DOUBLE_TYPE_SIZE == 96)
3652         {
3653           /* Clear the last two bytes of 12-byte Intel format */
3654           *(q+1) = 0;
3655         }
3656 #endif
3657     }
3658 #endif
3659
3660   /* combine sign and exponent */
3661   i = *p++;
3662 #ifdef IBM
3663   if (i)
3664     *q++ = *p++ | 0x8000;
3665   else
3666     *q++ = *p++;
3667   *q++ = 0;
3668 #endif
3669 #ifdef DEC
3670   if (i)
3671     *q-- = *p++ | 0x8000;
3672   else
3673     *q-- = *p++;
3674 #endif
3675 #ifdef IEEE
3676   if (REAL_WORDS_BIG_ENDIAN)
3677     {
3678 #ifdef ARM_EXTENDED_IEEE_FORMAT
3679       /* The exponent is in the lowest 15 bits of the first word.  */
3680       *q++ = i ? 0x8000 : 0;
3681       *q++ = *p++;
3682 #else
3683       if (i)
3684         *q++ = *p++ | 0x8000;
3685       else
3686         *q++ = *p++;
3687       *q++ = 0;
3688 #endif
3689     }
3690   else
3691     {
3692       if (i)
3693         *q-- = *p++ | 0x8000;
3694       else
3695         *q-- = *p++;
3696     }
3697 #endif
3698   /* skip over guard word */
3699   ++p;
3700   /* move the significand */
3701 #ifdef IBM
3702   for (i = 0; i < 4; i++)
3703     *q++ = *p++;
3704 #endif
3705 #ifdef DEC
3706   for (i = 0; i < 4; i++)
3707     *q-- = *p++;
3708 #endif
3709 #ifdef IEEE
3710   if (REAL_WORDS_BIG_ENDIAN)
3711     {
3712       for (i = 0; i < 4; i++)
3713         *q++ = *p++;
3714     }
3715   else
3716     {
3717 #ifdef INFINITY
3718       if (eiisinf (a))
3719         {
3720           /* Intel long double infinity significand.  */
3721           *q-- = 0x8000;
3722           *q-- = 0;
3723           *q-- = 0;
3724           *q = 0;
3725           return;
3726         }
3727 #endif
3728       for (i = 0; i < 4; i++)
3729         *q-- = *p++;
3730     }
3731 #endif
3732 }
3733
3734 /* e type to double precision.  */
3735
3736 #ifdef DEC
3737 /* Convert e-type X to DEC-format double E.  */
3738
3739 static void
3740 etoe53 (x, e)
3741      unsigned EMUSHORT *x, *e;
3742 {
3743   etodec (x, e);                /* see etodec.c */
3744 }
3745
3746 /* Convert exploded e-type X, that has already been rounded to
3747    56-bit double precision, to DEC double Y.  */
3748
3749 static void
3750 toe53 (x, y)
3751      unsigned EMUSHORT *x, *y;
3752 {
3753   todec (x, y);
3754 }
3755
3756 #else
3757 #ifdef IBM
3758 /* Convert e-type X to IBM 370-format double E.  */
3759
3760 static void
3761 etoe53 (x, e)
3762      unsigned EMUSHORT *x, *e;
3763 {
3764   etoibm (x, e, DFmode);
3765 }
3766
3767 /* Convert exploded e-type X, that has already been rounded to
3768    56-bit precision, to IBM 370 double Y.  */
3769
3770 static void
3771 toe53 (x, y)
3772      unsigned EMUSHORT *x, *y;
3773 {
3774   toibm (x, y, DFmode);
3775 }
3776
3777 #else /* it's neither DEC nor IBM */
3778 #ifdef C4X
3779 /* Convert e-type X to C4X-format long double E.  */
3780
3781 static void
3782 etoe53 (x, e)
3783      unsigned EMUSHORT *x, *e;
3784 {
3785   etoc4x (x, e, HFmode);
3786 }
3787
3788 /* Convert exploded e-type X, that has already been rounded to
3789    56-bit precision, to IBM 370 double Y.  */
3790
3791 static void
3792 toe53 (x, y)
3793      unsigned EMUSHORT *x, *y;
3794 {
3795   toc4x (x, y, HFmode);
3796 }
3797
3798 #else  /* it's neither DEC nor IBM nor C4X */
3799
3800 /* Convert e-type X to IEEE double E.  */
3801
3802 static void
3803 etoe53 (x, e)
3804      unsigned EMUSHORT *x, *e;
3805 {
3806   unsigned EMUSHORT xi[NI];
3807   EMULONG exp;
3808   int rndsav;
3809
3810 #ifdef NANS
3811   if (eisnan (x))
3812     {
3813       make_nan (e, eisneg (x), DFmode);
3814       return;
3815     }
3816 #endif
3817   emovi (x, xi);
3818   /* adjust exponent for offsets */
3819   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3820 #ifdef INFINITY
3821   if (eisinf (x))
3822     goto nonorm;
3823 #endif
3824   /* round off to nearest or even */
3825   rndsav = rndprc;
3826   rndprc = 53;
3827   emdnorm (xi, 0, 0, exp, 64);
3828   rndprc = rndsav;
3829 #ifdef INFINITY
3830  nonorm:
3831 #endif
3832   toe53 (xi, e);
3833 }
3834
3835 /* Convert exploded e-type X, that has already been rounded to
3836    53-bit precision, to IEEE double Y.  */
3837
3838 static void
3839 toe53 (x, y)
3840      unsigned EMUSHORT *x, *y;
3841 {
3842   unsigned EMUSHORT i;
3843   unsigned EMUSHORT *p;
3844
3845 #ifdef NANS
3846   if (eiisnan (x))
3847     {
3848       make_nan (y, eiisneg (x), DFmode);
3849       return;
3850     }
3851 #endif
3852   p = &x[0];
3853 #ifdef IEEE
3854   if (! REAL_WORDS_BIG_ENDIAN)
3855     y += 3;
3856 #endif
3857   *y = 0;                       /* output high order */
3858   if (*p++)
3859     *y = 0x8000;                /* output sign bit */
3860
3861   i = *p++;
3862   if (i >= (unsigned int) 2047)
3863     {
3864       /* Saturate at largest number less than infinity.  */
3865 #ifdef INFINITY
3866       *y |= 0x7ff0;
3867       if (! REAL_WORDS_BIG_ENDIAN)
3868         {
3869           *(--y) = 0;
3870           *(--y) = 0;
3871           *(--y) = 0;
3872         }
3873       else
3874         {
3875           ++y;
3876           *y++ = 0;
3877           *y++ = 0;
3878           *y++ = 0;
3879         }
3880 #else
3881       *y |= (unsigned EMUSHORT) 0x7fef;
3882       if (! REAL_WORDS_BIG_ENDIAN)
3883         {
3884           *(--y) = 0xffff;
3885           *(--y) = 0xffff;
3886           *(--y) = 0xffff;
3887         }
3888       else
3889         {
3890           ++y;
3891           *y++ = 0xffff;
3892           *y++ = 0xffff;
3893           *y++ = 0xffff;
3894         }
3895 #endif
3896       return;
3897     }
3898   if (i == 0)
3899     {
3900       eshift (x, 4);
3901     }
3902   else
3903     {
3904       i <<= 4;
3905       eshift (x, 5);
3906     }
3907   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3908   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
3909   if (! REAL_WORDS_BIG_ENDIAN)
3910     {
3911       *(--y) = *p++;
3912       *(--y) = *p++;
3913       *(--y) = *p;
3914     }
3915   else
3916     {
3917       ++y;
3918       *y++ = *p++;
3919       *y++ = *p++;
3920       *y++ = *p++;
3921     }
3922 }
3923
3924 #endif /* not C4X */
3925 #endif /* not IBM */
3926 #endif /* not DEC */
3927
3928
3929
3930 /* e type to single precision.  */
3931
3932 #ifdef IBM
3933 /* Convert e-type X to IBM 370 float E.  */
3934
3935 static void
3936 etoe24 (x, e)
3937      unsigned EMUSHORT *x, *e;
3938 {
3939   etoibm (x, e, SFmode);
3940 }
3941
3942 /* Convert exploded e-type X, that has already been rounded to
3943    float precision, to IBM 370 float Y.  */
3944
3945 static void
3946 toe24 (x, y)
3947      unsigned EMUSHORT *x, *y;
3948 {
3949   toibm (x, y, SFmode);
3950 }
3951
3952 #else
3953
3954 #ifdef C4X
3955 /* Convert e-type X to C4X float E.  */
3956
3957 static void
3958 etoe24 (x, e)
3959      unsigned EMUSHORT *x, *e;
3960 {
3961   etoc4x (x, e, QFmode);
3962 }
3963
3964 /* Convert exploded e-type X, that has already been rounded to
3965    float precision, to IBM 370 float Y.  */
3966
3967 static void
3968 toe24 (x, y)
3969      unsigned EMUSHORT *x, *y;
3970 {
3971   toc4x (x, y, QFmode);
3972 }
3973
3974 #else
3975
3976 /* Convert e-type X to IEEE float E.  DEC float is the same as IEEE float.  */
3977
3978 static void
3979 etoe24 (x, e)
3980      unsigned EMUSHORT *x, *e;
3981 {
3982   EMULONG exp;
3983   unsigned EMUSHORT xi[NI];
3984   int rndsav;
3985
3986 #ifdef NANS
3987   if (eisnan (x))
3988     {
3989       make_nan (e, eisneg (x), SFmode);
3990       return;
3991     }
3992 #endif
3993   emovi (x, xi);
3994   /* adjust exponent for offsets */
3995   exp = (EMULONG) xi[E] - (EXONE - 0177);
3996 #ifdef INFINITY
3997   if (eisinf (x))
3998     goto nonorm;
3999 #endif
4000   /* round off to nearest or even */
4001   rndsav = rndprc;
4002   rndprc = 24;
4003   emdnorm (xi, 0, 0, exp, 64);
4004   rndprc = rndsav;
4005 #ifdef INFINITY
4006  nonorm:
4007 #endif
4008   toe24 (xi, e);
4009 }
4010
4011 /* Convert exploded e-type X, that has already been rounded to
4012    float precision, to IEEE float Y.  */
4013
4014 static void
4015 toe24 (x, y)
4016      unsigned EMUSHORT *x, *y;
4017 {
4018   unsigned EMUSHORT i;
4019   unsigned EMUSHORT *p;
4020
4021 #ifdef NANS
4022   if (eiisnan (x))
4023     {
4024       make_nan (y, eiisneg (x), SFmode);
4025       return;
4026     }
4027 #endif
4028   p = &x[0];
4029 #ifdef IEEE
4030   if (! REAL_WORDS_BIG_ENDIAN)
4031     y += 1;
4032 #endif
4033 #ifdef DEC
4034   y += 1;
4035 #endif
4036   *y = 0;                       /* output high order */
4037   if (*p++)
4038     *y = 0x8000;                /* output sign bit */
4039
4040   i = *p++;
4041 /* Handle overflow cases.  */
4042   if (i >= 255)
4043     {
4044 #ifdef INFINITY
4045       *y |= (unsigned EMUSHORT) 0x7f80;
4046 #ifdef DEC
4047       *(--y) = 0;
4048 #endif
4049 #ifdef IEEE
4050       if (! REAL_WORDS_BIG_ENDIAN)
4051         *(--y) = 0;
4052       else
4053         {
4054           ++y;
4055           *y = 0;
4056         }
4057 #endif
4058 #else  /* no INFINITY */
4059       *y |= (unsigned EMUSHORT) 0x7f7f;
4060 #ifdef DEC
4061       *(--y) = 0xffff;
4062 #endif
4063 #ifdef IEEE
4064       if (! REAL_WORDS_BIG_ENDIAN)
4065         *(--y) = 0xffff;
4066       else
4067         {
4068           ++y;
4069           *y = 0xffff;
4070         }
4071 #endif
4072 #ifdef ERANGE
4073       errno = ERANGE;
4074 #endif
4075 #endif  /* no INFINITY */
4076       return;
4077     }
4078   if (i == 0)
4079     {
4080       eshift (x, 7);
4081     }
4082   else
4083     {
4084       i <<= 7;
4085       eshift (x, 8);
4086     }
4087   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
4088   /* High order output already has sign bit set.  */
4089   *y |= i;
4090 #ifdef DEC
4091   *(--y) = *p;
4092 #endif
4093 #ifdef IEEE
4094   if (! REAL_WORDS_BIG_ENDIAN)
4095     *(--y) = *p;
4096   else
4097     {
4098       ++y;
4099       *y = *p;
4100     }
4101 #endif
4102 }
4103 #endif  /* not C4X */
4104 #endif  /* not IBM */
4105
4106 /* Compare two e type numbers.
4107    Return +1 if a > b
4108            0 if a == b
4109           -1 if a < b
4110           -2 if either a or b is a NaN.  */
4111
4112 static int
4113 ecmp (a, b)
4114      unsigned EMUSHORT *a, *b;
4115 {
4116   unsigned EMUSHORT ai[NI], bi[NI];
4117   register unsigned EMUSHORT *p, *q;
4118   register int i;
4119   int msign;
4120
4121 #ifdef NANS
4122   if (eisnan (a)  || eisnan (b))
4123       return (-2);
4124 #endif
4125   emovi (a, ai);
4126   p = ai;
4127   emovi (b, bi);
4128   q = bi;
4129
4130   if (*p != *q)
4131     {                           /* the signs are different */
4132       /* -0 equals + 0 */
4133       for (i = 1; i < NI - 1; i++)
4134         {
4135           if (ai[i] != 0)
4136             goto nzro;
4137           if (bi[i] != 0)
4138             goto nzro;
4139         }
4140       return (0);
4141     nzro:
4142       if (*p == 0)
4143         return (1);
4144       else
4145         return (-1);
4146     }
4147   /* both are the same sign */
4148   if (*p == 0)
4149     msign = 1;
4150   else
4151     msign = -1;
4152   i = NI - 1;
4153   do
4154     {
4155       if (*p++ != *q++)
4156         {
4157           goto diff;
4158         }
4159     }
4160   while (--i > 0);
4161
4162   return (0);                   /* equality */
4163
4164  diff:
4165
4166   if (*(--p) > *(--q))
4167     return (msign);             /* p is bigger */
4168   else
4169     return (-msign);            /* p is littler */
4170 }
4171
4172 #if 0
4173 /* Find e-type nearest integer to X, as floor (X + 0.5).  */
4174
4175 static void
4176 eround (x, y)
4177      unsigned EMUSHORT *x, *y;
4178 {
4179   eadd (ehalf, x, y);
4180   efloor (y, y);
4181 }
4182 #endif /* 0 */
4183
4184 /* Convert HOST_WIDE_INT LP to e type Y.  */
4185
4186 static void
4187 ltoe (lp, y)
4188      HOST_WIDE_INT *lp;
4189      unsigned EMUSHORT *y;
4190 {
4191   unsigned EMUSHORT yi[NI];
4192   unsigned HOST_WIDE_INT ll;
4193   int k;
4194
4195   ecleaz (yi);
4196   if (*lp < 0)
4197     {
4198       /* make it positive */
4199       ll = (unsigned HOST_WIDE_INT) (-(*lp));
4200       yi[0] = 0xffff;           /* put correct sign in the e type number */
4201     }
4202   else
4203     {
4204       ll = (unsigned HOST_WIDE_INT) (*lp);
4205     }
4206   /* move the long integer to yi significand area */
4207 #if HOST_BITS_PER_WIDE_INT == 64
4208   yi[M] = (unsigned EMUSHORT) (ll >> 48);
4209   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
4210   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
4211   yi[M + 3] = (unsigned EMUSHORT) ll;
4212   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
4213 #else
4214   yi[M] = (unsigned EMUSHORT) (ll >> 16);
4215   yi[M + 1] = (unsigned EMUSHORT) ll;
4216   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4217 #endif
4218
4219   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4220     ecleaz (yi);                /* it was zero */
4221   else
4222     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
4223   emovo (yi, y);                /* output the answer */
4224 }
4225
4226 /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
4227
4228 static void
4229 ultoe (lp, y)
4230      unsigned HOST_WIDE_INT *lp;
4231      unsigned EMUSHORT *y;
4232 {
4233   unsigned EMUSHORT yi[NI];
4234   unsigned HOST_WIDE_INT ll;
4235   int k;
4236
4237   ecleaz (yi);
4238   ll = *lp;
4239
4240   /* move the long integer to ayi significand area */
4241 #if HOST_BITS_PER_WIDE_INT == 64
4242   yi[M] = (unsigned EMUSHORT) (ll >> 48);
4243   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
4244   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
4245   yi[M + 3] = (unsigned EMUSHORT) ll;
4246   yi[E] = EXONE + 47;           /* exponent if normalize shift count were 0 */
4247 #else
4248   yi[M] = (unsigned EMUSHORT) (ll >> 16);
4249   yi[M + 1] = (unsigned EMUSHORT) ll;
4250   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
4251 #endif
4252
4253   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4254     ecleaz (yi);                /* it was zero */
4255   else
4256     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
4257   emovo (yi, y);                /* output the answer */
4258 }
4259
4260
4261 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4262    part FRAC of e-type (packed internal format) floating point input X.
4263    The integer output I has the sign of the input, except that
4264    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4265    The output e-type fraction FRAC is the positive fractional
4266    part of abs (X).  */
4267
4268 static void
4269 eifrac (x, i, frac)
4270      unsigned EMUSHORT *x;
4271      HOST_WIDE_INT *i;
4272      unsigned EMUSHORT *frac;
4273 {
4274   unsigned EMUSHORT xi[NI];
4275   int j, k;
4276   unsigned HOST_WIDE_INT ll;
4277
4278   emovi (x, xi);
4279   k = (int) xi[E] - (EXONE - 1);
4280   if (k <= 0)
4281     {
4282       /* if exponent <= 0, integer = 0 and real output is fraction */
4283       *i = 0L;
4284       emovo (xi, frac);
4285       return;
4286     }
4287   if (k > (HOST_BITS_PER_WIDE_INT - 1))
4288     {
4289       /* long integer overflow: output large integer
4290          and correct fraction  */
4291       if (xi[0])
4292         *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4293       else
4294         {
4295 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4296           /* In this case, let it overflow and convert as if unsigned.  */
4297           euifrac (x, &ll, frac);
4298           *i = (HOST_WIDE_INT) ll;
4299           return;
4300 #else
4301           /* In other cases, return the largest positive integer.  */
4302           *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4303 #endif
4304         }
4305       eshift (xi, k);
4306       if (extra_warnings)
4307         warning ("overflow on truncation to integer");
4308     }
4309   else if (k > 16)
4310     {
4311       /* Shift more than 16 bits: first shift up k-16 mod 16,
4312          then shift up by 16's.  */
4313       j = k - ((k >> 4) << 4);
4314       eshift (xi, j);
4315       ll = xi[M];
4316       k -= j;
4317       do
4318         {
4319           eshup6 (xi);
4320           ll = (ll << 16) | xi[M];
4321         }
4322       while ((k -= 16) > 0);
4323       *i = ll;
4324       if (xi[0])
4325         *i = -(*i);
4326     }
4327   else
4328       {
4329         /* shift not more than 16 bits */
4330           eshift (xi, k);
4331         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4332         if (xi[0])
4333           *i = -(*i);
4334       }
4335   xi[0] = 0;
4336   xi[E] = EXONE - 1;
4337   xi[M] = 0;
4338   if ((k = enormlz (xi)) > NBITS)
4339     ecleaz (xi);
4340   else
4341     xi[E] -= (unsigned EMUSHORT) k;
4342
4343   emovo (xi, frac);
4344 }
4345
4346
4347 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4348    FRAC of e-type X.  A negative input yields integer output = 0 but
4349    correct fraction.  */
4350
4351 static void
4352 euifrac (x, i, frac)
4353      unsigned EMUSHORT *x;
4354      unsigned HOST_WIDE_INT *i;
4355      unsigned EMUSHORT *frac;
4356 {
4357   unsigned HOST_WIDE_INT ll;
4358   unsigned EMUSHORT xi[NI];
4359   int j, k;
4360
4361   emovi (x, xi);
4362   k = (int) xi[E] - (EXONE - 1);
4363   if (k <= 0)
4364     {
4365       /* if exponent <= 0, integer = 0 and argument is fraction */
4366       *i = 0L;
4367       emovo (xi, frac);
4368       return;
4369     }
4370   if (k > HOST_BITS_PER_WIDE_INT)
4371     {
4372       /* Long integer overflow: output large integer
4373          and correct fraction.
4374          Note, the BSD microvax compiler says that ~(0UL)
4375          is a syntax error.  */
4376       *i = ~(0L);
4377       eshift (xi, k);
4378       if (extra_warnings)
4379         warning ("overflow on truncation to unsigned integer");
4380     }
4381   else if (k > 16)
4382     {
4383       /* Shift more than 16 bits: first shift up k-16 mod 16,
4384          then shift up by 16's.  */
4385       j = k - ((k >> 4) << 4);
4386       eshift (xi, j);
4387       ll = xi[M];
4388       k -= j;
4389       do
4390         {
4391           eshup6 (xi);
4392           ll = (ll << 16) | xi[M];
4393         }
4394       while ((k -= 16) > 0);
4395       *i = ll;
4396     }
4397   else
4398     {
4399       /* shift not more than 16 bits */
4400       eshift (xi, k);
4401       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4402     }
4403
4404   if (xi[0])  /* A negative value yields unsigned integer 0.  */
4405     *i = 0L;
4406
4407   xi[0] = 0;
4408   xi[E] = EXONE - 1;
4409   xi[M] = 0;
4410   if ((k = enormlz (xi)) > NBITS)
4411     ecleaz (xi);
4412   else
4413     xi[E] -= (unsigned EMUSHORT) k;
4414
4415   emovo (xi, frac);
4416 }
4417
4418 /* Shift the significand of exploded e-type X up or down by SC bits.  */
4419
4420 static int
4421 eshift (x, sc)
4422      unsigned EMUSHORT *x;
4423      int sc;
4424 {
4425   unsigned EMUSHORT lost;
4426   unsigned EMUSHORT *p;
4427
4428   if (sc == 0)
4429     return (0);
4430
4431   lost = 0;
4432   p = x + NI - 1;
4433
4434   if (sc < 0)
4435     {
4436       sc = -sc;
4437       while (sc >= 16)
4438         {
4439           lost |= *p;           /* remember lost bits */
4440           eshdn6 (x);
4441           sc -= 16;
4442         }
4443
4444       while (sc >= 8)
4445         {
4446           lost |= *p & 0xff;
4447           eshdn8 (x);
4448           sc -= 8;
4449         }
4450
4451       while (sc > 0)
4452         {
4453           lost |= *p & 1;
4454           eshdn1 (x);
4455           sc -= 1;
4456         }
4457     }
4458   else
4459     {
4460       while (sc >= 16)
4461         {
4462           eshup6 (x);
4463           sc -= 16;
4464         }
4465
4466       while (sc >= 8)
4467         {
4468           eshup8 (x);
4469           sc -= 8;
4470         }
4471
4472       while (sc > 0)
4473         {
4474           eshup1 (x);
4475           sc -= 1;
4476         }
4477     }
4478   if (lost)
4479     lost = 1;
4480   return ((int) lost);
4481 }
4482
4483 /* Shift normalize the significand area of exploded e-type X.
4484    Return the shift count (up = positive).  */
4485
4486 static int
4487 enormlz (x)
4488      unsigned EMUSHORT x[];
4489 {
4490   register unsigned EMUSHORT *p;
4491   int sc;
4492
4493   sc = 0;
4494   p = &x[M];
4495   if (*p != 0)
4496     goto normdn;
4497   ++p;
4498   if (*p & 0x8000)
4499     return (0);                 /* already normalized */
4500   while (*p == 0)
4501     {
4502       eshup6 (x);
4503       sc += 16;
4504
4505       /* With guard word, there are NBITS+16 bits available.
4506        Return true if all are zero.  */
4507       if (sc > NBITS)
4508         return (sc);
4509     }
4510   /* see if high byte is zero */
4511   while ((*p & 0xff00) == 0)
4512     {
4513       eshup8 (x);
4514       sc += 8;
4515     }
4516   /* now shift 1 bit at a time */
4517   while ((*p & 0x8000) == 0)
4518     {
4519       eshup1 (x);
4520       sc += 1;
4521       if (sc > NBITS)
4522         {
4523           mtherr ("enormlz", UNDERFLOW);
4524           return (sc);
4525         }
4526     }
4527   return (sc);
4528
4529   /* Normalize by shifting down out of the high guard word
4530      of the significand */
4531  normdn:
4532
4533   if (*p & 0xff00)
4534     {
4535       eshdn8 (x);
4536       sc -= 8;
4537     }
4538   while (*p != 0)
4539     {
4540       eshdn1 (x);
4541       sc -= 1;
4542
4543       if (sc < -NBITS)
4544         {
4545           mtherr ("enormlz", OVERFLOW);
4546           return (sc);
4547         }
4548     }
4549   return (sc);
4550 }
4551
4552 /* Powers of ten used in decimal <-> binary conversions.  */
4553
4554 #define NTEN 12
4555 #define MAXP 4096
4556
4557 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128
4558 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4559 {
4560   {0x6576, 0x4a92, 0x804a, 0x153f,
4561    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4562   {0x6a32, 0xce52, 0x329a, 0x28ce,
4563    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4564   {0x526c, 0x50ce, 0xf18b, 0x3d28,
4565    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4566   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4567    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4568   {0x851e, 0xeab7, 0x98fe, 0x901b,
4569    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4570   {0x0235, 0x0137, 0x36b1, 0x336c,
4571    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4572   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4573    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4574   {0x0000, 0x0000, 0x0000, 0x0000,
4575    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4576   {0x0000, 0x0000, 0x0000, 0x0000,
4577    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4578   {0x0000, 0x0000, 0x0000, 0x0000,
4579    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4580   {0x0000, 0x0000, 0x0000, 0x0000,
4581    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4582   {0x0000, 0x0000, 0x0000, 0x0000,
4583    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4584   {0x0000, 0x0000, 0x0000, 0x0000,
4585    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4586 };
4587
4588 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4589 {
4590   {0x2030, 0xcffc, 0xa1c3, 0x8123,
4591    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4592   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4593    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4594   {0xf53f, 0xf698, 0x6bd3, 0x0158,
4595    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4596   {0xe731, 0x04d4, 0xe3f2, 0xd332,
4597    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4598   {0xa23e, 0x5308, 0xfefb, 0x1155,
4599    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4600   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4601    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4602   {0x2a20, 0x6224, 0x47b3, 0x98d7,
4603    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4604   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4605    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4606   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4607    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4608   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4609    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4610   {0xc155, 0xa4a8, 0x404e, 0x6113,
4611    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4612   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4613    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4614   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4615    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4616 };
4617 #else
4618 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4619 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4620 {
4621   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
4622   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
4623   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4624   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4625   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4626   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4627   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4628   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4629   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4630   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4631   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4632   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4633   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
4634 };
4635
4636 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4637 {
4638   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
4639   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
4640   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4641   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4642   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4643   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4644   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4645   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4646   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4647   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4648   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4649   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4650   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
4651 };
4652 #endif
4653
4654 #if 0
4655 /* Convert float value X to ASCII string STRING with NDIG digits after
4656    the decimal point.  */
4657
4658 static void
4659 e24toasc (x, string, ndigs)
4660      unsigned EMUSHORT x[];
4661      char *string;
4662      int ndigs;
4663 {
4664   unsigned EMUSHORT w[NI];
4665
4666   e24toe (x, w);
4667   etoasc (w, string, ndigs);
4668 }
4669
4670 /* Convert double value X to ASCII string STRING with NDIG digits after
4671    the decimal point.  */
4672
4673 static void
4674 e53toasc (x, string, ndigs)
4675      unsigned EMUSHORT x[];
4676      char *string;
4677      int ndigs;
4678 {
4679   unsigned EMUSHORT w[NI];
4680
4681   e53toe (x, w);
4682   etoasc (w, string, ndigs);
4683 }
4684
4685 /* Convert double extended value X to ASCII string STRING with NDIG digits
4686    after the decimal point.  */
4687
4688 static void
4689 e64toasc (x, string, ndigs)
4690      unsigned EMUSHORT x[];
4691      char *string;
4692      int ndigs;
4693 {
4694   unsigned EMUSHORT w[NI];
4695
4696   e64toe (x, w);
4697   etoasc (w, string, ndigs);
4698 }
4699
4700 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4701    after the decimal point.  */
4702
4703 static void
4704 e113toasc (x, string, ndigs)
4705      unsigned EMUSHORT x[];
4706      char *string;
4707      int ndigs;
4708 {
4709   unsigned EMUSHORT w[NI];
4710
4711   e113toe (x, w);
4712   etoasc (w, string, ndigs);
4713 }
4714 #endif /* 0 */
4715
4716 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4717    the decimal point.  */
4718
4719 static char wstring[80];        /* working storage for ASCII output */
4720
4721 static void
4722 etoasc (x, string, ndigs)
4723      unsigned EMUSHORT x[];
4724      char *string;
4725      int ndigs;
4726 {
4727   EMUSHORT digit;
4728   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4729   unsigned EMUSHORT *p, *r, *ten;
4730   unsigned EMUSHORT sign;
4731   int i, j, k, expon, rndsav;
4732   char *s, *ss;
4733   unsigned EMUSHORT m;
4734
4735
4736   rndsav = rndprc;
4737   ss = string;
4738   s = wstring;
4739   *ss = '\0';
4740   *s = '\0';
4741 #ifdef NANS
4742   if (eisnan (x))
4743     {
4744       sprintf (wstring, " NaN ");
4745       goto bxit;
4746     }
4747 #endif
4748   rndprc = NBITS;               /* set to full precision */
4749   emov (x, y);                  /* retain external format */
4750   if (y[NE - 1] & 0x8000)
4751     {
4752       sign = 0xffff;
4753       y[NE - 1] &= 0x7fff;
4754     }
4755   else
4756     {
4757       sign = 0;
4758     }
4759   expon = 0;
4760   ten = &etens[NTEN][0];
4761   emov (eone, t);
4762   /* Test for zero exponent */
4763   if (y[NE - 1] == 0)
4764     {
4765       for (k = 0; k < NE - 1; k++)
4766         {
4767           if (y[k] != 0)
4768             goto tnzro;         /* denormalized number */
4769         }
4770       goto isone;               /* valid all zeros */
4771     }
4772  tnzro:
4773
4774   /* Test for infinity.  */
4775   if (y[NE - 1] == 0x7fff)
4776     {
4777       if (sign)
4778         sprintf (wstring, " -Infinity ");
4779       else
4780         sprintf (wstring, " Infinity ");
4781       goto bxit;
4782     }
4783
4784   /* Test for exponent nonzero but significand denormalized.
4785    * This is an error condition.
4786    */
4787   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4788     {
4789       mtherr ("etoasc", DOMAIN);
4790       sprintf (wstring, "NaN");
4791       goto bxit;
4792     }
4793
4794   /* Compare to 1.0 */
4795   i = ecmp (eone, y);
4796   if (i == 0)
4797     goto isone;
4798
4799   if (i == -2)
4800     abort ();
4801
4802   if (i < 0)
4803     {                           /* Number is greater than 1 */
4804       /* Convert significand to an integer and strip trailing decimal zeros.  */
4805       emov (y, u);
4806       u[NE - 1] = EXONE + NBITS - 1;
4807
4808       p = &etens[NTEN - 4][0];
4809       m = 16;
4810       do
4811         {
4812           ediv (p, u, t);
4813           efloor (t, w);
4814           for (j = 0; j < NE - 1; j++)
4815             {
4816               if (t[j] != w[j])
4817                 goto noint;
4818             }
4819           emov (t, u);
4820           expon += (int) m;
4821         noint:
4822           p += NE;
4823           m >>= 1;
4824         }
4825       while (m != 0);
4826
4827       /* Rescale from integer significand */
4828       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4829       emov (u, y);
4830       /* Find power of 10 */
4831       emov (eone, t);
4832       m = MAXP;
4833       p = &etens[0][0];
4834       /* An unordered compare result shouldn't happen here.  */
4835       while (ecmp (ten, u) <= 0)
4836         {
4837           if (ecmp (p, u) <= 0)
4838             {
4839               ediv (p, u, u);
4840               emul (p, t, t);
4841               expon += (int) m;
4842             }
4843           m >>= 1;
4844           if (m == 0)
4845             break;
4846           p += NE;
4847         }
4848     }
4849   else
4850     {                           /* Number is less than 1.0 */
4851       /* Pad significand with trailing decimal zeros.  */
4852       if (y[NE - 1] == 0)
4853         {
4854           while ((y[NE - 2] & 0x8000) == 0)
4855             {
4856               emul (ten, y, y);
4857               expon -= 1;
4858             }
4859         }
4860       else
4861         {
4862           emovi (y, w);
4863           for (i = 0; i < NDEC + 1; i++)
4864             {
4865               if ((w[NI - 1] & 0x7) != 0)
4866                 break;
4867               /* multiply by 10 */
4868               emovz (w, u);
4869               eshdn1 (u);
4870               eshdn1 (u);
4871               eaddm (w, u);
4872               u[1] += 3;
4873               while (u[2] != 0)
4874                 {
4875                   eshdn1 (u);
4876                   u[1] += 1;
4877                 }
4878               if (u[NI - 1] != 0)
4879                 break;
4880               if (eone[NE - 1] <= u[1])
4881                 break;
4882               emovz (u, w);
4883               expon -= 1;
4884             }
4885           emovo (w, y);
4886         }
4887       k = -MAXP;
4888       p = &emtens[0][0];
4889       r = &etens[0][0];
4890       emov (y, w);
4891       emov (eone, t);
4892       while (ecmp (eone, w) > 0)
4893         {
4894           if (ecmp (p, w) >= 0)
4895             {
4896               emul (r, w, w);
4897               emul (r, t, t);
4898               expon += k;
4899             }
4900           k /= 2;
4901           if (k == 0)
4902             break;
4903           p += NE;
4904           r += NE;
4905         }
4906       ediv (t, eone, t);
4907     }
4908  isone:
4909   /* Find the first (leading) digit.  */
4910   emovi (t, w);
4911   emovz (w, t);
4912   emovi (y, w);
4913   emovz (w, y);
4914   eiremain (t, y);
4915   digit = equot[NI - 1];
4916   while ((digit == 0) && (ecmp (y, ezero) != 0))
4917     {
4918       eshup1 (y);
4919       emovz (y, u);
4920       eshup1 (u);
4921       eshup1 (u);
4922       eaddm (u, y);
4923       eiremain (t, y);
4924       digit = equot[NI - 1];
4925       expon -= 1;
4926     }
4927   s = wstring;
4928   if (sign)
4929     *s++ = '-';
4930   else
4931     *s++ = ' ';
4932   /* Examine number of digits requested by caller.  */
4933   if (ndigs < 0)
4934     ndigs = 0;
4935   if (ndigs > NDEC)
4936     ndigs = NDEC;
4937   if (digit == 10)
4938     {
4939       *s++ = '1';
4940       *s++ = '.';
4941       if (ndigs > 0)
4942         {
4943           *s++ = '0';
4944           ndigs -= 1;
4945         }
4946       expon += 1;
4947     }
4948   else
4949     {
4950       *s++ = (char)digit + '0';
4951       *s++ = '.';
4952     }
4953   /* Generate digits after the decimal point.  */
4954   for (k = 0; k <= ndigs; k++)
4955     {
4956       /* multiply current number by 10, without normalizing */
4957       eshup1 (y);
4958       emovz (y, u);
4959       eshup1 (u);
4960       eshup1 (u);
4961       eaddm (u, y);
4962       eiremain (t, y);
4963       *s++ = (char) equot[NI - 1] + '0';
4964     }
4965   digit = equot[NI - 1];
4966   --s;
4967   ss = s;
4968   /* round off the ASCII string */
4969   if (digit > 4)
4970     {
4971       /* Test for critical rounding case in ASCII output.  */
4972       if (digit == 5)
4973         {
4974           emovo (y, t);
4975           if (ecmp (t, ezero) != 0)
4976             goto roun;          /* round to nearest */
4977 #ifndef C4X
4978           if ((*(s - 1) & 1) == 0)
4979             goto doexp;         /* round to even */
4980 #endif
4981         }
4982       /* Round up and propagate carry-outs */
4983     roun:
4984       --s;
4985       k = *s & 0x7f;
4986       /* Carry out to most significant digit? */
4987       if (k == '.')
4988         {
4989           --s;
4990           k = *s;
4991           k += 1;
4992           *s = (char) k;
4993           /* Most significant digit carries to 10? */
4994           if (k > '9')
4995             {
4996               expon += 1;
4997               *s = '1';
4998             }
4999           goto doexp;
5000         }
5001       /* Round up and carry out from less significant digits */
5002       k += 1;
5003       *s = (char) k;
5004       if (k > '9')
5005         {
5006           *s = '0';
5007           goto roun;
5008         }
5009     }
5010  doexp:
5011   /*
5012      if (expon >= 0)
5013      sprintf (ss, "e+%d", expon);
5014      else
5015      sprintf (ss, "e%d", expon);
5016      */
5017   sprintf (ss, "e%d", expon);
5018  bxit:
5019   rndprc = rndsav;
5020   /* copy out the working string */
5021   s = string;
5022   ss = wstring;
5023   while (*ss == ' ')            /* strip possible leading space */
5024     ++ss;
5025   while ((*s++ = *ss++) != '\0')
5026     ;
5027 }
5028
5029
5030 /* Convert ASCII string to floating point.
5031
5032    Numeric input is a free format decimal number of any length, with
5033    or without decimal point.  Entering E after the number followed by an
5034    integer number causes the second number to be interpreted as a power of
5035    10 to be multiplied by the first number (i.e., "scientific" notation).  */
5036
5037 /* Convert ASCII string S to single precision float value Y.  */
5038
5039 static void
5040 asctoe24 (s, y)
5041      const char *s;
5042      unsigned EMUSHORT *y;
5043 {
5044   asctoeg (s, y, 24);
5045 }
5046
5047
5048 /* Convert ASCII string S to double precision value Y.  */
5049
5050 static void
5051 asctoe53 (s, y)
5052      const char *s;
5053      unsigned EMUSHORT *y;
5054 {
5055 #if defined(DEC) || defined(IBM)
5056   asctoeg (s, y, 56);
5057 #else
5058 #if defined(C4X)
5059   asctoeg (s, y, 32);
5060 #else
5061   asctoeg (s, y, 53);
5062 #endif
5063 #endif
5064 }
5065
5066
5067 /* Convert ASCII string S to double extended value Y.  */
5068
5069 static void
5070 asctoe64 (s, y)
5071      const char *s;
5072      unsigned EMUSHORT *y;
5073 {
5074   asctoeg (s, y, 64);
5075 }
5076
5077 /* Convert ASCII string S to 128-bit long double Y.  */
5078
5079 static void
5080 asctoe113 (s, y)
5081      const char *s;
5082      unsigned EMUSHORT *y;
5083 {
5084   asctoeg (s, y, 113);
5085 }
5086
5087 /* Convert ASCII string S to e type Y.  */
5088
5089 static void
5090 asctoe (s, y)
5091      const char *s;
5092      unsigned EMUSHORT *y;
5093 {
5094   asctoeg (s, y, NBITS);
5095 }
5096
5097 /* Convert ASCII string SS to e type Y, with a specified rounding precision
5098    of OPREC bits.  BASE is 16 for C9X hexadecimal floating constants.  */
5099
5100 static void
5101 asctoeg (ss, y, oprec)
5102      const char *ss;
5103      unsigned EMUSHORT *y;
5104      int oprec;
5105 {
5106   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
5107   int esign, decflg, sgnflg, nexp, exp, prec, lost;
5108   int k, trail, c, rndsav;
5109   EMULONG lexp;
5110   unsigned EMUSHORT nsign, *p;
5111   char *sp, *s, *lstr;
5112   int base = 10;
5113
5114   /* Copy the input string.  */
5115   lstr = (char *) alloca (strlen (ss) + 1);
5116
5117   while (*ss == ' ')            /* skip leading spaces */
5118     ++ss;
5119
5120   sp = lstr;
5121   while ((*sp++ = *ss++) != '\0')
5122     ;
5123   s = lstr;
5124
5125   if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
5126     {
5127       base = 16;
5128       s += 2;
5129     }
5130
5131   rndsav = rndprc;
5132   rndprc = NBITS;               /* Set to full precision */
5133   lost = 0;
5134   nsign = 0;
5135   decflg = 0;
5136   sgnflg = 0;
5137   nexp = 0;
5138   exp = 0;
5139   prec = 0;
5140   ecleaz (yy);
5141   trail = 0;
5142
5143  nxtcom:
5144   if (*s >= '0' && *s <= '9')
5145     k = *s - '0';
5146   else if (*s >= 'a')
5147     k = 10 + *s - 'a';
5148   else
5149     k = 10 + *s - 'A';
5150   if ((k >= 0) && (k < base))
5151     {
5152       /* Ignore leading zeros */
5153       if ((prec == 0) && (decflg == 0) && (k == 0))
5154         goto donchr;
5155       /* Identify and strip trailing zeros after the decimal point.  */
5156       if ((trail == 0) && (decflg != 0))
5157         {
5158           sp = s;
5159           while ((*sp >= '0' && *sp <= '9')
5160                  || (base == 16 && ((*sp >= 'a' && *sp <= 'f')
5161                                     || (*sp >= 'A' && *sp <= 'F'))))
5162             ++sp;
5163           /* Check for syntax error */
5164           c = *sp & 0x7f;
5165           if ((base != 10 || ((c != 'e') && (c != 'E')))
5166               && (base != 16 || ((c != 'p') && (c != 'P')))
5167               && (c != '\0')
5168               && (c != '\n') && (c != '\r') && (c != ' ')
5169               && (c != ','))
5170             goto error;
5171           --sp;
5172           while (*sp == '0')
5173             *sp-- = 'z';
5174           trail = 1;
5175           if (*s == 'z')
5176             goto donchr;
5177         }
5178
5179       /* If enough digits were given to more than fill up the yy register,
5180          continuing until overflow into the high guard word yy[2]
5181          guarantees that there will be a roundoff bit at the top
5182          of the low guard word after normalization.  */
5183
5184       if (yy[2] == 0)
5185         {
5186           if (base == 16)
5187             {
5188               if (decflg)
5189                 nexp += 4;      /* count digits after decimal point */
5190
5191               eshup1 (yy);      /* multiply current number by 16 */
5192               eshup1 (yy);
5193               eshup1 (yy);
5194               eshup1 (yy);
5195             }
5196           else
5197             {
5198               if (decflg)
5199                 nexp += 1;              /* count digits after decimal point */
5200
5201               eshup1 (yy);              /* multiply current number by 10 */
5202               emovz (yy, xt);
5203               eshup1 (xt);
5204               eshup1 (xt);
5205               eaddm (xt, yy);
5206             }
5207           /* Insert the current digit.  */
5208           ecleaz (xt);
5209           xt[NI - 2] = (unsigned EMUSHORT) k;
5210           eaddm (xt, yy);
5211         }
5212       else
5213         {
5214           /* Mark any lost non-zero digit.  */
5215           lost |= k;
5216           /* Count lost digits before the decimal point.  */
5217           if (decflg == 0)
5218             {
5219               if (base == 10)
5220                 nexp -= 1;
5221               else
5222                 nexp -= 4;
5223             }
5224         }
5225       prec += 1;
5226       goto donchr;
5227     }
5228
5229   switch (*s)
5230     {
5231     case 'z':
5232       break;
5233     case 'E':
5234     case 'e':
5235     case 'P':
5236     case 'p':
5237       goto expnt;
5238     case '.':                   /* decimal point */
5239       if (decflg)
5240         goto error;
5241       ++decflg;
5242       break;
5243     case '-':
5244       nsign = 0xffff;
5245       if (sgnflg)
5246         goto error;
5247       ++sgnflg;
5248       break;
5249     case '+':
5250       if (sgnflg)
5251         goto error;
5252       ++sgnflg;
5253       break;
5254     case ',':
5255     case ' ':
5256     case '\0':
5257     case '\n':
5258     case '\r':
5259       goto daldone;
5260     case 'i':
5261     case 'I':
5262       goto infinite;
5263     default:
5264     error:
5265 #ifdef NANS
5266       einan (yy);
5267 #else
5268       mtherr ("asctoe", DOMAIN);
5269       eclear (yy);
5270 #endif
5271       goto aexit;
5272     }
5273  donchr:
5274   ++s;
5275   goto nxtcom;
5276
5277   /* Exponent interpretation */
5278  expnt:
5279   /* 0.0eXXX is zero, regardless of XXX.  Check for the 0.0. */
5280   for (k = 0; k < NI; k++)
5281     {
5282       if (yy[k] != 0)
5283         goto read_expnt;
5284     }
5285   goto aexit;
5286
5287 read_expnt:
5288   esign = 1;
5289   exp = 0;
5290   ++s;
5291   /* check for + or - */
5292   if (*s == '-')
5293     {
5294       esign = -1;
5295       ++s;
5296     }
5297   if (*s == '+')
5298     ++s;
5299   while ((*s >= '0') && (*s <= '9'))
5300     {
5301       exp *= 10;
5302       exp += *s++ - '0';
5303       if (exp > 999999)
5304         break;
5305     }
5306   if (esign < 0)
5307     exp = -exp;
5308   if ((exp > MAXDECEXP) && (base == 10))
5309     {
5310  infinite:
5311       ecleaz (yy);
5312       yy[E] = 0x7fff;           /* infinity */
5313       goto aexit;
5314     }
5315   if ((exp < MINDECEXP) && (base == 10))
5316     {
5317  zero:
5318       ecleaz (yy);
5319       goto aexit;
5320     }
5321
5322  daldone:
5323   if (base == 16)
5324     {
5325       /* Base 16 hexadecimal floating constant.  */
5326       if ((k = enormlz (yy)) > NBITS)
5327         {
5328           ecleaz (yy);
5329           goto aexit;
5330         }
5331       /* Adjust the exponent.  NEXP is the number of hex digits,
5332          EXP is a power of 2.  */
5333       lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5334       if (lexp > 0x7fff)
5335         goto infinite;
5336       if (lexp < 0)
5337         goto zero;
5338       yy[E] = lexp;
5339       goto expdon;
5340     }
5341
5342   nexp = exp - nexp;
5343   /* Pad trailing zeros to minimize power of 10, per IEEE spec.  */
5344   while ((nexp > 0) && (yy[2] == 0))
5345     {
5346       emovz (yy, xt);
5347       eshup1 (xt);
5348       eshup1 (xt);
5349       eaddm (yy, xt);
5350       eshup1 (xt);
5351       if (xt[2] != 0)
5352         break;
5353       nexp -= 1;
5354       emovz (xt, yy);
5355     }
5356   if ((k = enormlz (yy)) > NBITS)
5357     {
5358       ecleaz (yy);
5359       goto aexit;
5360     }
5361   lexp = (EXONE - 1 + NBITS) - k;
5362   emdnorm (yy, lost, 0, lexp, 64);
5363   lost = 0;
5364
5365   /* Convert to external format:
5366
5367      Multiply by 10**nexp.  If precision is 64 bits,
5368      the maximum relative error incurred in forming 10**n
5369      for 0 <= n <= 324 is 8.2e-20, at 10**180.
5370      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5371      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
5372
5373   lexp = yy[E];
5374   if (nexp == 0)
5375     {
5376       k = 0;
5377       goto expdon;
5378     }
5379   esign = 1;
5380   if (nexp < 0)
5381     {
5382       nexp = -nexp;
5383       esign = -1;
5384       if (nexp > 4096)
5385         {
5386           /* Punt.  Can't handle this without 2 divides.  */
5387           emovi (etens[0], tt);
5388           lexp -= tt[E];
5389           k = edivm (tt, yy);
5390           lexp += EXONE;
5391           nexp -= 4096;
5392         }
5393     }
5394   p = &etens[NTEN][0];
5395   emov (eone, xt);
5396   exp = 1;
5397   do
5398     {
5399       if (exp & nexp)
5400         emul (p, xt, xt);
5401       p -= NE;
5402       exp = exp + exp;
5403     }
5404   while (exp <= MAXP);
5405
5406   emovi (xt, tt);
5407   if (esign < 0)
5408     {
5409       lexp -= tt[E];
5410       k = edivm (tt, yy);
5411       lexp += EXONE;
5412     }
5413   else
5414     {
5415       lexp += tt[E];
5416       k = emulm (tt, yy);
5417       lexp -= EXONE - 1;
5418     }
5419   lost = k;
5420
5421  expdon:
5422
5423   /* Round and convert directly to the destination type */
5424   if (oprec == 53)
5425     lexp -= EXONE - 0x3ff;
5426 #ifdef C4X
5427   else if (oprec == 24 || oprec == 32)
5428     lexp -= (EXONE - 0x7f);
5429 #else
5430 #ifdef IBM
5431   else if (oprec == 24 || oprec == 56)
5432     lexp -= EXONE - (0x41 << 2);
5433 #else
5434   else if (oprec == 24)
5435     lexp -= EXONE - 0177;
5436 #endif /* IBM */
5437 #endif /* C4X */
5438 #ifdef DEC
5439   else if (oprec == 56)
5440     lexp -= EXONE - 0201;
5441 #endif
5442   rndprc = oprec;
5443   emdnorm (yy, lost, 0, lexp, 64);
5444
5445  aexit:
5446
5447   rndprc = rndsav;
5448   yy[0] = nsign;
5449   switch (oprec)
5450     {
5451 #ifdef DEC
5452     case 56:
5453       todec (yy, y);            /* see etodec.c */
5454       break;
5455 #endif
5456 #ifdef IBM
5457     case 56:
5458       toibm (yy, y, DFmode);
5459       break;
5460 #endif
5461 #ifdef C4X
5462     case 32:
5463       toc4x (yy, y, HFmode);
5464       break;
5465 #endif
5466
5467     case 53:
5468       toe53 (yy, y);
5469       break;
5470     case 24:
5471       toe24 (yy, y);
5472       break;
5473     case 64:
5474       toe64 (yy, y);
5475       break;
5476     case 113:
5477       toe113 (yy, y);
5478       break;
5479     case NBITS:
5480       emovo (yy, y);
5481       break;
5482     }
5483 }
5484
5485
5486
5487 /* Return Y = largest integer not greater than X (truncated toward minus
5488    infinity).  */
5489
5490 static unsigned EMUSHORT bmask[] =
5491 {
5492   0xffff,
5493   0xfffe,
5494   0xfffc,
5495   0xfff8,
5496   0xfff0,
5497   0xffe0,
5498   0xffc0,
5499   0xff80,
5500   0xff00,
5501   0xfe00,
5502   0xfc00,
5503   0xf800,
5504   0xf000,
5505   0xe000,
5506   0xc000,
5507   0x8000,
5508   0x0000,
5509 };
5510
5511 static void
5512 efloor (x, y)
5513      unsigned EMUSHORT x[], y[];
5514 {
5515   register unsigned EMUSHORT *p;
5516   int e, expon, i;
5517   unsigned EMUSHORT f[NE];
5518
5519   emov (x, f);                  /* leave in external format */
5520   expon = (int) f[NE - 1];
5521   e = (expon & 0x7fff) - (EXONE - 1);
5522   if (e <= 0)
5523     {
5524       eclear (y);
5525       goto isitneg;
5526     }
5527   /* number of bits to clear out */
5528   e = NBITS - e;
5529   emov (f, y);
5530   if (e <= 0)
5531     return;
5532
5533   p = &y[0];
5534   while (e >= 16)
5535     {
5536       *p++ = 0;
5537       e -= 16;
5538     }
5539   /* clear the remaining bits */
5540   *p &= bmask[e];
5541   /* truncate negatives toward minus infinity */
5542  isitneg:
5543
5544   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5545     {
5546       for (i = 0; i < NE - 1; i++)
5547         {
5548           if (f[i] != y[i])
5549             {
5550               esub (eone, y, y);
5551               break;
5552             }
5553         }
5554     }
5555 }
5556
5557
5558 #if 0
5559 /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
5560    For example, 1.1 = 0.55 * 2^1.  */
5561
5562 static void
5563 efrexp (x, exp, s)
5564      unsigned EMUSHORT x[];
5565      int *exp;
5566      unsigned EMUSHORT s[];
5567 {
5568   unsigned EMUSHORT xi[NI];
5569   EMULONG li;
5570
5571   emovi (x, xi);
5572   /*  Handle denormalized numbers properly using long integer exponent.  */
5573   li = (EMULONG) ((EMUSHORT) xi[1]);
5574
5575   if (li == 0)
5576     {
5577       li -= enormlz (xi);
5578     }
5579   xi[1] = 0x3ffe;
5580   emovo (xi, s);
5581   *exp = (int) (li - 0x3ffe);
5582 }
5583 #endif
5584
5585 /* Return e type Y = X * 2^PWR2.  */
5586
5587 static void
5588 eldexp (x, pwr2, y)
5589      unsigned EMUSHORT x[];
5590      int pwr2;
5591      unsigned EMUSHORT y[];
5592 {
5593   unsigned EMUSHORT xi[NI];
5594   EMULONG li;
5595   int i;
5596
5597   emovi (x, xi);
5598   li = xi[1];
5599   li += pwr2;
5600   i = 0;
5601   emdnorm (xi, i, i, li, 64);
5602   emovo (xi, y);
5603 }
5604
5605
5606 #if 0
5607 /* C = remainder after dividing B by A, all e type values.
5608    Least significant integer quotient bits left in EQUOT.  */
5609
5610 static void
5611 eremain (a, b, c)
5612      unsigned EMUSHORT a[], b[], c[];
5613 {
5614   unsigned EMUSHORT den[NI], num[NI];
5615
5616 #ifdef NANS
5617   if (eisinf (b)
5618       || (ecmp (a, ezero) == 0)
5619       || eisnan (a)
5620       || eisnan (b))
5621     {
5622       enan (c, 0);
5623       return;
5624     }
5625 #endif
5626   if (ecmp (a, ezero) == 0)
5627     {
5628       mtherr ("eremain", SING);
5629       eclear (c);
5630       return;
5631     }
5632   emovi (a, den);
5633   emovi (b, num);
5634   eiremain (den, num);
5635   /* Sign of remainder = sign of quotient */
5636   if (a[0] == b[0])
5637     num[0] = 0;
5638   else
5639     num[0] = 0xffff;
5640   emovo (num, c);
5641 }
5642 #endif
5643
5644 /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
5645     remainder in NUM.  */
5646
5647 static void
5648 eiremain (den, num)
5649      unsigned EMUSHORT den[], num[];
5650 {
5651   EMULONG ld, ln;
5652   unsigned EMUSHORT j;
5653
5654   ld = den[E];
5655   ld -= enormlz (den);
5656   ln = num[E];
5657   ln -= enormlz (num);
5658   ecleaz (equot);
5659   while (ln >= ld)
5660     {
5661       if (ecmpm (den, num) <= 0)
5662         {
5663           esubm (den, num);
5664           j = 1;
5665         }
5666       else
5667           j = 0;
5668       eshup1 (equot);
5669       equot[NI - 1] |= j;
5670       eshup1 (num);
5671       ln -= 1;
5672     }
5673   emdnorm (num, 0, 0, ln, 0);
5674 }
5675
5676 /* Report an error condition CODE encountered in function NAME.
5677
5678     Mnemonic        Value          Significance
5679
5680      DOMAIN            1       argument domain error
5681      SING              2       function singularity
5682      OVERFLOW          3       overflow range error
5683      UNDERFLOW         4       underflow range error
5684      TLOSS             5       total loss of precision
5685      PLOSS             6       partial loss of precision
5686      INVALID           7       NaN - producing operation
5687      EDOM             33       Unix domain error code
5688      ERANGE           34       Unix range error code
5689
5690    The order of appearance of the following messages is bound to the
5691    error codes defined above.  */
5692
5693 int merror = 0;
5694 extern int merror;
5695
5696 static void
5697 mtherr (name, code)
5698      const char *name;
5699      int code;
5700 {
5701   /* The string passed by the calling program is supposed to be the
5702      name of the function in which the error occurred.
5703      The code argument selects which error message string will be printed.  */
5704
5705   if (strcmp (name, "esub") == 0)
5706     name = "subtraction";
5707   else if (strcmp (name, "ediv") == 0)
5708     name = "division";
5709   else if (strcmp (name, "emul") == 0)
5710     name = "multiplication";
5711   else if (strcmp (name, "enormlz") == 0)
5712     name = "normalization";
5713   else if (strcmp (name, "etoasc") == 0)
5714     name = "conversion to text";
5715   else if (strcmp (name, "asctoe") == 0)
5716     name = "parsing";
5717   else if (strcmp (name, "eremain") == 0)
5718     name = "modulus";
5719   else if (strcmp (name, "esqrt") == 0)
5720     name = "square root";
5721   if (extra_warnings)
5722     {
5723       switch (code)
5724         {
5725         case DOMAIN:    warning ("%s: argument domain error"    , name); break;
5726         case SING:      warning ("%s: function singularity"     , name); break;
5727         case OVERFLOW:  warning ("%s: overflow range error"     , name); break;
5728         case UNDERFLOW: warning ("%s: underflow range error"    , name); break;
5729         case TLOSS:     warning ("%s: total loss of precision"  , name); break;
5730         case PLOSS:     warning ("%s: partial loss of precision", name); break;
5731         case INVALID:   warning ("%s: NaN - producing operation", name); break;
5732         default:        abort ();
5733         }
5734     }
5735
5736   /* Set global error message word */
5737   merror = code + 1;
5738 }
5739
5740 #ifdef DEC
5741 /* Convert DEC double precision D to e type E.  */
5742
5743 static void
5744 dectoe (d, e)
5745      unsigned EMUSHORT *d;
5746      unsigned EMUSHORT *e;
5747 {
5748   unsigned EMUSHORT y[NI];
5749   register unsigned EMUSHORT r, *p;
5750
5751   ecleaz (y);                   /* start with a zero */
5752   p = y;                        /* point to our number */
5753   r = *d;                       /* get DEC exponent word */
5754   if (*d & (unsigned int) 0x8000)
5755     *p = 0xffff;                /* fill in our sign */
5756   ++p;                          /* bump pointer to our exponent word */
5757   r &= 0x7fff;                  /* strip the sign bit */
5758   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
5759     goto done;
5760
5761
5762   r >>= 7;                      /* shift exponent word down 7 bits */
5763   r += EXONE - 0201;            /* subtract DEC exponent offset */
5764   /* add our e type exponent offset */
5765   *p++ = r;                     /* to form our exponent */
5766
5767   r = *d++;                     /* now do the high order mantissa */
5768   r &= 0177;                    /* strip off the DEC exponent and sign bits */
5769   r |= 0200;                    /* the DEC understood high order mantissa bit */
5770   *p++ = r;                     /* put result in our high guard word */
5771
5772   *p++ = *d++;                  /* fill in the rest of our mantissa */
5773   *p++ = *d++;
5774   *p = *d;
5775
5776   eshdn8 (y);                   /* shift our mantissa down 8 bits */
5777  done:
5778   emovo (y, e);
5779 }
5780
5781 /* Convert e type X to DEC double precision D.  */
5782
5783 static void
5784 etodec (x, d)
5785      unsigned EMUSHORT *x, *d;
5786 {
5787   unsigned EMUSHORT xi[NI];
5788   EMULONG exp;
5789   int rndsav;
5790
5791   emovi (x, xi);
5792   /* Adjust exponent for offsets.  */
5793   exp = (EMULONG) xi[E] - (EXONE - 0201);
5794   /* Round off to nearest or even.  */
5795   rndsav = rndprc;
5796   rndprc = 56;
5797   emdnorm (xi, 0, 0, exp, 64);
5798   rndprc = rndsav;
5799   todec (xi, d);
5800 }
5801
5802 /* Convert exploded e-type X, that has already been rounded to
5803    56-bit precision, to DEC format double Y.  */
5804
5805 static void
5806 todec (x, y)
5807      unsigned EMUSHORT *x, *y;
5808 {
5809   unsigned EMUSHORT i;
5810   unsigned EMUSHORT *p;
5811
5812   p = x;
5813   *y = 0;
5814   if (*p++)
5815     *y = 0100000;
5816   i = *p++;
5817   if (i == 0)
5818     {
5819       *y++ = 0;
5820       *y++ = 0;
5821       *y++ = 0;
5822       *y++ = 0;
5823       return;
5824     }
5825   if (i > 0377)
5826     {
5827       *y++ |= 077777;
5828       *y++ = 0xffff;
5829       *y++ = 0xffff;
5830       *y++ = 0xffff;
5831 #ifdef ERANGE
5832       errno = ERANGE;
5833 #endif
5834       return;
5835     }
5836   i &= 0377;
5837   i <<= 7;
5838   eshup8 (x);
5839   x[M] &= 0177;
5840   i |= x[M];
5841   *y++ |= i;
5842   *y++ = x[M + 1];
5843   *y++ = x[M + 2];
5844   *y++ = x[M + 3];
5845 }
5846 #endif /* DEC */
5847
5848 #ifdef IBM
5849 /* Convert IBM single/double precision to e type.  */
5850
5851 static void
5852 ibmtoe (d, e, mode)
5853      unsigned EMUSHORT *d;
5854      unsigned EMUSHORT *e;
5855      enum machine_mode mode;
5856 {
5857   unsigned EMUSHORT y[NI];
5858   register unsigned EMUSHORT r, *p;
5859
5860   ecleaz (y);                   /* start with a zero */
5861   p = y;                        /* point to our number */
5862   r = *d;                       /* get IBM exponent word */
5863   if (*d & (unsigned int) 0x8000)
5864     *p = 0xffff;                /* fill in our sign */
5865   ++p;                          /* bump pointer to our exponent word */
5866   r &= 0x7f00;                  /* strip the sign bit */
5867   r >>= 6;                      /* shift exponent word down 6 bits */
5868                                 /* in fact shift by 8 right and 2 left */
5869   r += EXONE - (0x41 << 2);     /* subtract IBM exponent offset */
5870                                 /* add our e type exponent offset */
5871   *p++ = r;                     /* to form our exponent */
5872
5873   *p++ = *d++ & 0xff;           /* now do the high order mantissa */
5874                                 /* strip off the IBM exponent and sign bits */
5875   if (mode != SFmode)           /* there are only 2 words in SFmode */
5876     {
5877       *p++ = *d++;              /* fill in the rest of our mantissa */
5878       *p++ = *d++;
5879     }
5880   *p = *d;
5881
5882   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5883     y[0] = y[E] = 0;
5884   else
5885     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
5886                               /* handle change in RADIX */
5887   emovo (y, e);
5888 }
5889
5890
5891
5892 /* Convert e type to IBM single/double precision.  */
5893
5894 static void
5895 etoibm (x, d, mode)
5896      unsigned EMUSHORT *x, *d;
5897      enum machine_mode mode;
5898 {
5899   unsigned EMUSHORT xi[NI];
5900   EMULONG exp;
5901   int rndsav;
5902
5903   emovi (x, xi);
5904   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));        /* adjust exponent for offsets */
5905                                                         /* round off to nearest or even */
5906   rndsav = rndprc;
5907   rndprc = 56;
5908   emdnorm (xi, 0, 0, exp, 64);
5909   rndprc = rndsav;
5910   toibm (xi, d, mode);
5911 }
5912
5913 static void
5914 toibm (x, y, mode)
5915      unsigned EMUSHORT *x, *y;
5916      enum machine_mode mode;
5917 {
5918   unsigned EMUSHORT i;
5919   unsigned EMUSHORT *p;
5920   int r;
5921
5922   p = x;
5923   *y = 0;
5924   if (*p++)
5925     *y = 0x8000;
5926   i = *p++;
5927   if (i == 0)
5928     {
5929       *y++ = 0;
5930       *y++ = 0;
5931       if (mode != SFmode)
5932         {
5933           *y++ = 0;
5934           *y++ = 0;
5935         }
5936       return;
5937     }
5938   r = i & 0x3;
5939   i >>= 2;
5940   if (i > 0x7f)
5941     {
5942       *y++ |= 0x7fff;
5943       *y++ = 0xffff;
5944       if (mode != SFmode)
5945         {
5946           *y++ = 0xffff;
5947           *y++ = 0xffff;
5948         }
5949 #ifdef ERANGE
5950       errno = ERANGE;
5951 #endif
5952       return;
5953     }
5954   i &= 0x7f;
5955   *y |= (i << 8);
5956   eshift (x, r + 5);
5957   *y++ |= x[M];
5958   *y++ = x[M + 1];
5959   if (mode != SFmode)
5960     {
5961       *y++ = x[M + 2];
5962       *y++ = x[M + 3];
5963     }
5964 }
5965 #endif /* IBM */
5966
5967
5968 #ifdef C4X
5969 /* Convert C4X single/double precision to e type.  */
5970
5971 static void
5972 c4xtoe (d, e, mode)
5973      unsigned EMUSHORT *d;
5974      unsigned EMUSHORT *e;
5975      enum machine_mode mode;
5976 {
5977   unsigned EMUSHORT y[NI];
5978   int r;
5979   int isnegative;
5980   int size;
5981   int i;
5982   int carry;
5983
5984   /* Short-circuit the zero case. */
5985   if ((d[0] == 0x8000)
5986       && (d[1] == 0x0000)
5987       && ((mode == QFmode) || ((d[2] == 0x0000) && (d[3] == 0x0000))))
5988     {
5989       e[0] = 0;
5990       e[1] = 0;
5991       e[2] = 0;
5992       e[3] = 0;
5993       e[4] = 0;
5994       e[5] = 0;
5995       return;
5996     }
5997
5998   ecleaz (y);                   /* start with a zero */
5999   r = d[0];                     /* get sign/exponent part */
6000   if (r & (unsigned int) 0x0080)
6001   {
6002      y[0] = 0xffff;             /* fill in our sign */
6003      isnegative = TRUE;
6004   }
6005   else
6006   {
6007      isnegative = FALSE;
6008   }
6009
6010   r >>= 8;                      /* Shift exponent word down 8 bits.  */
6011   if (r & 0x80)                 /* Make the exponent negative if it is. */
6012   {
6013      r = r | (~0 & ~0xff);
6014   }
6015
6016   if (isnegative)
6017   {
6018      /* Now do the high order mantissa.  We don't "or" on the high bit
6019         because it is 2 (not 1) and is handled a little differently
6020         below.  */
6021      y[M] = d[0] & 0x7f;
6022
6023      y[M+1] = d[1];
6024      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
6025      {
6026         y[M+2] = d[2];          /* Fill in the rest of our mantissa.  */
6027         y[M+3] = d[3];
6028         size = 4;
6029      }
6030      else
6031      {
6032         size = 2;
6033      }
6034      eshift(y, -8);
6035
6036      /* Now do the two's complement on the data.  */
6037
6038      carry = 1; /* Initially add 1 for the two's complement. */
6039      for (i=size + M; i > M; i--)
6040      {
6041         if (carry && (y[i] == 0x0000))
6042         {
6043            /* We overflowed into the next word, carry is the same.  */
6044            y[i] = carry ? 0x0000 : 0xffff;
6045         }
6046         else
6047         {
6048            /* No overflow, just invert and add carry.  */
6049            y[i] = ((~y[i]) + carry) & 0xffff;
6050            carry = 0;
6051         }
6052      }
6053
6054      if (carry)
6055      {
6056         eshift(y, -1);
6057         y[M+1] |= 0x8000;
6058         r++;
6059      }
6060      y[1] = r + EXONE;
6061   }
6062   else
6063   {
6064     /* Add our e type exponent offset to form our exponent.  */
6065      r += EXONE;
6066      y[1] = r;
6067
6068      /* Now do the high order mantissa strip off the exponent and sign
6069         bits and add the high 1 bit.  */
6070      y[M] = (d[0] & 0x7f) | 0x80;
6071
6072      y[M+1] = d[1];
6073      if (mode != QFmode)        /* There are only 2 words in QFmode.  */
6074      {
6075         y[M+2] = d[2];          /* Fill in the rest of our mantissa.  */
6076         y[M+3] = d[3];
6077      }
6078      eshift(y, -8);
6079   }
6080
6081   emovo (y, e);
6082 }
6083
6084
6085 /* Convert e type to C4X single/double precision.  */
6086
6087 static void
6088 etoc4x (x, d, mode)
6089      unsigned EMUSHORT *x, *d;
6090      enum machine_mode mode;
6091 {
6092   unsigned EMUSHORT xi[NI];
6093   EMULONG exp;
6094   int rndsav;
6095
6096   emovi (x, xi);
6097
6098   /* Adjust exponent for offsets. */
6099   exp = (EMULONG) xi[E] - (EXONE - 0x7f);
6100
6101   /* Round off to nearest or even. */
6102   rndsav = rndprc;
6103   rndprc = mode == QFmode ? 24 : 32;
6104   emdnorm (xi, 0, 0, exp, 64);
6105   rndprc = rndsav;
6106   toc4x (xi, d, mode);
6107 }
6108
6109 static void
6110 toc4x (x, y, mode)
6111      unsigned EMUSHORT *x, *y;
6112      enum machine_mode mode;
6113 {
6114   int i;
6115   int v;
6116   int carry;
6117
6118   /* Short-circuit the zero case */
6119   if ((x[0] == 0)       /* Zero exponent and sign */
6120       && (x[1] == 0)
6121       && (x[M] == 0)    /* The rest is for zero mantissa */
6122       && (x[M+1] == 0)
6123       /* Only check for double if necessary */
6124       && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
6125     {
6126       /* We have a zero.  Put it into the output and return. */
6127       *y++ = 0x8000;
6128       *y++ = 0x0000;
6129       if (mode != QFmode)
6130         {
6131           *y++ = 0x0000;
6132           *y++ = 0x0000;
6133         }
6134       return;
6135     }
6136
6137   *y = 0;
6138
6139   /* Negative number require a two's complement conversion of the
6140      mantissa. */
6141   if (x[0])
6142     {
6143       *y = 0x0080;
6144
6145       i = ((int) x[1]) - 0x7f;
6146
6147       /* Now add 1 to the inverted data to do the two's complement. */
6148       if (mode != QFmode)
6149         v = 4 + M;
6150       else
6151         v = 2 + M;
6152       carry = 1;
6153       while (v > M)
6154         {
6155           if (x[v] == 0x0000)
6156             {
6157               x[v] = carry ? 0x0000 : 0xffff;
6158             }
6159           else
6160             {
6161               x[v] = ((~x[v]) + carry) & 0xffff;
6162               carry = 0;
6163             }
6164           v--;
6165         }
6166
6167       /* The following is a special case.  The C4X negative float requires
6168          a zero in the high bit (because the format is (2 - x) x 2^m), so
6169          if a one is in that bit, we have to shift left one to get rid
6170          of it.  This only occurs if the number is -1 x 2^m. */
6171       if (x[M+1] & 0x8000)
6172         {
6173           /* This is the case of -1 x 2^m, we have to rid ourselves of the
6174              high sign bit and shift the exponent. */
6175           eshift(x, 1);
6176           i--;
6177         }
6178     }
6179   else
6180     {
6181       i = ((int) x[1]) - 0x7f;
6182     }
6183
6184   if ((i < -128) || (i > 127))
6185     {
6186       y[0] |= 0xff7f;
6187       y[1] = 0xffff;
6188       if (mode != QFmode)
6189         {
6190           y[2] = 0xffff;
6191           y[3] = 0xffff;
6192         }
6193 #ifdef ERANGE
6194       errno = ERANGE;
6195 #endif
6196       return;
6197     }
6198
6199   y[0] |= ((i & 0xff) << 8);
6200
6201   eshift (x, 8);
6202
6203   y[0] |= x[M] & 0x7f;
6204   y[1] = x[M + 1];
6205   if (mode != QFmode)
6206     {
6207       y[2] = x[M + 2];
6208       y[3] = x[M + 3];
6209     }
6210 }
6211 #endif /* C4X */
6212
6213 /* Output a binary NaN bit pattern in the target machine's format.  */
6214
6215 /* If special NaN bit patterns are required, define them in tm.h
6216    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
6217    patterns here.  */
6218 #ifdef TFMODE_NAN
6219 TFMODE_NAN;
6220 #else
6221 #ifdef IEEE
6222 unsigned EMUSHORT TFbignan[8] =
6223  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6224 unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6225 #endif
6226 #endif
6227
6228 #ifdef XFMODE_NAN
6229 XFMODE_NAN;
6230 #else
6231 #ifdef IEEE
6232 unsigned EMUSHORT XFbignan[6] =
6233  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6234 unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6235 #endif
6236 #endif
6237
6238 #ifdef DFMODE_NAN
6239 DFMODE_NAN;
6240 #else
6241 #ifdef IEEE
6242 unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6243 unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
6244 #endif
6245 #endif
6246
6247 #ifdef SFMODE_NAN
6248 SFMODE_NAN;
6249 #else
6250 #ifdef IEEE
6251 unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
6252 unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
6253 #endif
6254 #endif
6255
6256
6257 #ifdef NANS
6258 static void
6259 make_nan (nan, sign, mode)
6260      unsigned EMUSHORT *nan;
6261      int sign;
6262      enum machine_mode mode;
6263 {
6264   int n;
6265   unsigned EMUSHORT *p;
6266
6267   switch (mode)
6268     {
6269 /* Possibly the `reserved operand' patterns on a VAX can be
6270    used like NaN's, but probably not in the same way as IEEE.  */
6271 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6272     case TFmode:
6273       n = 8;
6274       if (REAL_WORDS_BIG_ENDIAN)
6275         p = TFbignan;
6276       else
6277         p = TFlittlenan;
6278       break;
6279
6280     case XFmode:
6281       n = 6;
6282       if (REAL_WORDS_BIG_ENDIAN)
6283         p = XFbignan;
6284       else
6285         p = XFlittlenan;
6286       break;
6287
6288     case DFmode:
6289       n = 4;
6290       if (REAL_WORDS_BIG_ENDIAN)
6291         p = DFbignan;
6292       else
6293         p = DFlittlenan;
6294       break;
6295
6296     case SFmode:
6297     case HFmode:
6298       n = 2;
6299       if (REAL_WORDS_BIG_ENDIAN)
6300         p = SFbignan;
6301       else
6302         p = SFlittlenan;
6303       break;
6304 #endif
6305
6306     default:
6307       abort ();
6308     }
6309   if (REAL_WORDS_BIG_ENDIAN)
6310     *nan++ = (sign << 15) | (*p++ & 0x7fff);
6311   while (--n != 0)
6312     *nan++ = *p++;
6313   if (! REAL_WORDS_BIG_ENDIAN)
6314     *nan = (sign << 15) | (*p & 0x7fff);
6315 }
6316 #endif /* NANS */
6317
6318 /* This is the inverse of the function `etarsingle' invoked by
6319    REAL_VALUE_TO_TARGET_SINGLE.  */
6320
6321 REAL_VALUE_TYPE
6322 ereal_unto_float (f)
6323      long f;
6324 {
6325   REAL_VALUE_TYPE r;
6326   unsigned EMUSHORT s[2];
6327   unsigned EMUSHORT e[NE];
6328
6329   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6330    This is the inverse operation to what the function `endian' does.  */
6331   if (REAL_WORDS_BIG_ENDIAN)
6332     {
6333       s[0] = (unsigned EMUSHORT) (f >> 16);
6334       s[1] = (unsigned EMUSHORT) f;
6335     }
6336   else
6337     {
6338       s[0] = (unsigned EMUSHORT) f;
6339       s[1] = (unsigned EMUSHORT) (f >> 16);
6340     }
6341   /* Convert and promote the target float to E-type. */
6342   e24toe (s, e);
6343   /* Output E-type to REAL_VALUE_TYPE. */
6344   PUT_REAL (e, &r);
6345   return r;
6346 }
6347
6348
6349 /* This is the inverse of the function `etardouble' invoked by
6350    REAL_VALUE_TO_TARGET_DOUBLE.  */
6351
6352 REAL_VALUE_TYPE
6353 ereal_unto_double (d)
6354      long d[];
6355 {
6356   REAL_VALUE_TYPE r;
6357   unsigned EMUSHORT s[4];
6358   unsigned EMUSHORT e[NE];
6359
6360   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6361   if (REAL_WORDS_BIG_ENDIAN)
6362     {
6363       s[0] = (unsigned EMUSHORT) (d[0] >> 16);
6364       s[1] = (unsigned EMUSHORT) d[0];
6365       s[2] = (unsigned EMUSHORT) (d[1] >> 16);
6366       s[3] = (unsigned EMUSHORT) d[1];
6367     }
6368   else
6369     {
6370       /* Target float words are little-endian.  */
6371       s[0] = (unsigned EMUSHORT) d[0];
6372       s[1] = (unsigned EMUSHORT) (d[0] >> 16);
6373       s[2] = (unsigned EMUSHORT) d[1];
6374       s[3] = (unsigned EMUSHORT) (d[1] >> 16);
6375     }
6376   /* Convert target double to E-type. */
6377   e53toe (s, e);
6378   /* Output E-type to REAL_VALUE_TYPE. */
6379   PUT_REAL (e, &r);
6380   return r;
6381 }
6382
6383
6384 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6385    This is somewhat like ereal_unto_float, but the input types
6386    for these are different.  */
6387
6388 REAL_VALUE_TYPE
6389 ereal_from_float (f)
6390      HOST_WIDE_INT f;
6391 {
6392   REAL_VALUE_TYPE r;
6393   unsigned EMUSHORT s[2];
6394   unsigned EMUSHORT e[NE];
6395
6396   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6397    This is the inverse operation to what the function `endian' does.  */
6398   if (REAL_WORDS_BIG_ENDIAN)
6399     {
6400       s[0] = (unsigned EMUSHORT) (f >> 16);
6401       s[1] = (unsigned EMUSHORT) f;
6402     }
6403   else
6404     {
6405       s[0] = (unsigned EMUSHORT) f;
6406       s[1] = (unsigned EMUSHORT) (f >> 16);
6407     }
6408   /* Convert and promote the target float to E-type.  */
6409   e24toe (s, e);
6410   /* Output E-type to REAL_VALUE_TYPE.  */
6411   PUT_REAL (e, &r);
6412   return r;
6413 }
6414
6415
6416 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6417    This is somewhat like ereal_unto_double, but the input types
6418    for these are different.
6419
6420    The DFmode is stored as an array of HOST_WIDE_INT in the target's
6421    data format, with no holes in the bit packing.  The first element
6422    of the input array holds the bits that would come first in the
6423    target computer's memory.  */
6424
6425 REAL_VALUE_TYPE
6426 ereal_from_double (d)
6427      HOST_WIDE_INT d[];
6428 {
6429   REAL_VALUE_TYPE r;
6430   unsigned EMUSHORT s[4];
6431   unsigned EMUSHORT e[NE];
6432
6433   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
6434   if (REAL_WORDS_BIG_ENDIAN)
6435     {
6436 #if HOST_BITS_PER_WIDE_INT == 32
6437       s[0] = (unsigned EMUSHORT) (d[0] >> 16);
6438       s[1] = (unsigned EMUSHORT) d[0];
6439       s[2] = (unsigned EMUSHORT) (d[1] >> 16);
6440       s[3] = (unsigned EMUSHORT) d[1];
6441 #else
6442       /* In this case the entire target double is contained in the
6443          first array element.  The second element of the input is
6444          ignored.  */
6445       s[0] = (unsigned EMUSHORT) (d[0] >> 48);
6446       s[1] = (unsigned EMUSHORT) (d[0] >> 32);
6447       s[2] = (unsigned EMUSHORT) (d[0] >> 16);
6448       s[3] = (unsigned EMUSHORT) d[0];
6449 #endif
6450     }
6451   else
6452     {
6453       /* Target float words are little-endian.  */
6454       s[0] = (unsigned EMUSHORT) d[0];
6455       s[1] = (unsigned EMUSHORT) (d[0] >> 16);
6456 #if HOST_BITS_PER_WIDE_INT == 32
6457       s[2] = (unsigned EMUSHORT) d[1];
6458       s[3] = (unsigned EMUSHORT) (d[1] >> 16);
6459 #else
6460       s[2] = (unsigned EMUSHORT) (d[0] >> 32);
6461       s[3] = (unsigned EMUSHORT) (d[0] >> 48);
6462 #endif
6463     }
6464   /* Convert target double to E-type.  */
6465   e53toe (s, e);
6466   /* Output E-type to REAL_VALUE_TYPE.  */
6467   PUT_REAL (e, &r);
6468   return r;
6469 }
6470
6471
6472 #if 0
6473 /* Convert target computer unsigned 64-bit integer to e-type.
6474    The endian-ness of DImode follows the convention for integers,
6475    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
6476
6477 static void
6478 uditoe (di, e)
6479      unsigned EMUSHORT *di;  /* Address of the 64-bit int.  */
6480      unsigned EMUSHORT *e;
6481 {
6482   unsigned EMUSHORT yi[NI];
6483   int k;
6484
6485   ecleaz (yi);
6486   if (WORDS_BIG_ENDIAN)
6487     {
6488       for (k = M; k < M + 4; k++)
6489         yi[k] = *di++;
6490     }
6491   else
6492     {
6493       for (k = M + 3; k >= M; k--)
6494         yi[k] = *di++;
6495     }
6496   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6497   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6498     ecleaz (yi);                /* it was zero */
6499   else
6500     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
6501   emovo (yi, e);
6502 }
6503
6504 /* Convert target computer signed 64-bit integer to e-type.  */
6505
6506 static void
6507 ditoe (di, e)
6508      unsigned EMUSHORT *di;  /* Address of the 64-bit int.  */
6509      unsigned EMUSHORT *e;
6510 {
6511   unsigned EMULONG acc;
6512   unsigned EMUSHORT yi[NI];
6513   unsigned EMUSHORT carry;
6514   int k, sign;
6515
6516   ecleaz (yi);
6517   if (WORDS_BIG_ENDIAN)
6518     {
6519       for (k = M; k < M + 4; k++)
6520         yi[k] = *di++;
6521     }
6522   else
6523     {
6524       for (k = M + 3; k >= M; k--)
6525         yi[k] = *di++;
6526     }
6527   /* Take absolute value */
6528   sign = 0;
6529   if (yi[M] & 0x8000)
6530     {
6531       sign = 1;
6532       carry = 0;
6533       for (k = M + 3; k >= M; k--)
6534         {
6535           acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6536           yi[k] = acc;
6537           carry = 0;
6538           if (acc & 0x10000)
6539             carry = 1;
6540         }
6541     }
6542   yi[E] = EXONE + 47;   /* exponent if normalize shift count were 0 */
6543   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6544     ecleaz (yi);                /* it was zero */
6545   else
6546     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
6547   emovo (yi, e);
6548   if (sign)
6549         eneg (e);
6550 }
6551
6552
6553 /* Convert e-type to unsigned 64-bit int.  */
6554
6555 static void
6556 etoudi (x, i)
6557      unsigned EMUSHORT *x;
6558      unsigned EMUSHORT *i;
6559 {
6560   unsigned EMUSHORT xi[NI];
6561   int j, k;
6562
6563   emovi (x, xi);
6564   if (xi[0])
6565     {
6566       xi[M] = 0;
6567       goto noshift;
6568     }
6569   k = (int) xi[E] - (EXONE - 1);
6570   if (k <= 0)
6571     {
6572       for (j = 0; j < 4; j++)
6573         *i++ = 0;
6574       return;
6575     }
6576   if (k > 64)
6577     {
6578       for (j = 0; j < 4; j++)
6579         *i++ = 0xffff;
6580       if (extra_warnings)
6581         warning ("overflow on truncation to integer");
6582       return;
6583     }
6584   if (k > 16)
6585     {
6586       /* Shift more than 16 bits: first shift up k-16 mod 16,
6587          then shift up by 16's.  */
6588       j = k - ((k >> 4) << 4);
6589       if (j == 0)
6590         j = 16;
6591       eshift (xi, j);
6592       if (WORDS_BIG_ENDIAN)
6593         *i++ = xi[M];
6594       else
6595         {
6596           i += 3;
6597           *i-- = xi[M];
6598         }
6599       k -= j;
6600       do
6601         {
6602           eshup6 (xi);
6603           if (WORDS_BIG_ENDIAN)
6604             *i++ = xi[M];
6605           else
6606             *i-- = xi[M];
6607         }
6608       while ((k -= 16) > 0);
6609     }
6610   else
6611     {
6612         /* shift not more than 16 bits */
6613       eshift (xi, k);
6614
6615 noshift:
6616
6617       if (WORDS_BIG_ENDIAN)
6618         {
6619           i += 3;
6620           *i-- = xi[M];
6621           *i-- = 0;
6622           *i-- = 0;
6623           *i = 0;
6624         }
6625       else
6626         {
6627           *i++ = xi[M];
6628           *i++ = 0;
6629           *i++ = 0;
6630           *i = 0;
6631         }
6632     }
6633 }
6634
6635
6636 /* Convert e-type to signed 64-bit int.  */
6637
6638 static void
6639 etodi (x, i)
6640      unsigned EMUSHORT *x;
6641      unsigned EMUSHORT *i;
6642 {
6643   unsigned EMULONG acc;
6644   unsigned EMUSHORT xi[NI];
6645   unsigned EMUSHORT carry;
6646   unsigned EMUSHORT *isave;
6647   int j, k;
6648
6649   emovi (x, xi);
6650   k = (int) xi[E] - (EXONE - 1);
6651   if (k <= 0)
6652     {
6653       for (j = 0; j < 4; j++)
6654         *i++ = 0;
6655       return;
6656     }
6657   if (k > 64)
6658     {
6659       for (j = 0; j < 4; j++)
6660         *i++ = 0xffff;
6661       if (extra_warnings)
6662         warning ("overflow on truncation to integer");
6663       return;
6664     }
6665   isave = i;
6666   if (k > 16)
6667     {
6668       /* Shift more than 16 bits: first shift up k-16 mod 16,
6669          then shift up by 16's.  */
6670       j = k - ((k >> 4) << 4);
6671       if (j == 0)
6672         j = 16;
6673       eshift (xi, j);
6674       if (WORDS_BIG_ENDIAN)
6675         *i++ = xi[M];
6676       else
6677         {
6678           i += 3;
6679           *i-- = xi[M];
6680         }
6681       k -= j;
6682       do
6683         {
6684           eshup6 (xi);
6685           if (WORDS_BIG_ENDIAN)
6686             *i++ = xi[M];
6687           else
6688             *i-- = xi[M];
6689         }
6690       while ((k -= 16) > 0);
6691     }
6692   else
6693     {
6694         /* shift not more than 16 bits */
6695       eshift (xi, k);
6696
6697       if (WORDS_BIG_ENDIAN)
6698         {
6699           i += 3;
6700           *i = xi[M];
6701           *i-- = 0;
6702           *i-- = 0;
6703           *i = 0;
6704         }
6705       else
6706         {
6707           *i++ = xi[M];
6708           *i++ = 0;
6709           *i++ = 0;
6710           *i = 0;
6711         }
6712     }
6713   /* Negate if negative */
6714   if (xi[0])
6715     {
6716       carry = 0;
6717       if (WORDS_BIG_ENDIAN)
6718         isave += 3;
6719       for (k = 0; k < 4; k++)
6720         {
6721           acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6722           if (WORDS_BIG_ENDIAN)
6723             *isave-- = acc;
6724           else
6725             *isave++ = acc;
6726           carry = 0;
6727           if (acc & 0x10000)
6728             carry = 1;
6729         }
6730     }
6731 }
6732
6733
6734 /* Longhand square root routine.  */
6735
6736
6737 static int esqinited = 0;
6738 static unsigned short sqrndbit[NI];
6739
6740 static void
6741 esqrt (x, y)
6742      unsigned EMUSHORT *x, *y;
6743 {
6744   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6745   EMULONG m, exp;
6746   int i, j, k, n, nlups;
6747
6748   if (esqinited == 0)
6749     {
6750       ecleaz (sqrndbit);
6751       sqrndbit[NI - 2] = 1;
6752       esqinited = 1;
6753     }
6754   /* Check for arg <= 0 */
6755   i = ecmp (x, ezero);
6756   if (i <= 0)
6757     {
6758       if (i == -1)
6759         {
6760           mtherr ("esqrt", DOMAIN);
6761           eclear (y);
6762         }
6763       else
6764         emov (x, y);
6765       return;
6766     }
6767
6768 #ifdef INFINITY
6769   if (eisinf (x))
6770     {
6771       eclear (y);
6772       einfin (y);
6773       return;
6774     }
6775 #endif
6776   /* Bring in the arg and renormalize if it is denormal.  */
6777   emovi (x, xx);
6778   m = (EMULONG) xx[1];          /* local long word exponent */
6779   if (m == 0)
6780     m -= enormlz (xx);
6781
6782   /* Divide exponent by 2 */
6783   m -= 0x3ffe;
6784   exp = (unsigned short) ((m / 2) + 0x3ffe);
6785
6786   /* Adjust if exponent odd */
6787   if ((m & 1) != 0)
6788     {
6789       if (m > 0)
6790         exp += 1;
6791       eshdn1 (xx);
6792     }
6793
6794   ecleaz (sq);
6795   ecleaz (num);
6796   n = 8;                        /* get 8 bits of result per inner loop */
6797   nlups = rndprc;
6798   j = 0;
6799
6800   while (nlups > 0)
6801     {
6802       /* bring in next word of arg */
6803       if (j < NE)
6804         num[NI - 1] = xx[j + 3];
6805       /* Do additional bit on last outer loop, for roundoff.  */
6806       if (nlups <= 8)
6807         n = nlups + 1;
6808       for (i = 0; i < n; i++)
6809         {
6810           /* Next 2 bits of arg */
6811           eshup1 (num);
6812           eshup1 (num);
6813           /* Shift up answer */
6814           eshup1 (sq);
6815           /* Make trial divisor */
6816           for (k = 0; k < NI; k++)
6817             temp[k] = sq[k];
6818           eshup1 (temp);
6819           eaddm (sqrndbit, temp);
6820           /* Subtract and insert answer bit if it goes in */
6821           if (ecmpm (temp, num) <= 0)
6822             {
6823               esubm (temp, num);
6824               sq[NI - 2] |= 1;
6825             }
6826         }
6827       nlups -= n;
6828       j += 1;
6829     }
6830
6831   /* Adjust for extra, roundoff loop done.  */
6832   exp += (NBITS - 1) - rndprc;
6833
6834   /* Sticky bit = 1 if the remainder is nonzero.  */
6835   k = 0;
6836   for (i = 3; i < NI; i++)
6837     k |= (int) num[i];
6838
6839   /* Renormalize and round off.  */
6840   emdnorm (sq, k, 0, exp, 64);
6841   emovo (sq, y);
6842 }
6843 #endif
6844 #endif /* EMU_NON_COMPILE not defined */
6845 \f
6846 /* Return the binary precision of the significand for a given
6847    floating point mode.  The mode can hold an integer value
6848    that many bits wide, without losing any bits.  */
6849
6850 unsigned int
6851 significand_size (mode)
6852      enum machine_mode mode;
6853 {
6854
6855 /* Don't test the modes, but their sizes, lest this
6856    code won't work for BITS_PER_UNIT != 8 .  */
6857
6858 switch (GET_MODE_BITSIZE (mode))
6859   {
6860   case 32:
6861
6862 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6863     return 56;
6864 #endif
6865
6866     return 24;
6867
6868   case 64:
6869 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6870     return 53;
6871 #else
6872 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6873     return 56;
6874 #else
6875 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6876     return 56;
6877 #else
6878 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6879     return 56;
6880 #else
6881     abort ();
6882 #endif
6883 #endif
6884 #endif
6885 #endif
6886
6887   case 96:
6888     return 64;
6889   case 128:
6890     return 113;
6891
6892   default:
6893     abort ();
6894   }
6895 }