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