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