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