gcc/real.c

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