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