gcc/real.c

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