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