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