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