gcc/real.c

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