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