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