gcc/real.c

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