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