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