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