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