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