gcc/real.c

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