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