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