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