gcc/real.c

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