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