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