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