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