gcc/real.c

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