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