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