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 "config.h"
  25
  26 #include "tree.h"
  27
  28 /* To enable support of XFmode extended real floating point, define
  29 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  30
  31 To support cross compilation between IEEE and VAX floating
  32 point formats, define REAL_ARITHMETIC in the tm.h file.
  33
  34 In either case the machine files (tm.h) must not contain any code
  35 that tries to use host floating point arithmetic to convert
  36 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  37 etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  38 be intelligible to the host computer's native arithmetic.
  39
  40 The emulator defaults to the host's floating point format so that
  41 its decimal conversion functions can be used if desired (see
  42 real.h).
  43
  44 The first part of this file interfaces gcc to ieee.c, which is a
  45 floating point arithmetic suite that was not written with gcc in
  46 mind.  The interface is followed by ieee.c itself and related
  47 items. Avoid changing ieee.c unless you have suitable test
  48 programs available.  A special version of the PARANOIA floating
  49 point arithmetic tester, modified for this purpose, can be found
  50 on usc.edu : /pub/C-numanal/ieeetest.zoo.  Some tutorial
  51 information on ieee.c is given in my book: S. L. Moshier,
  52 _Methods and Programs for Mathematical Functions_, Prentice-Hall
  53 or Simon & Schuster Int'l, 1989.  A library of XFmode elementary
  54 transcendental functions can be obtained by ftp from
  55 research.att.com: netlib/cephes/ldouble.shar.Z  */
  56
  57 /* Type of computer arithmetic.
  58  * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
  59  * The following modification converts gcc macros into the ones
  60  * used by ieee.c.
  61  *
  62  * Note: long double formats differ between IBMPC and MIEEE
  63  * by more than just endian-ness.
  64  */
  65
  66 /* REAL_ARITHMETIC defined means that macros in real.h are
  67    defined to call emulator functions.  */
  68 #ifdef REAL_ARITHMETIC
  69
  70 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  71 /* PDP-11, Pro350, VAX: */
  72 #define DEC 1
  73 #else /* it's not VAX */
  74 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  75 #if WORDS_BIG_ENDIAN
  76 /* Motorola IEEE, high order words come first (Sun workstation): */
  77 #define MIEEE 1
  78 #else /* not big-endian */
  79 /* Intel IEEE, low order words come first:
  80  */
  81 #define IBMPC 1
  82 #endif /*  big-endian */
  83 #else /* it's not IEEE either */
  84 /* UNKnown arithmetic.  We don't support this and can't go on. */
  85 unknown arithmetic type
  86 #define UNK 1
  87 #endif /* not IEEE */
  88 #endif /* not VAX */
  89
  90 #else
  91 /* REAL_ARITHMETIC not defined means that the *host's* data
  92    structure will be used.  It may differ by endian-ness from the
  93    target machine's structure and will get its ends swapped
  94    accordingly (but not here).  Probably only the decimal <-> binary
  95    functions in this file will actually be used in this case.  */
  96 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  97 #define DEC 1
  98 #else /* it's not VAX */
  99 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 100 #ifdef HOST_WORDS_BIG_ENDIAN
 101 #define MIEEE 1
 102 #else /* not big-endian */
 103 #define IBMPC 1
 104 #endif /*  big-endian */
 105 #else /* it's not IEEE either */
 106 unknown arithmetic type
 107 #define UNK 1
 108 #endif /* not IEEE */
 109 #endif /* not VAX */
 110
 111 #endif /* REAL_ARITHMETIC not defined */
 112
 113 /* Define for support of infinity.  */
 114 #ifndef DEC
 115 #define INFINITY
 116 #endif
 117
 118
 119 /* ehead.h
 120  *
 121  * Include file for extended precision arithmetic programs.
 122  */
 123
 124 /* Number of 16 bit words in external e type format */
 125 #define NE 6
 126
 127 /* Number of 16 bit words in internal format */
 128 #define NI (NE+3)
 129
 130 /* Array offset to exponent */
 131 #define E 1
 132
 133 /* Array offset to high guard word */
 134 #define M 2
 135
 136 /* Number of bits of precision */
 137 #define NBITS ((NI-4)*16)
 138
 139 /* Maximum number of decimal digits in ASCII conversion
 140  * = NBITS*log10(2)
 141  */
 142 #define NDEC (NBITS*8/27)
 143
 144 /* The exponent of 1.0 */
 145 #define EXONE (0x3fff)
 146
 147 /* Find a host integer type that is at least 16 bits wide,
 148    and another type at least twice whatever that size is. */
 149
 150 #if HOST_BITS_PER_CHAR >= 16
 151 #define EMUSHORT char
 152 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
 153 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
 154 #else
 155 #if HOST_BITS_PER_SHORT >= 16
 156 #define EMUSHORT short
 157 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
 158 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
 159 #else
 160 #if HOST_BITS_PER_INT >= 16
 161 #define EMUSHORT int
 162 #define EMUSHORT_SIZE HOST_BITS_PER_INT
 163 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
 164 #else
 165 #if HOST_BITS_PER_LONG >= 16
 166 #define EMUSHORT long
 167 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 168 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 169 #else
 170 /*  You will have to modify this program to have a smaller unit size. */
 171 #define EMU_NON_COMPILE
 172 #endif
 173 #endif
 174 #endif
 175 #endif
 176
 177 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
 178 #define EMULONG short
 179 #else
 180 #if HOST_BITS_PER_INT >= EMULONG_SIZE
 181 #define EMULONG int
 182 #else
 183 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
 184 #define EMULONG long
 185 #else
 186 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
 187 #define EMULONG long long int
 188 #else
 189 /*  You will have to modify this program to have a smaller unit size. */
 190 #define EMU_NON_COMPILE
 191 #endif
 192 #endif
 193 #endif
 194 #endif
 195
 196
 197 /* The host interface doesn't work if no 16-bit size exists. */
 198 #if EMUSHORT_SIZE != 16
 199 #define EMU_NON_COMPILE
 200 #endif
 201
 202 /* OK to continue compilation. */
 203 #ifndef EMU_NON_COMPILE
 204
 205 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
 206    In GET_REAL and PUT_REAL, r and e are pointers.
 207    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
 208    in memory, with no holes.  */
 209
 210 #if LONG_DOUBLE_TYPE_SIZE == 96
 211 #define GET_REAL(r,e) bcopy (r, e, 2*NE)
 212 #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
 213 #else /* no XFmode */
 214
 215 #ifdef REAL_ARITHMETIC
 216 /* Emulator uses target format internally
 217    but host stores it in host endian-ness. */
 218
 219 #if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
 220 #define GET_REAL(r,e) e53toe ((r), (e))
 221 #define PUT_REAL(e,r) etoe53 ((e), (r))
 222
 223 #else /* endian-ness differs */
 224 /* emulator uses target endian-ness internally */
 225 #define GET_REAL(r,e)           \
 226 do { EMUSHORT w[4];             \
 227  w[3] = ((EMUSHORT *) r)[0];    \
 228  w[2] = ((EMUSHORT *) r)[1];    \
 229  w[1] = ((EMUSHORT *) r)[2];    \
 230  w[0] = ((EMUSHORT *) r)[3];    \
 231  e53toe (w, (e)); } while (0)
 232
 233 #define PUT_REAL(e,r)           \
 234 do { EMUSHORT w[4];             \
 235  etoe53 ((e), w);               \
 236  *((EMUSHORT *) r) = w[3];      \
 237  *((EMUSHORT *) r + 1) = w[2];  \
 238  *((EMUSHORT *) r + 2) = w[1];  \
 239  *((EMUSHORT *) r + 3) = w[0]; } while (0)
 240
 241 #endif /* endian-ness differs */
 242
 243 #else /* not REAL_ARITHMETIC */
 244
 245 /* emulator uses host format */
 246 #define GET_REAL(r,e) e53toe ((r), (e))
 247 #define PUT_REAL(e,r) etoe53 ((e), (r))
 248
 249 #endif /* not REAL_ARITHMETIC */
 250 #endif /* no XFmode */
 251
 252 void pedwarn ();
 253 int ecmp (), enormlz (), eshift (), eisneg (), eisinf ();
 254 void eadd (), esub (), emul (), ediv ();
 255 void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
 256 void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
 257 void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
 258 void eround (), ereal_to_decimal ();
 259 void esqrt (), elog (), eexp (), etanh (), epow ();
 260 void asctoe (), asctoe24 (), asctoe53 (), asctoe64 ();
 261 void etoasc (), e24toasc (), e53toasc (), e64toasc ();
 262 void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
 263 void mtherr ();
 264 unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 265 unsigned EMUSHORT elog2[], esqrt2[];
 266
 267 /* Pack output array with 32-bit numbers obtained from
 268    array containing 16-bit numbers, swapping ends if required. */
 269 void
 270 endian (e, x, mode)
 271      unsigned EMUSHORT e[];
 272      long x[];
 273      enum machine_mode mode;
 274 {
 275   unsigned long th, t;
 276
 277 #if WORDS_BIG_ENDIAN
 278   switch (mode)
 279     {
 280
 281     case XFmode:
 282
 283       /* Swap halfwords in the third long. */
 284       th = (unsigned long) e[4] & 0xffff;
 285       t = (unsigned long) e[5] & 0xffff;
 286       t |= th << 16;
 287       x[2] = (long) t;
 288       /* fall into the double case */
 289
 290     case DFmode:
 291
 292       /* swap halfwords in the second word */
 293       th = (unsigned long) e[2] & 0xffff;
 294       t = (unsigned long) e[3] & 0xffff;
 295       t |= th << 16;
 296       x[1] = (long) t;
 297       /* fall into the float case */
 298
 299     case SFmode:
 300
 301       /* swap halfwords in the first word */
 302       th = (unsigned long) e[0] & 0xffff;
 303       t = (unsigned long) e[1] & 0xffff;
 304       t |= th << 16;
 305       x[0] = t;
 306       break;
 307
 308     default:
 309       abort ();
 310     }
 311
 312 #else
 313
 314   /* Pack the output array without swapping. */
 315
 316   switch (mode)
 317     {
 318
 319     case XFmode:
 320
 321       /* Pack the third long.
 322          Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
 323          in it.  */
 324       th = (unsigned long) e[5] & 0xffff;
 325       t = (unsigned long) e[4] & 0xffff;
 326       t |= th << 16;
 327       x[2] = (long) t;
 328       /* fall into the double case */
 329
 330     case DFmode:
 331
 332       /* pack the second long */
 333       th = (unsigned long) e[3] & 0xffff;
 334       t = (unsigned long) e[2] & 0xffff;
 335       t |= th << 16;
 336       x[1] = (long) t;
 337       /* fall into the float case */
 338
 339     case SFmode:
 340
 341       /* pack the first long */
 342       th = (unsigned long) e[1] & 0xffff;
 343       t = (unsigned long) e[0] & 0xffff;
 344       t |= th << 16;
 345       x[0] = t;
 346       break;
 347
 348     default:
 349       abort ();
 350     }
 351
 352 #endif
 353 }
 354
 355
 356 /* This is the implementation of the REAL_ARITHMETIC macro.
 357  */
 358 void
 359 earith (value, icode, r1, r2)
 360      REAL_VALUE_TYPE *value;
 361      int icode;
 362      REAL_VALUE_TYPE *r1;
 363      REAL_VALUE_TYPE *r2;
 364 {
 365   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
 366   enum tree_code code;
 367
 368   GET_REAL (r1, d1);
 369   GET_REAL (r2, d2);
 370   code = (enum tree_code) icode;
 371   switch (code)
 372     {
 373     case PLUS_EXPR:
 374       eadd (d2, d1, v);
 375       break;
 376
 377     case MINUS_EXPR:
 378       esub (d2, d1, v);         /* d1 - d2 */
 379       break;
 380
 381     case MULT_EXPR:
 382       emul (d2, d1, v);
 383       break;
 384
 385     case RDIV_EXPR:
 386 #ifndef REAL_INFINITY
 387       if (ecmp (d2, ezero) == 0)
 388         abort ();
 389 #endif
 390       ediv (d2, d1, v); /* d1/d2 */
 391       break;
 392
 393     case MIN_EXPR:              /* min (d1,d2) */
 394       if (ecmp (d1, d2) < 0)
 395         emov (d1, v);
 396       else
 397         emov (d2, v);
 398       break;
 399
 400     case MAX_EXPR:              /* max (d1,d2) */
 401       if (ecmp (d1, d2) > 0)
 402         emov (d1, v);
 403       else
 404         emov (d2, v);
 405       break;
 406     default:
 407       emov (ezero, v);
 408       break;
 409     }
 410 PUT_REAL (v, value);
 411 }
 412
 413
 414 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
 415  * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x))
 416  */
 417 REAL_VALUE_TYPE
 418 etrunci (x)
 419      REAL_VALUE_TYPE x;
 420 {
 421   unsigned EMUSHORT f[NE], g[NE];
 422   REAL_VALUE_TYPE r;
 423   long l;
 424
 425   GET_REAL (&x, g);
 426   eifrac (g, &l, f);
 427   ltoe (&l, g);
 428   PUT_REAL (g, &r);
 429   return (r);
 430 }
 431
 432
 433 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
 434  * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x))
 435  */
 436 REAL_VALUE_TYPE
 437 etruncui (x)
 438      REAL_VALUE_TYPE x;
 439 {
 440   unsigned EMUSHORT f[NE], g[NE];
 441   REAL_VALUE_TYPE r;
 442   unsigned long l;
 443
 444   GET_REAL (&x, g);
 445   euifrac (g, &l, f);
 446   ultoe (&l, g);
 447   PUT_REAL (g, &r);
 448   return (r);
 449 }
 450
 451
 452 /* This is the REAL_VALUE_ATOF function.
 453  * It converts a decimal string to binary, rounding off
 454  * as indicated by the machine_mode argument.  Then it
 455  * promotes the rounded value to REAL_VALUE_TYPE.
 456  */
 457 REAL_VALUE_TYPE
 458 ereal_atof (s, t)
 459      char *s;
 460      enum machine_mode t;
 461 {
 462   unsigned EMUSHORT tem[NE], e[NE];
 463   REAL_VALUE_TYPE r;
 464
 465   switch (t)
 466     {
 467     case SFmode:
 468       asctoe24 (s, tem);
 469       e24toe (tem, e);
 470       break;
 471     case DFmode:
 472       asctoe53 (s, tem);
 473       e53toe (tem, e);
 474       break;
 475     case XFmode:
 476       asctoe64 (s, tem);
 477       e64toe (tem, e);
 478       break;
 479     default:
 480       asctoe (s, e);
 481     }
 482   PUT_REAL (e, &r);
 483   return (r);
 484 }
 485
 486
 487 /* Expansion of REAL_NEGATE.
 488  */
 489 REAL_VALUE_TYPE
 490 ereal_negate (x)
 491      REAL_VALUE_TYPE x;
 492 {
 493   unsigned EMUSHORT e[NE];
 494   REAL_VALUE_TYPE r;
 495
 496   GET_REAL (&x, e);
 497   eneg (e);
 498   PUT_REAL (e, &r);
 499   return (r);
 500 }
 501
 502
 503 /* Round real to int
 504  * implements REAL_VALUE_FIX (x) (eroundi (x))
 505  * The type of rounding is left unspecified by real.h.
 506  * It is implemented here as round to nearest (add .5 and chop).
 507  */
 508 int
 509 eroundi (x)
 510      REAL_VALUE_TYPE x;
 511 {
 512   unsigned EMUSHORT f[NE], g[NE];
 513   EMULONG l;
 514
 515   GET_REAL (&x, f);
 516   eround (f, g);
 517   eifrac (g, &l, f);
 518   return ((int) l);
 519 }
 520
 521 /* Round real to nearest unsigned int
 522  * implements  REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
 523  * Negative input returns zero.
 524  * The type of rounding is left unspecified by real.h.
 525  * It is implemented here as round to nearest (add .5 and chop).
 526  */
 527 unsigned int
 528 eroundui (x)
 529      REAL_VALUE_TYPE x;
 530 {
 531   unsigned EMUSHORT f[NE], g[NE];
 532   unsigned EMULONG l;
 533
 534   GET_REAL (&x, f);
 535   eround (f, g);
 536   euifrac (g, &l, f);
 537   return ((unsigned int)l);
 538 }
 539
 540
 541 /* REAL_VALUE_FROM_INT macro.
 542  */
 543 void
 544 ereal_from_int (d, i, j)
 545      REAL_VALUE_TYPE *d;
 546      long i, j;
 547 {
 548   unsigned EMUSHORT df[NE], dg[NE];
 549   long low, high;
 550   int sign;
 551
 552   sign = 0;
 553   low = i;
 554   if ((high = j) < 0)
 555     {
 556       sign = 1;
 557       /* complement and add 1 */
 558       high = ~high;
 559       if (low)
 560         low = -low;
 561       else
 562         high += 1;
 563     }
 564   eldexp (eone, HOST_BITS_PER_LONG, df);
 565   ultoe (&high, dg);
 566   emul (dg, df, dg);
 567   ultoe (&low, df);
 568   eadd (df, dg, dg);
 569   if (sign)
 570     eneg (dg);
 571   PUT_REAL (dg, d);
 572 }
 573
 574
 575 /* REAL_VALUE_FROM_UNSIGNED_INT macro.
 576  */
 577 void
 578 ereal_from_uint (d, i, j)
 579      REAL_VALUE_TYPE *d;
 580      unsigned long i, j;
 581 {
 582   unsigned EMUSHORT df[NE], dg[NE];
 583   unsigned long low, high;
 584
 585   low = i;
 586   high = j;
 587   eldexp (eone, HOST_BITS_PER_LONG, df);
 588   ultoe (&high, dg);
 589   emul (dg, df, dg);
 590   ultoe (&low, df);
 591   eadd (df, dg, dg);
 592   PUT_REAL (dg, d);
 593 }
 594
 595
 596 /* REAL_VALUE_TO_INT macro
 597  */
 598 void
 599 ereal_to_int (low, high, rr)
 600      long *low, *high;
 601      REAL_VALUE_TYPE rr;
 602 {
 603   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
 604   int s;
 605
 606   GET_REAL (&rr, d);
 607   /* convert positive value */
 608   s = 0;
 609   if (eisneg (d))
 610     {
 611       eneg (d);
 612       s = 1;
 613     }
 614   eldexp (eone, HOST_BITS_PER_LONG, df);
 615   ediv (df, d, dg);             /* dg = d / 2^32 is the high word */
 616   euifrac (dg, high, dh);
 617   emul (df, dh, dg);            /* fractional part is the low word */
 618   euifrac (dg, low, dh);
 619   if (s)
 620     {
 621       /* complement and add 1 */
 622       *high = ~(*high);
 623       if (*low)
 624         *low = -(*low);
 625       else
 626         *high += 1;
 627     }
 628 }
 629
 630
 631 /* REAL_VALUE_LDEXP macro.
 632  */
 633 REAL_VALUE_TYPE
 634 ereal_ldexp (x, n)
 635      REAL_VALUE_TYPE x;
 636      int n;
 637 {
 638   unsigned EMUSHORT e[NE], y[NE];
 639   REAL_VALUE_TYPE r;
 640
 641   GET_REAL (&x, e);
 642   eldexp (e, n, y);
 643   PUT_REAL (y, &r);
 644   return (r);
 645 }
 646
 647 /* These routines are conditionally compiled because functions
 648  * of the same names may be defined in fold-const.c.  */
 649 #ifdef REAL_ARITHMETIC
 650
 651 /* Check for infinity in a REAL_VALUE_TYPE. */
 652 int
 653 target_isinf (x)
 654      REAL_VALUE_TYPE x;
 655 {
 656   unsigned EMUSHORT e[NE];
 657
 658 #ifdef INFINITY
 659   GET_REAL (&x, e);
 660   return (eisinf (e));
 661 #else
 662   return 0;
 663 #endif
 664 }
 665
 666
 667 /* Check whether an IEEE double precision number is a NaN. */
 668
 669 int
 670 target_isnan (x)
 671      REAL_VALUE_TYPE x;
 672 {
 673   return (0);
 674 }
 675
 676
 677 /* Check for a negative IEEE double precision number.
 678  * this means strictly less than zero, not -0.
 679  */
 680
 681 int
 682 target_negative (x)
 683      REAL_VALUE_TYPE x;
 684 {
 685   unsigned EMUSHORT e[NE];
 686
 687   GET_REAL (&x, e);
 688   if (ecmp (e, ezero) < 0)
 689     return (1);
 690   return (0);
 691 }
 692
 693 /* Expansion of REAL_VALUE_TRUNCATE.
 694  * The result is in floating point, rounded to nearest or even.
 695  */
 696 REAL_VALUE_TYPE
 697 real_value_truncate (mode, arg)
 698      enum machine_mode mode;
 699      REAL_VALUE_TYPE arg;
 700 {
 701   unsigned EMUSHORT e[NE], t[NE];
 702   REAL_VALUE_TYPE r;
 703
 704   GET_REAL (&arg, e);
 705   eclear (t);
 706   switch (mode)
 707     {
 708     case XFmode:
 709       etoe64 (e, t);
 710       e64toe (t, t);
 711       break;
 712
 713     case DFmode:
 714       etoe53 (e, t);
 715       e53toe (t, t);
 716       break;
 717
 718     case SFmode:
 719       etoe24 (e, t);
 720       e24toe (t, t);
 721       break;
 722
 723     case SImode:
 724       r = etrunci (e);
 725       return (r);
 726
 727     default:
 728       abort ();
 729     }
 730   PUT_REAL (t, &r);
 731   return (r);
 732 }
 733
 734 #endif /* REAL_ARITHMETIC defined */
 735
 736 /* Target values are arrays of host longs. A long is guaranteed
 737    to be at least 32 bits wide. */
 738 void
 739 etarldouble (r, l)
 740      REAL_VALUE_TYPE r;
 741      long l[];
 742 {
 743   unsigned EMUSHORT e[NE];
 744
 745   GET_REAL (&r, e);
 746   etoe64 (e, e);
 747   endian (e, l, XFmode);
 748 }
 749
 750 void
 751 etardouble (r, l)
 752      REAL_VALUE_TYPE r;
 753      long l[];
 754 {
 755   unsigned EMUSHORT e[NE];
 756
 757   GET_REAL (&r, e);
 758   etoe53 (e, e);
 759   endian (e, l, DFmode);
 760 }
 761
 762 long
 763 etarsingle (r)
 764      REAL_VALUE_TYPE r;
 765 {
 766   unsigned EMUSHORT e[NE];
 767   unsigned long l;
 768
 769   GET_REAL (&r, e);
 770   etoe24 (e, e);
 771   endian (e, &l, SFmode);
 772   return ((long) l);
 773 }
 774
 775 void
 776 ereal_to_decimal (x, s)
 777      REAL_VALUE_TYPE x;
 778      char *s;
 779 {
 780   unsigned EMUSHORT e[NE];
 781
 782   GET_REAL (&x, e);
 783   etoasc (e, s, 20);
 784 }
 785
 786 int
 787 ereal_cmp (x, y)
 788      REAL_VALUE_TYPE x, y;
 789 {
 790   unsigned EMUSHORT ex[NE], ey[NE];
 791
 792   GET_REAL (&x, ex);
 793   GET_REAL (&y, ey);
 794   return (ecmp (ex, ey));
 795 }
 796
 797 int
 798 ereal_isneg (x)
 799      REAL_VALUE_TYPE x;
 800 {
 801   unsigned EMUSHORT ex[NE];
 802
 803   GET_REAL (&x, ex);
 804   return (eisneg (ex));
 805 }
 806
 807 /* End of REAL_ARITHMETIC interface */
 808
 809 /*                                                      ieee.c
 810  *
 811  *    Extended precision IEEE binary floating point arithmetic routines
 812  *
 813  * Numbers are stored in C language as arrays of 16-bit unsigned
 814  * short integers.  The arguments of the routines are pointers to
 815  * the arrays.
 816  *
 817  *
 818  * External e type data structure, simulates Intel 8087 chip
 819  * temporary real format but possibly with a larger significand:
 820  *
 821  *      NE-1 significand words  (least significant word first,
 822  *                               most significant bit is normally set)
 823  *      exponent                (value = EXONE for 1.0,
 824  *                              top bit is the sign)
 825  *
 826  *
 827  * Internal data structure of a number (a "word" is 16 bits):
 828  *
 829  * ei[0]        sign word       (0 for positive, 0xffff for negative)
 830  * ei[1]        biased exponent (value = EXONE for the number 1.0)
 831  * ei[2]        high guard word (always zero after normalization)
 832  * ei[3]
 833  * to ei[NI-2]  significand     (NI-4 significand words,
 834  *                               most significant word first,
 835  *                               most significant bit is set)
 836  * ei[NI-1]     low guard word  (0x8000 bit is rounding place)
 837  *
 838  *
 839  *
 840  *              Routines for external format numbers
 841  *
 842  *      asctoe (string, e)      ASCII string to extended double e type
 843  *      asctoe64 (string, &d)   ASCII string to long double
 844  *      asctoe53 (string, &d)   ASCII string to double
 845  *      asctoe24 (string, &f)   ASCII string to single
 846  *      asctoeg (string, e, prec) ASCII string to specified precision
 847  *      e24toe (&f, e)          IEEE single precision to e type
 848  *      e53toe (&d, e)          IEEE double precision to e type
 849  *      e64toe (&d, e)          IEEE long double precision to e type
 850  *      eabs (e)                        absolute value
 851  *      eadd (a, b, c)          c = b + a
 852  *      eclear (e)              e = 0
 853  *      ecmp (a, b)             compare a to b, return 1, 0, or -1
 854  *      ediv (a, b, c)          c = b / a
 855  *      efloor (a, b)           truncate to integer, toward -infinity
 856  *      efrexp (a, exp, s)      extract exponent and significand
 857  *      eifrac (e, &l, frac)    e to long integer and e type fraction
 858  *      euifrac (e, &l, frac)   e to unsigned long integer and e type fraction
 859  *      einfin (e)              set e to infinity, leaving its sign alone
 860  *      eldexp (a, n, b)        multiply by 2**n
 861  *      emov (a, b)             b = a
 862  *      emul (a, b, c)          c = b * a
 863  *      eneg (e)                        e = -e
 864  *      eround (a, b)           b = nearest integer value to a
 865  *      esub (a, b, c)          c = b - a
 866  *      e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
 867  *      e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
 868  *      e64toasc (&d, str, n)   long double to ASCII string
 869  *      etoasc (e, str, n)      e to ASCII string, n digits after decimal
 870  *      etoe24 (e, &f)          convert e type to IEEE single precision
 871  *      etoe53 (e, &d)          convert e type to IEEE double precision
 872  *      etoe64 (e, &d)          convert e type to IEEE long double precision
 873  *      ltoe (&l, e)            long (32 bit) integer to e type
 874  *      ultoe (&l, e)           unsigned long (32 bit) integer to e type
 875  *      eisneg (e)              1 if sign bit of e != 0, else 0
 876  *      eisinf (e)              1 if e has maximum exponent
 877  *
 878  *
 879  *              Routines for internal format numbers
 880  *
 881  *      eaddm (ai, bi)          add significands, bi = bi + ai
 882  *      ecleaz (ei)             ei = 0
 883  *      ecleazs (ei)            set ei = 0 but leave its sign alone
 884  *      ecmpm (ai, bi)          compare significands, return 1, 0, or -1
 885  *      edivm (ai, bi)          divide  significands, bi = bi / ai
 886  *      emdnorm (ai,l,s,exp)    normalize and round off
 887  *      emovi (a, ai)           convert external a to internal ai
 888  *      emovo (ai, a)           convert internal ai to external a
 889  *      emovz (ai, bi)          bi = ai, low guard word of bi = 0
 890  *      emulm (ai, bi)          multiply significands, bi = bi * ai
 891  *      enormlz (ei)            left-justify the significand
 892  *      eshdn1 (ai)             shift significand and guards down 1 bit
 893  *      eshdn8 (ai)             shift down 8 bits
 894  *      eshdn6 (ai)             shift down 16 bits
 895  *      eshift (ai, n)          shift ai n bits up (or down if n < 0)
 896  *      eshup1 (ai)             shift significand and guards up 1 bit
 897  *      eshup8 (ai)             shift up 8 bits
 898  *      eshup6 (ai)             shift up 16 bits
 899  *      esubm (ai, bi)          subtract significands, bi = bi - ai
 900  *
 901  *
 902  * The result is always normalized and rounded to NI-4 word precision
 903  * after each arithmetic operation.
 904  *
 905  * Exception flags and NaNs are NOT fully supported.
 906  * This arithmetic should never produce a NaN output, but it might
 907  * be confused by a NaN input.
 908  * Define INFINITY in mconf.h for support of infinity; otherwise a
 909  * saturation arithmetic is implemented.
 910  * Denormals are always supported here where appropriate (e.g., not
 911  * for conversion to DEC numbers).
 912  *
 913  */
 914
 915 /*
 916  * Revision history:
 917  *
 918  *  5 Jan 84    PDP-11 assembly language version
 919  *  2 Mar 86    fixed bug in asctoq
 920  *  6 Dec 86    C language version
 921  * 30 Aug 88    100 digit version, improved rounding
 922  * 15 May 92    80-bit long double support
 923  *
 924  * Author:  S. L. Moshier.
 925  */
 926
 927
 928 /*                                                      mconf.h
 929  *
 930  *      Common include file for math routines
 931  *
 932  *
 933  *
 934  * SYNOPSIS:
 935  *
 936  * #include "mconf.h"
 937  *
 938  *
 939  *
 940  * DESCRIPTION:
 941  *
 942  * This file contains definitions for error codes that are
 943  * passed to the common error handling routine mtherr
 944  * (which see).
 945  *
 946  * The file also includes a conditional assembly definition
 947  * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
 948  * IEEE, or UNKnown).
 949  *
 950  * For Digital Equipment PDP-11 and VAX computers, certain
 951  * IBM systems, and others that use numbers with a 56-bit
 952  * significand, the symbol DEC should be defined.  In this
 953  * mode, most floating point constants are given as arrays
 954  * of octal integers to eliminate decimal to binary conversion
 955  * errors that might be introduced by the compiler.
 956  *
 957  * For computers, such as IBM PC, that follow the IEEE
 958  * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
 959  * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
 960  * These numbers have 53-bit significands.  In this mode, constants
 961  * are provided as arrays of hexadecimal 16 bit integers.
 962  *
 963  * To accommodate other types of computer arithmetic, all
 964  * constants are also provided in a normal decimal radix
 965  * which one can hope are correctly converted to a suitable
 966  * format by the available C language compiler.  To invoke
 967  * this mode, the symbol UNK is defined.
 968  *
 969  * An important difference among these modes is a predefined
 970  * set of machine arithmetic constants for each.  The numbers
 971  * MACHEP (the machine roundoff error), MAXNUM (largest number
 972  * represented), and several other parameters are preset by
 973  * the configuration symbol.  Check the file const.c to
 974  * ensure that these values are correct for your computer.
 975  *
 976  * For ANSI C compatibility, define ANSIC equal to 1.  Currently
 977  * this affects only the atan2 function and others that use it.
 978  */
 979 \f
 980 /*
 981 Cephes Math Library Release 2.1:  January, 1989
 982 Copyright 1984, 1987, 1989 by Stephen L. Moshier
 983 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 984 */
 985
 986
 987 /* Constant definitions for math error conditions
 988  */
 989
 990 #define DOMAIN          1       /* argument domain error */
 991 #define SING            2       /* argument singularity */
 992 #define OVERFLOW        3       /* overflow range error */
 993 #define UNDERFLOW       4       /* underflow range error */
 994 #define TLOSS           5       /* total loss of precision */
 995 #define PLOSS           6       /* partial loss of precision */
 996
 997 #define EDOM            33
 998 #define ERANGE          34
 999
1000 /*  e type constants used by high precision check routines */
1001
1002 /*include "ehead.h"*/
1003 /* 0.0 */
1004 unsigned EMUSHORT ezero[NE] =
1005 {
1006   0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1007 extern unsigned EMUSHORT ezero[];
1008
1009 /* 5.0E-1 */
1010 unsigned EMUSHORT ehalf[NE] =
1011 {
1012   0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1013 extern unsigned EMUSHORT ehalf[];
1014
1015 /* 1.0E0 */
1016 unsigned EMUSHORT eone[NE] =
1017 {
1018   0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1019 extern unsigned EMUSHORT eone[];
1020
1021 /* 2.0E0 */
1022 unsigned EMUSHORT etwo[NE] =
1023 {
1024   0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1025 extern unsigned EMUSHORT etwo[];
1026
1027 /* 3.2E1 */
1028 unsigned EMUSHORT e32[NE] =
1029 {
1030   0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1031 extern unsigned EMUSHORT e32[];
1032
1033 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1034 unsigned EMUSHORT elog2[NE] =
1035 {
1036   0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1037 extern unsigned EMUSHORT elog2[];
1038
1039 /* 1.41421356237309504880168872420969807856967187537695E0 */
1040 unsigned EMUSHORT esqrt2[NE] =
1041 {
1042   0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1043 extern unsigned EMUSHORT esqrt2[];
1044
1045 /* 2/sqrt (PI) =
1046  * 1.12837916709551257389615890312154517168810125865800E0 */
1047 unsigned EMUSHORT eoneopi[NE] =
1048 {
1049   0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1050 extern unsigned EMUSHORT eoneopi[];
1051
1052 /* 3.14159265358979323846264338327950288419716939937511E0 */
1053 unsigned EMUSHORT epi[NE] =
1054 {
1055   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1056 extern unsigned EMUSHORT epi[];
1057
1058 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1059 unsigned EMUSHORT eeul[NE] =
1060 {
1061   0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1062 extern unsigned EMUSHORT eeul[];
1063
1064 /*
1065 include "ehead.h"
1066 include "mconf.h"
1067 */
1068
1069
1070
1071 /* Control register for rounding precision.
1072  * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1073  */
1074 int rndprc = NBITS;
1075 extern int rndprc;
1076
1077 void eaddm (), esubm (), emdnorm (), asctoeg ();
1078 static void toe24 (), toe53 (), toe64 ();
1079 void eremain (), einit (), eiremain ();
1080 int ecmpm (), edivm (), emulm ();
1081 void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1082 void etodec (), todec (), dectoe ();
1083
1084
1085
1086
1087 void
1088 einit ()
1089 {
1090 }
1091
1092 /*
1093 ; Clear out entire external format number.
1094 ;
1095 ; unsigned EMUSHORT x[];
1096 ; eclear (x);
1097 */
1098
1099 void
1100 eclear (x)
1101      register unsigned EMUSHORT *x;
1102 {
1103   register int i;
1104
1105   for (i = 0; i < NE; i++)
1106     *x++ = 0;
1107 }
1108
1109
1110
1111 /* Move external format number from a to b.
1112  *
1113  * emov (a, b);
1114  */
1115
1116 void
1117 emov (a, b)
1118      register unsigned EMUSHORT *a, *b;
1119 {
1120   register int i;
1121
1122   for (i = 0; i < NE; i++)
1123     *b++ = *a++;
1124 }
1125
1126
1127 /*
1128 ;       Absolute value of external format number
1129 ;
1130 ;       EMUSHORT x[NE];
1131 ;       eabs (x);
1132 */
1133
1134 void
1135 eabs (x)
1136      unsigned EMUSHORT x[];     /* x is the memory address of a short */
1137 {
1138
1139   x[NE - 1] &= 0x7fff;          /* sign is top bit of last word of external format */
1140 }
1141
1142
1143
1144
1145 /*
1146 ;       Negate external format number
1147 ;
1148 ;       unsigned EMUSHORT x[NE];
1149 ;       eneg (x);
1150 */
1151
1152 void
1153 eneg (x)
1154      unsigned EMUSHORT x[];
1155 {
1156
1157   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1158 }
1159
1160
1161
1162 /* Return 1 if external format number is negative,
1163  * else return zero.
1164  */
1165 int
1166 eisneg (x)
1167      unsigned EMUSHORT x[];
1168 {
1169
1170   if (x[NE - 1] & 0x8000)
1171     return (1);
1172   else
1173     return (0);
1174 }
1175
1176
1177 /* Return 1 if external format number has maximum possible exponent,
1178  * else return zero.
1179  */
1180 int
1181 eisinf (x)
1182      unsigned EMUSHORT x[];
1183 {
1184
1185   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1186     return (1);
1187   else
1188     return (0);
1189 }
1190
1191
1192 /*
1193 ; Fill entire number, including exponent and significand, with
1194 ; largest possible number.  These programs implement a saturation
1195 ; value that is an ordinary, legal number.  A special value
1196 ; "infinity" may also be implemented; this would require tests
1197 ; for that value and implementation of special rules for arithmetic
1198 ; operations involving inifinity.
1199 */
1200
1201 void
1202 einfin (x)
1203      register unsigned EMUSHORT *x;
1204 {
1205   register int i;
1206
1207 #ifdef INFINITY
1208   for (i = 0; i < NE - 1; i++)
1209     *x++ = 0;
1210   *x |= 32767;
1211 #else
1212   for (i = 0; i < NE - 1; i++)
1213     *x++ = 0xffff;
1214   *x |= 32766;
1215   if (rndprc < NBITS)
1216     {
1217       if (rndprc == 64)
1218         {
1219           *(x - 5) = 0;
1220         }
1221       if (rndprc == 53)
1222         {
1223           *(x - 4) = 0xf800;
1224         }
1225       else
1226         {
1227           *(x - 4) = 0;
1228           *(x - 3) = 0;
1229           *(x - 2) = 0xff00;
1230         }
1231     }
1232 #endif
1233 }
1234
1235
1236
1237 /* Move in external format number,
1238  * converting it to internal format.
1239  */
1240 void
1241 emovi (a, b)
1242      unsigned EMUSHORT *a, *b;
1243 {
1244   register unsigned EMUSHORT *p, *q;
1245   int i;
1246
1247   q = b;
1248   p = a + (NE - 1);             /* point to last word of external number */
1249   /* get the sign bit */
1250   if (*p & 0x8000)
1251     *q++ = 0xffff;
1252   else
1253     *q++ = 0;
1254   /* get the exponent */
1255   *q = *p--;
1256   *q++ &= 0x7fff;               /* delete the sign bit */
1257 #ifdef INFINITY
1258   if ((*(q - 1) & 0x7fff) == 0x7fff)
1259     {
1260       for (i = 2; i < NI; i++)
1261         *q++ = 0;
1262       return;
1263     }
1264 #endif
1265   /* clear high guard word */
1266   *q++ = 0;
1267   /* move in the significand */
1268   for (i = 0; i < NE - 1; i++)
1269     *q++ = *p--;
1270   /* clear low guard word */
1271   *q = 0;
1272 }
1273
1274
1275 /* Move internal format number out,
1276  * converting it to external format.
1277  */
1278 void
1279 emovo (a, b)
1280      unsigned EMUSHORT *a, *b;
1281 {
1282   register unsigned EMUSHORT *p, *q;
1283   unsigned EMUSHORT i;
1284
1285   p = a;
1286   q = b + (NE - 1);             /* point to output exponent */
1287   /* combine sign and exponent */
1288   i = *p++;
1289   if (i)
1290     *q-- = *p++ | 0x8000;
1291   else
1292     *q-- = *p++;
1293 #ifdef INFINITY
1294   if (*(p - 1) == 0x7fff)
1295     {
1296       einfin (b);
1297       return;
1298     }
1299 #endif
1300   /* skip over guard word */
1301   ++p;
1302   /* move the significand */
1303   for (i = 0; i < NE - 1; i++)
1304     *q-- = *p++;
1305 }
1306
1307
1308
1309
1310 /* Clear out internal format number.
1311  */
1312
1313 void
1314 ecleaz (xi)
1315      register unsigned EMUSHORT *xi;
1316 {
1317   register int i;
1318
1319   for (i = 0; i < NI; i++)
1320     *xi++ = 0;
1321 }
1322
1323
1324 /* same, but don't touch the sign. */
1325
1326 void
1327 ecleazs (xi)
1328      register unsigned EMUSHORT *xi;
1329 {
1330   register int i;
1331
1332   ++xi;
1333   for (i = 0; i < NI - 1; i++)
1334     *xi++ = 0;
1335 }
1336
1337
1338
1339 /* Move internal format number from a to b.
1340  */
1341 void
1342 emovz (a, b)
1343      register unsigned EMUSHORT *a, *b;
1344 {
1345   register int i;
1346
1347   for (i = 0; i < NI - 1; i++)
1348     *b++ = *a++;
1349   /* clear low guard word */
1350   *b = 0;
1351 }
1352
1353
1354 /*
1355 ;       Compare significands of numbers in internal format.
1356 ;       Guard words are included in the comparison.
1357 ;
1358 ;       unsigned EMUSHORT a[NI], b[NI];
1359 ;       cmpm (a, b);
1360 ;
1361 ;       for the significands:
1362 ;       returns +1 if a > b
1363 ;                0 if a == b
1364 ;               -1 if a < b
1365 */
1366 int
1367 ecmpm (a, b)
1368      register unsigned EMUSHORT *a, *b;
1369 {
1370   int i;
1371
1372   a += M;                       /* skip up to significand area */
1373   b += M;
1374   for (i = M; i < NI; i++)
1375     {
1376       if (*a++ != *b++)
1377         goto difrnt;
1378     }
1379   return (0);
1380
1381  difrnt:
1382   if (*(--a) > *(--b))
1383     return (1);
1384   else
1385     return (-1);
1386 }
1387
1388
1389 /*
1390 ;       Shift significand down by 1 bit
1391 */
1392
1393 void
1394 eshdn1 (x)
1395      register unsigned EMUSHORT *x;
1396 {
1397   register unsigned EMUSHORT bits;
1398   int i;
1399
1400   x += M;                       /* point to significand area */
1401
1402   bits = 0;
1403   for (i = M; i < NI; i++)
1404     {
1405       if (*x & 1)
1406         bits |= 1;
1407       *x >>= 1;
1408       if (bits & 2)
1409         *x |= 0x8000;
1410       bits <<= 1;
1411       ++x;
1412     }
1413 }
1414
1415
1416
1417 /*
1418 ;       Shift significand up by 1 bit
1419 */
1420
1421 void
1422 eshup1 (x)
1423      register unsigned EMUSHORT *x;
1424 {
1425   register unsigned EMUSHORT bits;
1426   int i;
1427
1428   x += NI - 1;
1429   bits = 0;
1430
1431   for (i = M; i < NI; i++)
1432     {
1433       if (*x & 0x8000)
1434         bits |= 1;
1435       *x <<= 1;
1436       if (bits & 2)
1437         *x |= 1;
1438       bits <<= 1;
1439       --x;
1440     }
1441 }
1442
1443
1444
1445 /*
1446 ;       Shift significand down by 8 bits
1447 */
1448
1449 void
1450 eshdn8 (x)
1451      register unsigned EMUSHORT *x;
1452 {
1453   register unsigned EMUSHORT newbyt, oldbyt;
1454   int i;
1455
1456   x += M;
1457   oldbyt = 0;
1458   for (i = M; i < NI; i++)
1459     {
1460       newbyt = *x << 8;
1461       *x >>= 8;
1462       *x |= oldbyt;
1463       oldbyt = newbyt;
1464       ++x;
1465     }
1466 }
1467
1468 /*
1469 ;       Shift significand up by 8 bits
1470 */
1471
1472 void
1473 eshup8 (x)
1474      register unsigned EMUSHORT *x;
1475 {
1476   int i;
1477   register unsigned EMUSHORT newbyt, oldbyt;
1478
1479   x += NI - 1;
1480   oldbyt = 0;
1481
1482   for (i = M; i < NI; i++)
1483     {
1484       newbyt = *x >> 8;
1485       *x <<= 8;
1486       *x |= oldbyt;
1487       oldbyt = newbyt;
1488       --x;
1489     }
1490 }
1491
1492 /*
1493 ;       Shift significand up by 16 bits
1494 */
1495
1496 void
1497 eshup6 (x)
1498      register unsigned EMUSHORT *x;
1499 {
1500   int i;
1501   register unsigned EMUSHORT *p;
1502
1503   p = x + M;
1504   x += M + 1;
1505
1506   for (i = M; i < NI - 1; i++)
1507     *p++ = *x++;
1508
1509   *p = 0;
1510 }
1511
1512 /*
1513 ;       Shift significand down by 16 bits
1514 */
1515
1516 void
1517 eshdn6 (x)
1518      register unsigned EMUSHORT *x;
1519 {
1520   int i;
1521   register unsigned EMUSHORT *p;
1522
1523   x += NI - 1;
1524   p = x + 1;
1525
1526   for (i = M; i < NI - 1; i++)
1527     *(--p) = *(--x);
1528
1529   *(--p) = 0;
1530 }
1531 \f
1532 /*
1533 ;       Add significands
1534 ;       x + y replaces y
1535 */
1536
1537 void
1538 eaddm (x, y)
1539      unsigned EMUSHORT *x, *y;
1540 {
1541   register unsigned EMULONG a;
1542   int i;
1543   unsigned int carry;
1544
1545   x += NI - 1;
1546   y += NI - 1;
1547   carry = 0;
1548   for (i = M; i < NI; i++)
1549     {
1550       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1551       if (a & 0x10000)
1552         carry = 1;
1553       else
1554         carry = 0;
1555       *y = (unsigned EMUSHORT) a;
1556       --x;
1557       --y;
1558     }
1559 }
1560
1561 /*
1562 ;       Subtract significands
1563 ;       y - x replaces y
1564 */
1565
1566 void
1567 esubm (x, y)
1568      unsigned EMUSHORT *x, *y;
1569 {
1570   unsigned EMULONG a;
1571   int i;
1572   unsigned int carry;
1573
1574   x += NI - 1;
1575   y += NI - 1;
1576   carry = 0;
1577   for (i = M; i < NI; i++)
1578     {
1579       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1580       if (a & 0x10000)
1581         carry = 1;
1582       else
1583         carry = 0;
1584       *y = (unsigned EMUSHORT) a;
1585       --x;
1586       --y;
1587     }
1588 }
1589
1590
1591 /* Divide significands */
1592
1593 static unsigned EMUSHORT equot[NI];
1594
1595 int
1596 edivm (den, num)
1597      unsigned EMUSHORT den[], num[];
1598 {
1599   int i;
1600   register unsigned EMUSHORT *p, *q;
1601   unsigned EMUSHORT j;
1602
1603   p = &equot[0];
1604   *p++ = num[0];
1605   *p++ = num[1];
1606
1607   for (i = M; i < NI; i++)
1608     {
1609       *p++ = 0;
1610     }
1611
1612   /* Use faster compare and subtraction if denominator
1613    * has only 15 bits of significance.
1614    */
1615   p = &den[M + 2];
1616   if (*p++ == 0)
1617     {
1618       for (i = M + 3; i < NI; i++)
1619         {
1620           if (*p++ != 0)
1621             goto fulldiv;
1622         }
1623       if ((den[M + 1] & 1) != 0)
1624         goto fulldiv;
1625       eshdn1 (num);
1626       eshdn1 (den);
1627
1628       p = &den[M + 1];
1629       q = &num[M + 1];
1630
1631       for (i = 0; i < NBITS + 2; i++)
1632         {
1633           if (*p <= *q)
1634             {
1635               *q -= *p;
1636               j = 1;
1637             }
1638           else
1639             {
1640               j = 0;
1641             }
1642           eshup1 (equot);
1643           equot[NI - 2] |= j;
1644           eshup1 (num);
1645         }
1646       goto divdon;
1647     }
1648
1649   /* The number of quotient bits to calculate is
1650    * NBITS + 1 scaling guard bit + 1 roundoff bit.
1651    */
1652  fulldiv:
1653
1654   p = &equot[NI - 2];
1655   for (i = 0; i < NBITS + 2; i++)
1656     {
1657       if (ecmpm (den, num) <= 0)
1658         {
1659           esubm (den, num);
1660           j = 1;                /* quotient bit = 1 */
1661         }
1662       else
1663         j = 0;
1664       eshup1 (equot);
1665       *p |= j;
1666       eshup1 (num);
1667     }
1668
1669  divdon:
1670
1671   eshdn1 (equot);
1672   eshdn1 (equot);
1673
1674   /* test for nonzero remainder after roundoff bit */
1675   p = &num[M];
1676   j = 0;
1677   for (i = M; i < NI; i++)
1678     {
1679       j |= *p++;
1680     }
1681   if (j)
1682     j = 1;
1683
1684
1685   for (i = 0; i < NI; i++)
1686     num[i] = equot[i];
1687   return ((int) j);
1688 }
1689
1690
1691 /* Multiply significands */
1692 int
1693 emulm (a, b)
1694      unsigned EMUSHORT a[], b[];
1695 {
1696   unsigned EMUSHORT *p, *q;
1697   int i, j, k;
1698
1699   equot[0] = b[0];
1700   equot[1] = b[1];
1701   for (i = M; i < NI; i++)
1702     equot[i] = 0;
1703
1704   p = &a[NI - 2];
1705   k = NBITS;
1706   while (*p == 0)               /* significand is not supposed to be all zero */
1707     {
1708       eshdn6 (a);
1709       k -= 16;
1710     }
1711   if ((*p & 0xff) == 0)
1712     {
1713       eshdn8 (a);
1714       k -= 8;
1715     }
1716
1717   q = &equot[NI - 1];
1718   j = 0;
1719   for (i = 0; i < k; i++)
1720     {
1721       if (*p & 1)
1722         eaddm (b, equot);
1723       /* remember if there were any nonzero bits shifted out */
1724       if (*q & 1)
1725         j |= 1;
1726       eshdn1 (a);
1727       eshdn1 (equot);
1728     }
1729
1730   for (i = 0; i < NI; i++)
1731     b[i] = equot[i];
1732
1733   /* return flag for lost nonzero bits */
1734   return (j);
1735 }
1736
1737
1738
1739 /*
1740  * Normalize and round off.
1741  *
1742  * The internal format number to be rounded is "s".
1743  * Input "lost" indicates whether or not the number is exact.
1744  * This is the so-called sticky bit.
1745  *
1746  * Input "subflg" indicates whether the number was obtained
1747  * by a subtraction operation.  In that case if lost is nonzero
1748  * then the number is slightly smaller than indicated.
1749  *
1750  * Input "exp" is the biased exponent, which may be negative.
1751  * the exponent field of "s" is ignored but is replaced by
1752  * "exp" as adjusted by normalization and rounding.
1753  *
1754  * Input "rcntrl" is the rounding control.
1755  */
1756
1757 static int rlast = -1;
1758 static int rw = 0;
1759 static unsigned EMUSHORT rmsk = 0;
1760 static unsigned EMUSHORT rmbit = 0;
1761 static unsigned EMUSHORT rebit = 0;
1762 static int re = 0;
1763 static unsigned EMUSHORT rbit[NI];
1764
1765 void
1766 emdnorm (s, lost, subflg, exp, rcntrl)
1767      unsigned EMUSHORT s[];
1768      int lost;
1769      int subflg;
1770      EMULONG exp;
1771      int rcntrl;
1772 {
1773   int i, j;
1774   unsigned EMUSHORT r;
1775
1776   /* Normalize */
1777   j = enormlz (s);
1778
1779   /* a blank significand could mean either zero or infinity. */
1780 #ifndef INFINITY
1781   if (j > NBITS)
1782     {
1783       ecleazs (s);
1784       return;
1785     }
1786 #endif
1787   exp -= j;
1788 #ifndef INFINITY
1789   if (exp >= 32767L)
1790     goto overf;
1791 #else
1792   if ((j > NBITS) && (exp < 32767))
1793     {
1794       ecleazs (s);
1795       return;
1796     }
1797 #endif
1798   if (exp < 0L)
1799     {
1800       if (exp > (EMULONG) (-NBITS - 1))
1801         {
1802           j = (int) exp;
1803           i = eshift (s, j);
1804           if (i)
1805             lost = 1;
1806         }
1807       else
1808         {
1809           ecleazs (s);
1810           return;
1811         }
1812     }
1813   /* Round off, unless told not to by rcntrl. */
1814   if (rcntrl == 0)
1815     goto mdfin;
1816   /* Set up rounding parameters if the control register changed. */
1817   if (rndprc != rlast)
1818     {
1819       ecleaz (rbit);
1820       switch (rndprc)
1821         {
1822         default:
1823         case NBITS:
1824           rw = NI - 1;          /* low guard word */
1825           rmsk = 0xffff;
1826           rmbit = 0x8000;
1827           rbit[rw - 1] = 1;
1828           re = NI - 2;
1829           rebit = 1;
1830           break;
1831         case 64:
1832           rw = 7;
1833           rmsk = 0xffff;
1834           rmbit = 0x8000;
1835           rbit[rw - 1] = 1;
1836           re = rw - 1;
1837           rebit = 1;
1838           break;
1839           /* For DEC arithmetic */
1840         case 56:
1841           rw = 6;
1842           rmsk = 0xff;
1843           rmbit = 0x80;
1844           rbit[rw] = 0x100;
1845           re = rw;
1846           rebit = 0x100;
1847           break;
1848         case 53:
1849           rw = 6;
1850           rmsk = 0x7ff;
1851           rmbit = 0x0400;
1852           rbit[rw] = 0x800;
1853           re = rw;
1854           rebit = 0x800;
1855           break;
1856         case 24:
1857           rw = 4;
1858           rmsk = 0xff;
1859           rmbit = 0x80;
1860           rbit[rw] = 0x100;
1861           re = rw;
1862           rebit = 0x100;
1863           break;
1864         }
1865       rlast = rndprc;
1866     }
1867
1868   if (rndprc >= 64)
1869     {
1870       r = s[rw] & rmsk;
1871       if (rndprc == 64)
1872         {
1873           i = rw + 1;
1874           while (i < NI)
1875             {
1876               if (s[i])
1877                 r |= 1;
1878               s[i] = 0;
1879               ++i;
1880             }
1881         }
1882     }
1883   else
1884     {
1885       if (exp <= 0)
1886         eshdn1 (s);
1887       r = s[rw] & rmsk;
1888       /* These tests assume NI = 8 */
1889       i = rw + 1;
1890       while (i < NI)
1891         {
1892           if (s[i])
1893             r |= 1;
1894           s[i] = 0;
1895           ++i;
1896         }
1897       /*
1898          if (rndprc == 24)
1899          {
1900          if (s[5] || s[6])
1901          r |= 1;
1902          s[5] = 0;
1903          s[6] = 0;
1904          }
1905          */
1906       s[rw] &= ~rmsk;
1907     }
1908   if ((r & rmbit) != 0)
1909     {
1910       if (r == rmbit)
1911         {
1912           if (lost == 0)
1913             {                   /* round to even */
1914               if ((s[re] & rebit) == 0)
1915                 goto mddone;
1916             }
1917           else
1918             {
1919               if (subflg != 0)
1920                 goto mddone;
1921             }
1922         }
1923       eaddm (rbit, s);
1924     }
1925  mddone:
1926   if ((rndprc < 64) && (exp <= 0))
1927     {
1928       eshup1 (s);
1929     }
1930   if (s[2] != 0)
1931     {                           /* overflow on roundoff */
1932       eshdn1 (s);
1933       exp += 1;
1934     }
1935  mdfin:
1936   s[NI - 1] = 0;
1937   if (exp >= 32767L)
1938     {
1939 #ifndef INFINITY
1940     overf:
1941 #endif
1942 #ifdef INFINITY
1943       s[1] = 32767;
1944       for (i = 2; i < NI - 1; i++)
1945         s[i] = 0;
1946       pedwarn ("floating point overflow");
1947 #else
1948       s[1] = 32766;
1949       s[2] = 0;
1950       for (i = M + 1; i < NI - 1; i++)
1951         s[i] = 0xffff;
1952       s[NI - 1] = 0;
1953       if (rndprc < 64)
1954         {
1955           s[rw] &= ~rmsk;
1956           if (rndprc == 24)
1957             {
1958               s[5] = 0;
1959               s[6] = 0;
1960             }
1961         }
1962 #endif
1963       return;
1964     }
1965   if (exp < 0)
1966     s[1] = 0;
1967   else
1968     s[1] = (unsigned EMUSHORT) exp;
1969 }
1970
1971
1972
1973 /*
1974 ;       Subtract external format numbers.
1975 ;
1976 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1977 ;       esub (a, b, c);  c = b - a
1978 */
1979
1980 static int subflg = 0;
1981
1982 void
1983 esub (a, b, c)
1984      unsigned EMUSHORT *a, *b, *c;
1985 {
1986
1987   subflg = 1;
1988   eadd1 (a, b, c);
1989 }
1990
1991
1992 /*
1993 ;       Add.
1994 ;
1995 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1996 ;       eadd (a, b, c);  c = b + a
1997 */
1998 void
1999 eadd (a, b, c)
2000      unsigned EMUSHORT *a, *b, *c;
2001 {
2002
2003   subflg = 0;
2004   eadd1 (a, b, c);
2005 }
2006
2007 void
2008 eadd1 (a, b, c)
2009      unsigned EMUSHORT *a, *b, *c;
2010 {
2011   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2012   int i, lost, j, k;
2013   EMULONG lt, lta, ltb;
2014
2015 #ifdef INFINITY
2016   if (eisinf (a))
2017     {
2018       emov (a, c);
2019       if (subflg)
2020         eneg (c);
2021       return;
2022     }
2023   if (eisinf (b))
2024     {
2025       emov (b, c);
2026       return;
2027     }
2028 #endif
2029   emovi (a, ai);
2030   emovi (b, bi);
2031   if (subflg)
2032     ai[0] = ~ai[0];
2033
2034   /* compare exponents */
2035   lta = ai[E];
2036   ltb = bi[E];
2037   lt = lta - ltb;
2038   if (lt > 0L)
2039     {                           /* put the larger number in bi */
2040       emovz (bi, ci);
2041       emovz (ai, bi);
2042       emovz (ci, ai);
2043       ltb = bi[E];
2044       lt = -lt;
2045     }
2046   lost = 0;
2047   if (lt != 0L)
2048     {
2049       if (lt < (EMULONG) (-NBITS - 1))
2050         goto done;              /* answer same as larger addend */
2051       k = (int) lt;
2052       lost = eshift (ai, k);    /* shift the smaller number down */
2053     }
2054   else
2055     {
2056       /* exponents were the same, so must compare significands */
2057       i = ecmpm (ai, bi);
2058       if (i == 0)
2059         {                       /* the numbers are identical in magnitude */
2060           /* if different signs, result is zero */
2061           if (ai[0] != bi[0])
2062             {
2063               eclear (c);
2064               return;
2065             }
2066           /* if same sign, result is double */
2067           /* double denomalized tiny number */
2068           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2069             {
2070               eshup1 (bi);
2071               goto done;
2072             }
2073           /* add 1 to exponent unless both are zero! */
2074           for (j = 1; j < NI - 1; j++)
2075             {
2076               if (bi[j] != 0)
2077                 {
2078                   /* This could overflow, but let emovo take care of that. */
2079                   ltb += 1;
2080                   break;
2081                 }
2082             }
2083           bi[E] = (unsigned EMUSHORT) ltb;
2084           goto done;
2085         }
2086       if (i > 0)
2087         {                       /* put the larger number in bi */
2088           emovz (bi, ci);
2089           emovz (ai, bi);
2090           emovz (ci, ai);
2091         }
2092     }
2093   if (ai[0] == bi[0])
2094     {
2095       eaddm (ai, bi);
2096       subflg = 0;
2097     }
2098   else
2099     {
2100       esubm (ai, bi);
2101       subflg = 1;
2102     }
2103   emdnorm (bi, lost, subflg, ltb, 64);
2104
2105  done:
2106   emovo (bi, c);
2107 }
2108
2109
2110
2111 /*
2112 ;       Divide.
2113 ;
2114 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2115 ;       ediv (a, b, c); c = b / a
2116 */
2117 void
2118 ediv (a, b, c)
2119      unsigned EMUSHORT *a, *b, *c;
2120 {
2121   unsigned EMUSHORT ai[NI], bi[NI];
2122   int i;
2123   EMULONG lt, lta, ltb;
2124
2125 #ifdef INFINITY
2126   if (eisinf (b))
2127     {
2128       if (eisneg (a) ^ eisneg (b))
2129         *(c + (NE - 1)) = 0x8000;
2130       else
2131         *(c + (NE - 1)) = 0;
2132       einfin (c);
2133       return;
2134     }
2135   if (eisinf (a))
2136     {
2137       eclear (c);
2138       return;
2139     }
2140 #endif
2141   emovi (a, ai);
2142   emovi (b, bi);
2143   lta = ai[E];
2144   ltb = bi[E];
2145   if (bi[E] == 0)
2146     {                           /* See if numerator is zero. */
2147       for (i = 1; i < NI - 1; i++)
2148         {
2149           if (bi[i] != 0)
2150             {
2151               ltb -= enormlz (bi);
2152               goto dnzro1;
2153             }
2154         }
2155       eclear (c);
2156       return;
2157     }
2158  dnzro1:
2159
2160   if (ai[E] == 0)
2161     {                           /* possible divide by zero */
2162       for (i = 1; i < NI - 1; i++)
2163         {
2164           if (ai[i] != 0)
2165             {
2166               lta -= enormlz (ai);
2167               goto dnzro2;
2168             }
2169         }
2170       if (ai[0] == bi[0])
2171         *(c + (NE - 1)) = 0;
2172       else
2173         *(c + (NE - 1)) = 0x8000;
2174       einfin (c);
2175       mtherr ("ediv", SING);
2176       return;
2177     }
2178  dnzro2:
2179
2180   i = edivm (ai, bi);
2181   /* calculate exponent */
2182   lt = ltb - lta + EXONE;
2183   emdnorm (bi, i, 0, lt, 64);
2184   /* set the sign */
2185   if (ai[0] == bi[0])
2186     bi[0] = 0;
2187   else
2188     bi[0] = 0Xffff;
2189   emovo (bi, c);
2190 }
2191
2192
2193
2194 /*
2195 ;       Multiply.
2196 ;
2197 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2198 ;       emul (a, b, c); c = b * a
2199 */
2200 void
2201 emul (a, b, c)
2202      unsigned EMUSHORT *a, *b, *c;
2203 {
2204   unsigned EMUSHORT ai[NI], bi[NI];
2205   int i, j;
2206   EMULONG lt, lta, ltb;
2207
2208 #ifdef INFINITY
2209   if (eisinf (a) || eisinf (b))
2210     {
2211       if (eisneg (a) ^ eisneg (b))
2212         *(c + (NE - 1)) = 0x8000;
2213       else
2214         *(c + (NE - 1)) = 0;
2215       einfin (c);
2216       return;
2217     }
2218 #endif
2219   emovi (a, ai);
2220   emovi (b, bi);
2221   lta = ai[E];
2222   ltb = bi[E];
2223   if (ai[E] == 0)
2224     {
2225       for (i = 1; i < NI - 1; i++)
2226         {
2227           if (ai[i] != 0)
2228             {
2229               lta -= enormlz (ai);
2230               goto mnzer1;
2231             }
2232         }
2233       eclear (c);
2234       return;
2235     }
2236  mnzer1:
2237
2238   if (bi[E] == 0)
2239     {
2240       for (i = 1; i < NI - 1; i++)
2241         {
2242           if (bi[i] != 0)
2243             {
2244               ltb -= enormlz (bi);
2245               goto mnzer2;
2246             }
2247         }
2248       eclear (c);
2249       return;
2250     }
2251  mnzer2:
2252
2253   /* Multiply significands */
2254   j = emulm (ai, bi);
2255   /* calculate exponent */
2256   lt = lta + ltb - (EXONE - 1);
2257   emdnorm (bi, j, 0, lt, 64);
2258   /* calculate sign of product */
2259   if (ai[0] == bi[0])
2260     bi[0] = 0;
2261   else
2262     bi[0] = 0xffff;
2263   emovo (bi, c);
2264 }
2265
2266
2267
2268
2269 /*
2270 ; Convert IEEE double precision to e type
2271 ;       double d;
2272 ;       unsigned EMUSHORT x[N+2];
2273 ;       e53toe (&d, x);
2274 */
2275 void
2276 e53toe (e, y)
2277      unsigned EMUSHORT *e, *y;
2278 {
2279 #ifdef DEC
2280
2281   dectoe (e, y);                /* see etodec.c */
2282
2283 #else
2284
2285   register unsigned EMUSHORT r;
2286   register unsigned EMUSHORT *p;
2287   unsigned EMUSHORT yy[NI];
2288   int denorm, k;
2289
2290   denorm = 0;                   /* flag if denormalized number */
2291   ecleaz (yy);
2292 #ifdef IBMPC
2293   e += 3;
2294 #endif
2295   r = *e;
2296   yy[0] = 0;
2297   if (r & 0x8000)
2298     yy[0] = 0xffff;
2299   yy[M] = (r & 0x0f) | 0x10;
2300   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2301 #ifdef INFINITY
2302   if (r == 0x7ff0)
2303     {
2304       einfin (y);
2305       if (r & 0x8000)
2306         eneg (y);
2307       return;
2308     }
2309 #endif
2310   r >>= 4;
2311   /* If zero exponent, then the significand is denormalized.
2312    * So, take back the understood high significand bit. */
2313   if (r == 0)
2314     {
2315       denorm = 1;
2316       yy[M] &= ~0x10;
2317     }
2318   r += EXONE - 01777;
2319   yy[E] = r;
2320   p = &yy[M + 1];
2321 #ifdef IBMPC
2322   *p++ = *(--e);
2323   *p++ = *(--e);
2324   *p++ = *(--e);
2325 #endif
2326 #ifdef MIEEE
2327   ++e;
2328   *p++ = *e++;
2329   *p++ = *e++;
2330   *p++ = *e++;
2331 #endif
2332   (void) eshift (yy, -5);
2333   if (denorm)
2334     {                           /* if zero exponent, then normalize the significand */
2335       if ((k = enormlz (yy)) > NBITS)
2336         ecleazs (yy);
2337       else
2338         yy[E] -= (unsigned EMUSHORT) (k - 1);
2339     }
2340   emovo (yy, y);
2341 #endif /* not DEC */
2342 }
2343
2344 void
2345 e64toe (e, y)
2346      unsigned EMUSHORT *e, *y;
2347 {
2348   unsigned EMUSHORT yy[NI];
2349   unsigned EMUSHORT *p, *q;
2350   int i;
2351
2352   p = yy;
2353   for (i = 0; i < NE - 5; i++)
2354     *p++ = 0;
2355 #ifdef IBMPC
2356   for (i = 0; i < 5; i++)
2357     *p++ = *e++;
2358 #endif
2359 #ifdef DEC
2360   for (i = 0; i < 5; i++)
2361     *p++ = *e++;
2362 #endif
2363 #ifdef MIEEE
2364   p = &yy[0] + (NE - 1);
2365   *p-- = *e++;
2366   ++e;
2367   for (i = 0; i < 4; i++)
2368     *p-- = *e++;
2369 #endif
2370   p = yy;
2371   q = y;
2372 #ifdef INFINITY
2373   if (*p == 0x7fff)
2374     {
2375       einfin (y);
2376       if (*p & 0x8000)
2377         eneg (y);
2378       return;
2379     }
2380 #endif
2381   for (i = 0; i < NE; i++)
2382     *q++ = *p++;
2383 }
2384
2385
2386 /*
2387 ; Convert IEEE single precision to e type
2388 ;       float d;
2389 ;       unsigned EMUSHORT x[N+2];
2390 ;       dtox (&d, x);
2391 */
2392 void
2393 e24toe (e, y)
2394      unsigned EMUSHORT *e, *y;
2395 {
2396   register unsigned EMUSHORT r;
2397   register unsigned EMUSHORT *p;
2398   unsigned EMUSHORT yy[NI];
2399   int denorm, k;
2400
2401   denorm = 0;                   /* flag if denormalized number */
2402   ecleaz (yy);
2403 #ifdef IBMPC
2404   e += 1;
2405 #endif
2406 #ifdef DEC
2407   e += 1;
2408 #endif
2409   r = *e;
2410   yy[0] = 0;
2411   if (r & 0x8000)
2412     yy[0] = 0xffff;
2413   yy[M] = (r & 0x7f) | 0200;
2414   r &= ~0x807f;                 /* strip sign and 7 significand bits */
2415 #ifdef INFINITY
2416   if (r == 0x7f80)
2417     {
2418       einfin (y);
2419       if (r & 0x8000)
2420         eneg (y);
2421       return;
2422     }
2423 #endif
2424   r >>= 7;
2425   /* If zero exponent, then the significand is denormalized.
2426    * So, take back the understood high significand bit. */
2427   if (r == 0)
2428     {
2429       denorm = 1;
2430       yy[M] &= ~0200;
2431     }
2432   r += EXONE - 0177;
2433   yy[E] = r;
2434   p = &yy[M + 1];
2435 #ifdef IBMPC
2436   *p++ = *(--e);
2437 #endif
2438 #ifdef DEC
2439   *p++ = *(--e);
2440 #endif
2441 #ifdef MIEEE
2442   ++e;
2443   *p++ = *e++;
2444 #endif
2445   (void) eshift (yy, -8);
2446   if (denorm)
2447     {                           /* if zero exponent, then normalize the significand */
2448       if ((k = enormlz (yy)) > NBITS)
2449         ecleazs (yy);
2450       else
2451         yy[E] -= (unsigned EMUSHORT) (k - 1);
2452     }
2453   emovo (yy, y);
2454 }
2455
2456
2457 void
2458 etoe64 (x, e)
2459      unsigned EMUSHORT *x, *e;
2460 {
2461   unsigned EMUSHORT xi[NI];
2462   EMULONG exp;
2463   int rndsav;
2464
2465   emovi (x, xi);
2466   /* adjust exponent for offset */
2467   exp = (EMULONG) xi[E];
2468 #ifdef INFINITY
2469   if (eisinf (x))
2470     goto nonorm;
2471 #endif
2472   /* round off to nearest or even */
2473   rndsav = rndprc;
2474   rndprc = 64;
2475   emdnorm (xi, 0, 0, exp, 64);
2476   rndprc = rndsav;
2477  nonorm:
2478   toe64 (xi, e);
2479 }
2480
2481 /* move out internal format to ieee long double */
2482 static void
2483 toe64 (a, b)
2484      unsigned EMUSHORT *a, *b;
2485 {
2486   register unsigned EMUSHORT *p, *q;
2487   unsigned EMUSHORT i;
2488
2489   p = a;
2490 #ifdef MIEEE
2491   q = b;
2492 #else
2493   q = b + 4;                    /* point to output exponent */
2494 #if LONG_DOUBLE_TYPE_SIZE == 96
2495   /* Clear the last two bytes of 12-byte Intel format */
2496   *(q+1) = 0;
2497 #endif
2498 #endif
2499
2500   /* combine sign and exponent */
2501   i = *p++;
2502 #ifdef MIEEE
2503   if (i)
2504     *q++ = *p++ | 0x8000;
2505   else
2506     *q++ = *p++;
2507   *q++ = 0;
2508 #else
2509   if (i)
2510     *q-- = *p++ | 0x8000;
2511   else
2512     *q-- = *p++;
2513 #endif
2514   /* skip over guard word */
2515   ++p;
2516   /* move the significand */
2517 #ifdef MIEEE
2518   for (i = 0; i < 4; i++)
2519     *q++ = *p++;
2520 #else
2521   for (i = 0; i < 4; i++)
2522     *q-- = *p++;
2523 #endif
2524 }
2525
2526
2527 /*
2528 ; e type to IEEE double precision
2529 ;       double d;
2530 ;       unsigned EMUSHORT x[NE];
2531 ;       etoe53 (x, &d);
2532 */
2533
2534 #ifdef DEC
2535
2536 void
2537 etoe53 (x, e)
2538      unsigned EMUSHORT *x, *e;
2539 {
2540   etodec (x, e);                /* see etodec.c */
2541 }
2542
2543 static void
2544 toe53 (x, y)
2545      unsigned EMUSHORT *x, *y;
2546 {
2547   todec (x, y);
2548 }
2549
2550 #else
2551
2552 void
2553 etoe53 (x, e)
2554      unsigned EMUSHORT *x, *e;
2555 {
2556   unsigned EMUSHORT xi[NI];
2557   EMULONG exp;
2558   int rndsav;
2559
2560   emovi (x, xi);
2561   /* adjust exponent for offsets */
2562   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2563 #ifdef INFINITY
2564   if (eisinf (x))
2565     goto nonorm;
2566 #endif
2567   /* round off to nearest or even */
2568   rndsav = rndprc;
2569   rndprc = 53;
2570   emdnorm (xi, 0, 0, exp, 64);
2571   rndprc = rndsav;
2572  nonorm:
2573   toe53 (xi, e);
2574 }
2575
2576
2577 static void
2578 toe53 (x, y)
2579      unsigned EMUSHORT *x, *y;
2580 {
2581   unsigned EMUSHORT i;
2582   unsigned EMUSHORT *p;
2583
2584
2585   p = &x[0];
2586 #ifdef IBMPC
2587   y += 3;
2588 #endif
2589   *y = 0;                       /* output high order */
2590   if (*p++)
2591     *y = 0x8000;                /* output sign bit */
2592
2593   i = *p++;
2594   if (i >= (unsigned int) 2047)
2595     {                           /* Saturate at largest number less than infinity. */
2596 #ifdef INFINITY
2597       *y |= 0x7ff0;
2598 #ifdef IBMPC
2599       *(--y) = 0;
2600       *(--y) = 0;
2601       *(--y) = 0;
2602 #endif
2603 #ifdef MIEEE
2604       ++y;
2605       *y++ = 0;
2606       *y++ = 0;
2607       *y++ = 0;
2608 #endif
2609 #else
2610       *y |= (unsigned EMUSHORT) 0x7fef;
2611 #ifdef IBMPC
2612       *(--y) = 0xffff;
2613       *(--y) = 0xffff;
2614       *(--y) = 0xffff;
2615 #endif
2616 #ifdef MIEEE
2617       ++y;
2618       *y++ = 0xffff;
2619       *y++ = 0xffff;
2620       *y++ = 0xffff;
2621 #endif
2622 #endif
2623       return;
2624     }
2625   if (i == 0)
2626     {
2627       (void) eshift (x, 4);
2628     }
2629   else
2630     {
2631       i <<= 4;
2632       (void) eshift (x, 5);
2633     }
2634   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
2635   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
2636 #ifdef IBMPC
2637   *(--y) = *p++;
2638   *(--y) = *p++;
2639   *(--y) = *p;
2640 #endif
2641 #ifdef MIEEE
2642   ++y;
2643   *y++ = *p++;
2644   *y++ = *p++;
2645   *y++ = *p++;
2646 #endif
2647 }
2648
2649 #endif /* not DEC */
2650
2651
2652
2653 /*
2654 ; e type to IEEE single precision
2655 ;       float d;
2656 ;       unsigned EMUSHORT x[N+2];
2657 ;       xtod (x, &d);
2658 */
2659 void
2660 etoe24 (x, e)
2661      unsigned EMUSHORT *x, *e;
2662 {
2663   EMULONG exp;
2664   unsigned EMUSHORT xi[NI];
2665   int rndsav;
2666
2667   emovi (x, xi);
2668   /* adjust exponent for offsets */
2669   exp = (EMULONG) xi[E] - (EXONE - 0177);
2670 #ifdef INFINITY
2671   if (eisinf (x))
2672     goto nonorm;
2673 #endif
2674   /* round off to nearest or even */
2675   rndsav = rndprc;
2676   rndprc = 24;
2677   emdnorm (xi, 0, 0, exp, 64);
2678   rndprc = rndsav;
2679  nonorm:
2680   toe24 (xi, e);
2681 }
2682
2683 static void
2684 toe24 (x, y)
2685      unsigned EMUSHORT *x, *y;
2686 {
2687   unsigned EMUSHORT i;
2688   unsigned EMUSHORT *p;
2689
2690   p = &x[0];
2691 #ifdef IBMPC
2692   y += 1;
2693 #endif
2694 #ifdef DEC
2695   y += 1;
2696 #endif
2697   *y = 0;                       /* output high order */
2698   if (*p++)
2699     *y = 0x8000;                /* output sign bit */
2700
2701   i = *p++;
2702   if (i >= 255)
2703     {                           /* Saturate at largest number less than infinity. */
2704 #ifdef INFINITY
2705       *y |= (unsigned EMUSHORT) 0x7f80;
2706 #ifdef IBMPC
2707       *(--y) = 0;
2708 #endif
2709 #ifdef DEC
2710       *(--y) = 0;
2711 #endif
2712 #ifdef MIEEE
2713       ++y;
2714       *y = 0;
2715 #endif
2716 #else
2717       *y |= (unsigned EMUSHORT) 0x7f7f;
2718 #ifdef IBMPC
2719       *(--y) = 0xffff;
2720 #endif
2721 #ifdef DEC
2722       *(--y) = 0xffff;
2723 #endif
2724 #ifdef MIEEE
2725       ++y;
2726       *y = 0xffff;
2727 #endif
2728 #endif
2729       return;
2730     }
2731   if (i == 0)
2732     {
2733       (void) eshift (x, 7);
2734     }
2735   else
2736     {
2737       i <<= 7;
2738       (void) eshift (x, 8);
2739     }
2740   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
2741   *y |= i;                      /* high order output already has sign bit set */
2742 #ifdef IBMPC
2743   *(--y) = *p;
2744 #endif
2745 #ifdef DEC
2746   *(--y) = *p;
2747 #endif
2748 #ifdef MIEEE
2749   ++y;
2750   *y = *p;
2751 #endif
2752 }
2753
2754
2755 /* Compare two e type numbers.
2756  *
2757  * unsigned EMUSHORT a[NE], b[NE];
2758  * ecmp (a, b);
2759  *
2760  *  returns +1 if a > b
2761  *           0 if a == b
2762  *          -1 if a < b
2763  */
2764 int
2765 ecmp (a, b)
2766      unsigned EMUSHORT *a, *b;
2767 {
2768   unsigned EMUSHORT ai[NI], bi[NI];
2769   register unsigned EMUSHORT *p, *q;
2770   register int i;
2771   int msign;
2772
2773   emovi (a, ai);
2774   p = ai;
2775   emovi (b, bi);
2776   q = bi;
2777
2778   if (*p != *q)
2779     {                           /* the signs are different */
2780       /* -0 equals + 0 */
2781       for (i = 1; i < NI - 1; i++)
2782         {
2783           if (ai[i] != 0)
2784             goto nzro;
2785           if (bi[i] != 0)
2786             goto nzro;
2787         }
2788       return (0);
2789     nzro:
2790       if (*p == 0)
2791         return (1);
2792       else
2793         return (-1);
2794     }
2795   /* both are the same sign */
2796   if (*p == 0)
2797     msign = 1;
2798   else
2799     msign = -1;
2800   i = NI - 1;
2801   do
2802     {
2803       if (*p++ != *q++)
2804         {
2805           goto diff;
2806         }
2807     }
2808   while (--i > 0);
2809
2810   return (0);                   /* equality */
2811
2812
2813
2814  diff:
2815
2816   if (*(--p) > *(--q))
2817     return (msign);             /* p is bigger */
2818   else
2819     return (-msign);            /* p is littler */
2820 }
2821
2822
2823
2824
2825 /* Find nearest integer to x = floor (x + 0.5)
2826  *
2827  * unsigned EMUSHORT x[NE], y[NE]
2828  * eround (x, y);
2829  */
2830 void
2831 eround (x, y)
2832      unsigned EMUSHORT *x, *y;
2833 {
2834   eadd (ehalf, x, y);
2835   efloor (y, y);
2836 }
2837
2838
2839
2840
2841 /*
2842 ; convert long integer to e type
2843 ;
2844 ;       long l;
2845 ;       unsigned EMUSHORT x[NE];
2846 ;       ltoe (&l, x);
2847 ; note &l is the memory address of l
2848 */
2849 void
2850 ltoe (lp, y)
2851      long *lp;                  /* lp is the memory address of a long integer */
2852      unsigned EMUSHORT *y;              /* y is the address of a short */
2853 {
2854   unsigned EMUSHORT yi[NI];
2855   unsigned long ll;
2856   int k;
2857
2858   ecleaz (yi);
2859   if (*lp < 0)
2860     {
2861       /* make it positive */
2862       ll = (unsigned long) (-(*lp));
2863       yi[0] = 0xffff;           /* put correct sign in the e type number */
2864     }
2865   else
2866     {
2867       ll = (unsigned long) (*lp);
2868     }
2869   /* move the long integer to yi significand area */
2870   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2871   yi[M + 1] = (unsigned EMUSHORT) ll;
2872
2873   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2874   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2875     ecleaz (yi);                /* it was zero */
2876   else
2877     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
2878   emovo (yi, y);                /* output the answer */
2879 }
2880
2881 /*
2882 ; convert unsigned long integer to e type
2883 ;
2884 ;       unsigned long l;
2885 ;       unsigned EMUSHORT x[NE];
2886 ;       ltox (&l, x);
2887 ; note &l is the memory address of l
2888 */
2889 void
2890 ultoe (lp, y)
2891      unsigned long *lp;         /* lp is the memory address of a long integer */
2892      unsigned EMUSHORT *y;              /* y is the address of a short */
2893 {
2894   unsigned EMUSHORT yi[NI];
2895   unsigned long ll;
2896   int k;
2897
2898   ecleaz (yi);
2899   ll = *lp;
2900
2901   /* move the long integer to ayi significand area */
2902   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2903   yi[M + 1] = (unsigned EMUSHORT) ll;
2904
2905   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2906   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2907     ecleaz (yi);                /* it was zero */
2908   else
2909     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
2910   emovo (yi, y);                /* output the answer */
2911 }
2912
2913
2914 /*
2915 ;       Find long integer and fractional parts
2916
2917 ;       long i;
2918 ;       unsigned EMUSHORT x[NE], frac[NE];
2919 ;       xifrac (x, &i, frac);
2920
2921   The integer output has the sign of the input.  The fraction is
2922 the positive fractional part of abs (x).
2923 */
2924 void
2925 eifrac (x, i, frac)
2926      unsigned EMUSHORT *x;
2927      long *i;
2928      unsigned EMUSHORT *frac;
2929 {
2930   unsigned EMUSHORT xi[NI];
2931   int k;
2932
2933   emovi (x, xi);
2934   k = (int) xi[E] - (EXONE - 1);
2935   if (k <= 0)
2936     {
2937       /* if exponent <= 0, integer = 0 and real output is fraction */
2938       *i = 0L;
2939       emovo (xi, frac);
2940       return;
2941     }
2942   if (k > (HOST_BITS_PER_LONG - 1))
2943     {
2944       /*
2945          ;      long integer overflow: output large integer
2946          ;      and correct fraction
2947          */
2948       if (xi[0])
2949         *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
2950       else
2951         *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
2952       (void) eshift (xi, k);
2953       pedwarn ("Overflow on truncation to integer");
2954       goto lab11;
2955     }
2956
2957   if (k > 16)
2958     {
2959       /*
2960          ; shift more than 16 bits: shift up k-16, output the integer,
2961          ; then complete the shift to get the fraction.
2962          */
2963       k -= 16;
2964       (void) eshift (xi, k);
2965
2966       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
2967       eshup6 (xi);
2968       goto lab10;
2969     }
2970
2971   /* shift not more than 16 bits */
2972   (void) eshift (xi, k);
2973   *i = (long) xi[M] & 0xffff;
2974
2975  lab10:
2976
2977   if (xi[0])
2978     *i = -(*i);
2979  lab11:
2980
2981   xi[0] = 0;
2982   xi[E] = EXONE - 1;
2983   xi[M] = 0;
2984   if ((k = enormlz (xi)) > NBITS)
2985     ecleaz (xi);
2986   else
2987     xi[E] -= (unsigned EMUSHORT) k;
2988
2989   emovo (xi, frac);
2990 }
2991
2992
2993 /*
2994 ;       Find unsigned long integer and fractional parts
2995
2996 ;       unsigned long i;
2997 ;       unsigned EMUSHORT x[NE], frac[NE];
2998 ;       xifrac (x, &i, frac);
2999
3000   A negative e type input yields integer output = 0
3001   but correct fraction.
3002 */
3003 void
3004 euifrac (x, i, frac)
3005      unsigned EMUSHORT *x;
3006      long *i;
3007      unsigned EMUSHORT *frac;
3008 {
3009   unsigned EMUSHORT xi[NI];
3010   int k;
3011
3012   emovi (x, xi);
3013   k = (int) xi[E] - (EXONE - 1);
3014   if (k <= 0)
3015     {
3016       /* if exponent <= 0, integer = 0 and argument is fraction */
3017       *i = 0L;
3018       emovo (xi, frac);
3019       return;
3020     }
3021   if (k > 32)
3022     {
3023       /*
3024          ;      long integer overflow: output large integer
3025          ;      and correct fraction
3026          */
3027       *i = ~(0L);
3028       (void) eshift (xi, k);
3029       pedwarn ("Overflow on truncation to unsigned integer");
3030       goto lab10;
3031     }
3032
3033   if (k > 16)
3034     {
3035       /*
3036          ; shift more than 16 bits: shift up k-16, output the integer,
3037          ; then complete the shift to get the fraction.
3038          */
3039       k -= 16;
3040       (void) eshift (xi, k);
3041
3042       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
3043       eshup6 (xi);
3044       goto lab10;
3045     }
3046
3047   /* shift not more than 16 bits */
3048   (void) eshift (xi, k);
3049   *i = (long) xi[M] & 0xffff;
3050
3051  lab10:
3052
3053   if (xi[0])
3054     *i = 0L;
3055
3056   xi[0] = 0;
3057   xi[E] = EXONE - 1;
3058   xi[M] = 0;
3059   if ((k = enormlz (xi)) > NBITS)
3060     ecleaz (xi);
3061   else
3062     xi[E] -= (unsigned EMUSHORT) k;
3063
3064   emovo (xi, frac);
3065 }
3066
3067
3068
3069 /*
3070 ;       Shift significand
3071 ;
3072 ;       Shifts significand area up or down by the number of bits
3073 ;       given by the variable sc.
3074 */
3075 int
3076 eshift (x, sc)
3077      unsigned EMUSHORT *x;
3078      int sc;
3079 {
3080   unsigned EMUSHORT lost;
3081   unsigned EMUSHORT *p;
3082
3083   if (sc == 0)
3084     return (0);
3085
3086   lost = 0;
3087   p = x + NI - 1;
3088
3089   if (sc < 0)
3090     {
3091       sc = -sc;
3092       while (sc >= 16)
3093         {
3094           lost |= *p;           /* remember lost bits */
3095           eshdn6 (x);
3096           sc -= 16;
3097         }
3098
3099       while (sc >= 8)
3100         {
3101           lost |= *p & 0xff;
3102           eshdn8 (x);
3103           sc -= 8;
3104         }
3105
3106       while (sc > 0)
3107         {
3108           lost |= *p & 1;
3109           eshdn1 (x);
3110           sc -= 1;
3111         }
3112     }
3113   else
3114     {
3115       while (sc >= 16)
3116         {
3117           eshup6 (x);
3118           sc -= 16;
3119         }
3120
3121       while (sc >= 8)
3122         {
3123           eshup8 (x);
3124           sc -= 8;
3125         }
3126
3127       while (sc > 0)
3128         {
3129           eshup1 (x);
3130           sc -= 1;
3131         }
3132     }
3133   if (lost)
3134     lost = 1;
3135   return ((int) lost);
3136 }
3137
3138
3139
3140 /*
3141 ;       normalize
3142 ;
3143 ; Shift normalizes the significand area pointed to by argument
3144 ; shift count (up = positive) is returned.
3145 */
3146 int
3147 enormlz (x)
3148      unsigned EMUSHORT x[];
3149 {
3150   register unsigned EMUSHORT *p;
3151   int sc;
3152
3153   sc = 0;
3154   p = &x[M];
3155   if (*p != 0)
3156     goto normdn;
3157   ++p;
3158   if (*p & 0x8000)
3159     return (0);                 /* already normalized */
3160   while (*p == 0)
3161     {
3162       eshup6 (x);
3163       sc += 16;
3164       /* With guard word, there are NBITS+16 bits available.
3165        * return true if all are zero.
3166        */
3167       if (sc > NBITS)
3168         return (sc);
3169     }
3170   /* see if high byte is zero */
3171   while ((*p & 0xff00) == 0)
3172     {
3173       eshup8 (x);
3174       sc += 8;
3175     }
3176   /* now shift 1 bit at a time */
3177   while ((*p & 0x8000) == 0)
3178     {
3179       eshup1 (x);
3180       sc += 1;
3181       if (sc > NBITS)
3182         {
3183           mtherr ("enormlz", UNDERFLOW);
3184           return (sc);
3185         }
3186     }
3187   return (sc);
3188
3189   /* Normalize by shifting down out of the high guard word
3190      of the significand */
3191  normdn:
3192
3193   if (*p & 0xff00)
3194     {
3195       eshdn8 (x);
3196       sc -= 8;
3197     }
3198   while (*p != 0)
3199     {
3200       eshdn1 (x);
3201       sc -= 1;
3202
3203       if (sc < -NBITS)
3204         {
3205           mtherr ("enormlz", OVERFLOW);
3206           return (sc);
3207         }
3208     }
3209   return (sc);
3210 }
3211
3212
3213
3214
3215 /* Convert e type number to decimal format ASCII string.
3216  * The constants are for 64 bit precision.
3217  */
3218
3219 #define NTEN 12
3220 #define MAXP 4096
3221
3222 static unsigned EMUSHORT etens[NTEN + 1][NE] =
3223 {
3224   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
3225   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
3226   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3227   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3228   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3229   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3230   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3231   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3232   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3233   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3234   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3235   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3236   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
3237 };
3238
3239 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
3240 {
3241   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
3242   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
3243   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3244   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3245   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3246   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3247   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3248   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3249   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3250   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3251   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3252   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3253   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
3254 };
3255
3256 void
3257 e24toasc (x, string, ndigs)
3258      unsigned EMUSHORT x[];
3259      char *string;
3260      int ndigs;
3261 {
3262   unsigned EMUSHORT w[NI];
3263
3264 #ifdef INFINITY
3265 #ifdef IBMPC
3266   if ((x[1] & 0x7f80) == 0x7f80)
3267 #else
3268   if ((x[0] & 0x7f80) == 0x7f80)
3269 #endif
3270     {
3271 #ifdef IBMPC
3272       if (x[1] & 0x8000)
3273 #else
3274       if (x[0] & 0x8000)
3275 #endif
3276         sprintf (string, " -Infinity ");
3277       else
3278         sprintf (string, " Infinity ");
3279       return;
3280     }
3281 #endif
3282   e24toe (x, w);
3283   etoasc (w, string, ndigs);
3284 }
3285
3286
3287 void
3288 e53toasc (x, string, ndigs)
3289      unsigned EMUSHORT x[];
3290      char *string;
3291      int ndigs;
3292 {
3293   unsigned EMUSHORT w[NI];
3294
3295 #ifdef INFINITY
3296 #ifdef IBMPC
3297   if ((x[3] & 0x7ff0) == 0x7ff0)
3298 #else
3299   if ((x[0] & 0x7ff0) == 0x7ff0)
3300 #endif
3301     {
3302 #ifdef IBMPC
3303       if (x[3] & 0x8000)
3304 #else
3305       if (x[0] & 0x8000)
3306 #endif
3307         sprintf (string, " -Infinity ");
3308       else
3309         sprintf (string, " Infinity ");
3310       return;
3311     }
3312 #endif
3313   e53toe (x, w);
3314   etoasc (w, string, ndigs);
3315 }
3316
3317
3318 void
3319 e64toasc (x, string, ndigs)
3320      unsigned EMUSHORT x[];
3321      char *string;
3322      int ndigs;
3323 {
3324   unsigned EMUSHORT w[NI];
3325
3326 #ifdef INFINITY
3327 #ifdef IBMPC
3328   if ((x[4] & 0x7fff) == 0x7fff)
3329 #else
3330   if ((x[0] & 0x7fff) == 0x7fff)
3331 #endif
3332     {
3333 #ifdef IBMPC
3334       if (x[4] & 0x8000)
3335 #else
3336       if (x[0] & 0x8000)
3337 #endif
3338         sprintf (string, " -Infinity ");
3339       else
3340         sprintf (string, " Infinity ");
3341       return;
3342     }
3343 #endif
3344   e64toe (x, w);
3345   etoasc (w, string, ndigs);
3346 }
3347
3348
3349 static char wstring[80];        /* working storage for ASCII output */
3350
3351 void
3352 etoasc (x, string, ndigs)
3353      unsigned EMUSHORT x[];
3354      char *string;
3355      int ndigs;
3356 {
3357   EMUSHORT digit;
3358   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
3359   unsigned EMUSHORT *p, *r, *ten;
3360   unsigned EMUSHORT sign;
3361   int i, j, k, expon, rndsav;
3362   char *s, *ss;
3363   unsigned EMUSHORT m;
3364
3365   ss = string;
3366   s = wstring;
3367   while ((*s++ = *ss++) != '\0')
3368     ;
3369   rndsav = rndprc;
3370   rndprc = NBITS;               /* set to full precision */
3371   emov (x, y);                  /* retain external format */
3372   if (y[NE - 1] & 0x8000)
3373     {
3374       sign = 0xffff;
3375       y[NE - 1] &= 0x7fff;
3376     }
3377   else
3378     {
3379       sign = 0;
3380     }
3381   expon = 0;
3382   ten = &etens[NTEN][0];
3383   emov (eone, t);
3384   /* Test for zero exponent */
3385   if (y[NE - 1] == 0)
3386     {
3387       for (k = 0; k < NE - 1; k++)
3388         {
3389           if (y[k] != 0)
3390             goto tnzro;         /* denormalized number */
3391         }
3392       goto isone;               /* legal all zeros */
3393     }
3394  tnzro:
3395
3396   /* Test for infinity.  Don't bother with illegal infinities.
3397    */
3398   if (y[NE - 1] == 0x7fff)
3399     {
3400       if (sign)
3401         sprintf (wstring, " -Infinity ");
3402       else
3403         sprintf (wstring, " Infinity ");
3404       goto bxit;
3405     }
3406
3407   /* Test for exponent nonzero but significand denormalized.
3408    * This is an error condition.
3409    */
3410   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3411     {
3412       mtherr ("etoasc", DOMAIN);
3413       sprintf (wstring, "NaN");
3414       goto bxit;
3415     }
3416
3417   /* Compare to 1.0 */
3418   i = ecmp (eone, y);
3419   if (i == 0)
3420     goto isone;
3421
3422   if (i < 0)
3423     {                           /* Number is greater than 1 */
3424       /* Convert significand to an integer and strip trailing decimal zeros. */
3425       emov (y, u);
3426       u[NE - 1] = EXONE + NBITS - 1;
3427
3428       p = &etens[NTEN - 4][0];
3429       m = 16;
3430       do
3431         {
3432           ediv (p, u, t);
3433           efloor (t, w);
3434           for (j = 0; j < NE - 1; j++)
3435             {
3436               if (t[j] != w[j])
3437                 goto noint;
3438             }
3439           emov (t, u);
3440           expon += (int) m;
3441         noint:
3442           p += NE;
3443           m >>= 1;
3444         }
3445       while (m != 0);
3446
3447       /* Rescale from integer significand */
3448       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3449       emov (u, y);
3450       /* Find power of 10 */
3451       emov (eone, t);
3452       m = MAXP;
3453       p = &etens[0][0];
3454       while (ecmp (ten, u) <= 0)
3455         {
3456           if (ecmp (p, u) <= 0)
3457             {
3458               ediv (p, u, u);
3459               emul (p, t, t);
3460               expon += (int) m;
3461             }
3462           m >>= 1;
3463           if (m == 0)
3464             break;
3465           p += NE;
3466         }
3467     }
3468   else
3469     {                           /* Number is less than 1.0 */
3470       /* Pad significand with trailing decimal zeros. */
3471       if (y[NE - 1] == 0)
3472         {
3473           while ((y[NE - 2] & 0x8000) == 0)
3474             {
3475               emul (ten, y, y);
3476               expon -= 1;
3477             }
3478         }
3479       else
3480         {
3481           emovi (y, w);
3482           for (i = 0; i < NDEC + 1; i++)
3483             {
3484               if ((w[NI - 1] & 0x7) != 0)
3485                 break;
3486               /* multiply by 10 */
3487               emovz (w, u);
3488               eshdn1 (u);
3489               eshdn1 (u);
3490               eaddm (w, u);
3491               u[1] += 3;
3492               while (u[2] != 0)
3493                 {
3494                   eshdn1 (u);
3495                   u[1] += 1;
3496                 }
3497               if (u[NI - 1] != 0)
3498                 break;
3499               if (eone[NE - 1] <= u[1])
3500                 break;
3501               emovz (u, w);
3502               expon -= 1;
3503             }
3504           emovo (w, y);
3505         }
3506       k = -MAXP;
3507       p = &emtens[0][0];
3508       r = &etens[0][0];
3509       emov (y, w);
3510       emov (eone, t);
3511       while (ecmp (eone, w) > 0)
3512         {
3513           if (ecmp (p, w) >= 0)
3514             {
3515               emul (r, w, w);
3516               emul (r, t, t);
3517               expon += k;
3518             }
3519           k /= 2;
3520           if (k == 0)
3521             break;
3522           p += NE;
3523           r += NE;
3524         }
3525       ediv (t, eone, t);
3526     }
3527  isone:
3528   /* Find the first (leading) digit. */
3529   emovi (t, w);
3530   emovz (w, t);
3531   emovi (y, w);
3532   emovz (w, y);
3533   eiremain (t, y);
3534   digit = equot[NI - 1];
3535   while ((digit == 0) && (ecmp (y, ezero) != 0))
3536     {
3537       eshup1 (y);
3538       emovz (y, u);
3539       eshup1 (u);
3540       eshup1 (u);
3541       eaddm (u, y);
3542       eiremain (t, y);
3543       digit = equot[NI - 1];
3544       expon -= 1;
3545     }
3546   s = wstring;
3547   if (sign)
3548     *s++ = '-';
3549   else
3550     *s++ = ' ';
3551   *s++ = (char) digit + '0';
3552   *s++ = '.';
3553   /* Examine number of digits requested by caller. */
3554   if (ndigs < 0)
3555     ndigs = 0;
3556   if (ndigs > NDEC)
3557     ndigs = NDEC;
3558   /* Generate digits after the decimal point. */
3559   for (k = 0; k <= ndigs; k++)
3560     {
3561       /* multiply current number by 10, without normalizing */
3562       eshup1 (y);
3563       emovz (y, u);
3564       eshup1 (u);
3565       eshup1 (u);
3566       eaddm (u, y);
3567       eiremain (t, y);
3568       *s++ = (char) equot[NI - 1] + '0';
3569     }
3570   digit = equot[NI - 1];
3571   --s;
3572   ss = s;
3573   /* round off the ASCII string */
3574   if (digit > 4)
3575     {
3576       /* Test for critical rounding case in ASCII output. */
3577       if (digit == 5)
3578         {
3579           emovo (y, t);
3580           if (ecmp (t, ezero) != 0)
3581             goto roun;          /* round to nearest */
3582           if ((*(s - 1) & 1) == 0)
3583             goto doexp;         /* round to even */
3584         }
3585       /* Round up and propagate carry-outs */
3586     roun:
3587       --s;
3588       k = *s & 0x7f;
3589       /* Carry out to most significant digit? */
3590       if (k == '.')
3591         {
3592           --s;
3593           k = *s;
3594           k += 1;
3595           *s = (char) k;
3596           /* Most significant digit carries to 10? */
3597           if (k > '9')
3598             {
3599               expon += 1;
3600               *s = '1';
3601             }
3602           goto doexp;
3603         }
3604       /* Round up and carry out from less significant digits */
3605       k += 1;
3606       *s = (char) k;
3607       if (k > '9')
3608         {
3609           *s = '0';
3610           goto roun;
3611         }
3612     }
3613  doexp:
3614   /*
3615      if (expon >= 0)
3616      sprintf (ss, "e+%d", expon);
3617      else
3618      sprintf (ss, "e%d", expon);
3619      */
3620   sprintf (ss, "e%d", expon);
3621  bxit:
3622   rndprc = rndsav;
3623   /* copy out the working string */
3624   s = string;
3625   ss = wstring;
3626   while (*ss == ' ')            /* strip possible leading space */
3627     ++ss;
3628   while ((*s++ = *ss++) != '\0')
3629     ;
3630 }
3631
3632
3633
3634
3635 /*
3636 ;                                                               ASCTOQ
3637 ;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
3638 ;                                       SLM, 3 JAN 78
3639 ;
3640 ;       Convert ASCII string to quadruple precision floating point
3641 ;
3642 ;               Numeric input is free field decimal number
3643 ;               with max of 15 digits with or without
3644 ;               decimal point entered as ASCII from teletype.
3645 ;       Entering E after the number followed by a second
3646 ;       number causes the second number to be interpreted
3647 ;       as a power of 10 to be multiplied by the first number
3648 ;       (i.e., "scientific" notation).
3649 ;
3650 ;       Usage:
3651 ;               asctoq (string, q);
3652 */
3653
3654 /* ASCII to single */
3655 void
3656 asctoe24 (s, y)
3657      char *s;
3658      unsigned EMUSHORT *y;
3659 {
3660   asctoeg (s, y, 24);
3661 }
3662
3663
3664 /* ASCII to double */
3665 void
3666 asctoe53 (s, y)
3667      char *s;
3668      unsigned EMUSHORT *y;
3669 {
3670 #ifdef DEC
3671   asctoeg (s, y, 56);
3672 #else
3673   asctoeg (s, y, 53);
3674 #endif
3675 }
3676
3677
3678 /* ASCII to long double */
3679 void
3680 asctoe64 (s, y)
3681      char *s;
3682      unsigned EMUSHORT *y;
3683 {
3684   asctoeg (s, y, 64);
3685 }
3686
3687 /* ASCII to super double */
3688 void
3689 asctoe (s, y)
3690      char *s;
3691      unsigned EMUSHORT *y;
3692 {
3693   asctoeg (s, y, NBITS);
3694 }
3695
3696 /* Space to make a copy of the input string: */
3697 static char lstr[82];
3698
3699 void
3700 asctoeg (ss, y, oprec)
3701      char *ss;
3702      unsigned EMUSHORT *y;
3703      int oprec;
3704 {
3705   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
3706   int esign, decflg, sgnflg, nexp, exp, prec, lost;
3707   int k, trail, c, rndsav;
3708   EMULONG lexp;
3709   unsigned EMUSHORT nsign, *p;
3710   char *sp, *s;
3711
3712   /* Copy the input string. */
3713   s = ss;
3714   while (*s == ' ')             /* skip leading spaces */
3715     ++s;
3716   sp = lstr;
3717   for (k = 0; k < 79; k++)
3718     {
3719       if ((*sp++ = *s++) == '\0')
3720         break;
3721     }
3722   *sp = '\0';
3723   s = lstr;
3724
3725   rndsav = rndprc;
3726   rndprc = NBITS;               /* Set to full precision */
3727   lost = 0;
3728   nsign = 0;
3729   decflg = 0;
3730   sgnflg = 0;
3731   nexp = 0;
3732   exp = 0;
3733   prec = 0;
3734   ecleaz (yy);
3735   trail = 0;
3736
3737  nxtcom:
3738   k = *s - '0';
3739   if ((k >= 0) && (k <= 9))
3740     {
3741       /* Ignore leading zeros */
3742       if ((prec == 0) && (decflg == 0) && (k == 0))
3743         goto donchr;
3744       /* Identify and strip trailing zeros after the decimal point. */
3745       if ((trail == 0) && (decflg != 0))
3746         {
3747           sp = s;
3748           while ((*sp >= '0') && (*sp <= '9'))
3749             ++sp;
3750           /* Check for syntax error */
3751           c = *sp & 0x7f;
3752           if ((c != 'e') && (c != 'E') && (c != '\0')
3753               && (c != '\n') && (c != '\r') && (c != ' ')
3754               && (c != ','))
3755             goto error;
3756           --sp;
3757           while (*sp == '0')
3758             *sp-- = 'z';
3759           trail = 1;
3760           if (*s == 'z')
3761             goto donchr;
3762         }
3763       /* If enough digits were given to more than fill up the yy register,
3764        * continuing until overflow into the high guard word yy[2]
3765        * guarantees that there will be a roundoff bit at the top
3766        * of the low guard word after normalization.
3767        */
3768       if (yy[2] == 0)
3769         {
3770           if (decflg)
3771             nexp += 1;          /* count digits after decimal point */
3772           eshup1 (yy);          /* multiply current number by 10 */
3773           emovz (yy, xt);
3774           eshup1 (xt);
3775           eshup1 (xt);
3776           eaddm (xt, yy);
3777           ecleaz (xt);
3778           xt[NI - 2] = (unsigned EMUSHORT) k;
3779           eaddm (xt, yy);
3780         }
3781       else
3782         {
3783           lost |= k;
3784         }
3785       prec += 1;
3786       goto donchr;
3787     }
3788
3789   switch (*s)
3790     {
3791     case 'z':
3792       break;
3793     case 'E':
3794     case 'e':
3795       goto expnt;
3796     case '.':                   /* decimal point */
3797       if (decflg)
3798         goto error;
3799       ++decflg;
3800       break;
3801     case '-':
3802       nsign = 0xffff;
3803       if (sgnflg)
3804         goto error;
3805       ++sgnflg;
3806       break;
3807     case '+':
3808       if (sgnflg)
3809         goto error;
3810       ++sgnflg;
3811       break;
3812     case ',':
3813     case ' ':
3814     case '\0':
3815     case '\n':
3816     case '\r':
3817       goto daldone;
3818     case 'i':
3819     case 'I':
3820       ecleaz (yy);
3821       yy[E] = 0x7fff;           /* infinity */
3822       goto aexit;
3823     default:
3824     error:
3825       mtherr ("asctoe", DOMAIN);
3826       eclear (y);
3827       goto aexit;
3828     }
3829  donchr:
3830   ++s;
3831   goto nxtcom;
3832
3833   /* Exponent interpretation */
3834  expnt:
3835
3836   esign = 1;
3837   exp = 0;
3838   ++s;
3839   /* check for + or - */
3840   if (*s == '-')
3841     {
3842       esign = -1;
3843       ++s;
3844     }
3845   if (*s == '+')
3846     ++s;
3847   while ((*s >= '0') && (*s <= '9'))
3848     {
3849       exp *= 10;
3850       exp += *s++ - '0';
3851     }
3852   if (esign < 0)
3853     exp = -exp;
3854
3855  daldone:
3856   nexp = exp - nexp;
3857   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3858   while ((nexp > 0) && (yy[2] == 0))
3859     {
3860       emovz (yy, xt);
3861       eshup1 (xt);
3862       eshup1 (xt);
3863       eaddm (yy, xt);
3864       eshup1 (xt);
3865       if (xt[2] != 0)
3866         break;
3867       nexp -= 1;
3868       emovz (xt, yy);
3869     }
3870   if ((k = enormlz (yy)) > NBITS)
3871     {
3872       ecleaz (yy);
3873       goto aexit;
3874     }
3875   lexp = (EXONE - 1 + NBITS) - k;
3876   emdnorm (yy, lost, 0, lexp, 64);
3877   /* convert to external format */
3878
3879
3880   /* Multiply by 10**nexp.  If precision is 64 bits,
3881    * the maximum relative error incurred in forming 10**n
3882    * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3883    * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3884    * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3885    */
3886   lexp = yy[E];
3887   if (nexp == 0)
3888     {
3889       k = 0;
3890       goto expdon;
3891     }
3892   esign = 1;
3893   if (nexp < 0)
3894     {
3895       nexp = -nexp;
3896       esign = -1;
3897       if (nexp > 4096)
3898         {                       /* Punt.  Can't handle this without 2 divides. */
3899           emovi (etens[0], tt);
3900           lexp -= tt[E];
3901           k = edivm (tt, yy);
3902           lexp += EXONE;
3903           nexp -= 4096;
3904         }
3905     }
3906   p = &etens[NTEN][0];
3907   emov (eone, xt);
3908   exp = 1;
3909   do
3910     {
3911       if (exp & nexp)
3912         emul (p, xt, xt);
3913       p -= NE;
3914       exp = exp + exp;
3915     }
3916   while (exp <= MAXP);
3917
3918   emovi (xt, tt);
3919   if (esign < 0)
3920     {
3921       lexp -= tt[E];
3922       k = edivm (tt, yy);
3923       lexp += EXONE;
3924     }
3925   else
3926     {
3927       lexp += tt[E];
3928       k = emulm (tt, yy);
3929       lexp -= EXONE - 1;
3930     }
3931
3932  expdon:
3933
3934   /* Round and convert directly to the destination type */
3935   if (oprec == 53)
3936     lexp -= EXONE - 0x3ff;
3937   else if (oprec == 24)
3938     lexp -= EXONE - 0177;
3939 #ifdef DEC
3940   else if (oprec == 56)
3941     lexp -= EXONE - 0201;
3942 #endif
3943   rndprc = oprec;
3944   emdnorm (yy, k, 0, lexp, 64);
3945
3946  aexit:
3947
3948   rndprc = rndsav;
3949   yy[0] = nsign;
3950   switch (oprec)
3951     {
3952 #ifdef DEC
3953     case 56:
3954       todec (yy, y);            /* see etodec.c */
3955       break;
3956 #endif
3957     case 53:
3958       toe53 (yy, y);
3959       break;
3960     case 24:
3961       toe24 (yy, y);
3962       break;
3963     case 64:
3964       toe64 (yy, y);
3965       break;
3966     case NBITS:
3967       emovo (yy, y);
3968       break;
3969     }
3970 }
3971
3972
3973
3974 /* y = largest integer not greater than x
3975  * (truncated toward minus infinity)
3976  *
3977  * unsigned EMUSHORT x[NE], y[NE]
3978  *
3979  * efloor (x, y);
3980  */
3981 static unsigned EMUSHORT bmask[] =
3982 {
3983   0xffff,
3984   0xfffe,
3985   0xfffc,
3986   0xfff8,
3987   0xfff0,
3988   0xffe0,
3989   0xffc0,
3990   0xff80,
3991   0xff00,
3992   0xfe00,
3993   0xfc00,
3994   0xf800,
3995   0xf000,
3996   0xe000,
3997   0xc000,
3998   0x8000,
3999   0x0000,
4000 };
4001
4002 void
4003 efloor (x, y)
4004      unsigned EMUSHORT x[], y[];
4005 {
4006   register unsigned EMUSHORT *p;
4007   int e, expon, i;
4008   unsigned EMUSHORT f[NE];
4009
4010   emov (x, f);                  /* leave in external format */
4011   expon = (int) f[NE - 1];
4012   e = (expon & 0x7fff) - (EXONE - 1);
4013   if (e <= 0)
4014     {
4015       eclear (y);
4016       goto isitneg;
4017     }
4018   /* number of bits to clear out */
4019   e = NBITS - e;
4020   emov (f, y);
4021   if (e <= 0)
4022     return;
4023
4024   p = &y[0];
4025   while (e >= 16)
4026     {
4027       *p++ = 0;
4028       e -= 16;
4029     }
4030   /* clear the remaining bits */
4031   *p &= bmask[e];
4032   /* truncate negatives toward minus infinity */
4033  isitneg:
4034
4035   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4036     {
4037       for (i = 0; i < NE - 1; i++)
4038         {
4039           if (f[i] != y[i])
4040             {
4041               esub (eone, y, y);
4042               break;
4043             }
4044         }
4045     }
4046 }
4047
4048
4049 /* unsigned EMUSHORT x[], s[];
4050  * int *exp;
4051  *
4052  * efrexp (x, exp, s);
4053  *
4054  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
4055  * For example, 1.1 = 0.55 * 2**1
4056  * Handles denormalized numbers properly using long integer exp.
4057  */
4058 void
4059 efrexp (x, exp, s)
4060      unsigned EMUSHORT x[];
4061      int *exp;
4062      unsigned EMUSHORT s[];
4063 {
4064   unsigned EMUSHORT xi[NI];
4065   EMULONG li;
4066
4067   emovi (x, xi);
4068   li = (EMULONG) ((EMUSHORT) xi[1]);
4069
4070   if (li == 0)
4071     {
4072       li -= enormlz (xi);
4073     }
4074   xi[1] = 0x3ffe;
4075   emovo (xi, s);
4076   *exp = (int) (li - 0x3ffe);
4077 }
4078
4079
4080
4081 /* unsigned EMUSHORT x[], y[];
4082  * long pwr2;
4083  *
4084  * eldexp (x, pwr2, y);
4085  *
4086  * Returns y = x * 2**pwr2.
4087  */
4088 void
4089 eldexp (x, pwr2, y)
4090      unsigned EMUSHORT x[];
4091      int pwr2;
4092      unsigned EMUSHORT y[];
4093 {
4094   unsigned EMUSHORT xi[NI];
4095   EMULONG li;
4096   int i;
4097
4098   emovi (x, xi);
4099   li = xi[1];
4100   li += pwr2;
4101   i = 0;
4102   emdnorm (xi, i, i, li, 64);
4103   emovo (xi, y);
4104 }
4105
4106
4107 /* c = remainder after dividing b by a
4108  * Least significant integer quotient bits left in equot[].
4109  */
4110 void
4111 eremain (a, b, c)
4112      unsigned EMUSHORT a[], b[], c[];
4113 {
4114   unsigned EMUSHORT den[NI], num[NI];
4115
4116   if (ecmp (a, ezero) == 0)
4117     {
4118       mtherr ("eremain", SING);
4119       eclear (c);
4120       return;
4121     }
4122   emovi (a, den);
4123   emovi (b, num);
4124   eiremain (den, num);
4125   /* Sign of remainder = sign of quotient */
4126   if (a[0] == b[0])
4127     num[0] = 0;
4128   else
4129     num[0] = 0xffff;
4130   emovo (num, c);
4131 }
4132
4133 void
4134 eiremain (den, num)
4135      unsigned EMUSHORT den[], num[];
4136 {
4137   EMULONG ld, ln;
4138   unsigned EMUSHORT j;
4139
4140   ld = den[E];
4141   ld -= enormlz (den);
4142   ln = num[E];
4143   ln -= enormlz (num);
4144   ecleaz (equot);
4145   while (ln >= ld)
4146     {
4147       if (ecmpm (den, num) <= 0)
4148         {
4149           esubm (den, num);
4150           j = 1;
4151         }
4152       else
4153         {
4154           j = 0;
4155         }
4156       eshup1 (equot);
4157       equot[NI - 1] |= j;
4158       eshup1 (num);
4159       ln -= 1;
4160     }
4161   emdnorm (num, 0, 0, ln, 0);
4162 }
4163
4164 /*                                                      mtherr.c
4165  *
4166  *      Library common error handling routine
4167  *
4168  *
4169  *
4170  * SYNOPSIS:
4171  *
4172  * char *fctnam;
4173  * int code;
4174  * void mtherr ();
4175  *
4176  * mtherr (fctnam, code);
4177  *
4178  *
4179  *
4180  * DESCRIPTION:
4181  *
4182  * This routine may be called to report one of the following
4183  * error conditions (in the include file mconf.h).
4184  *
4185  *   Mnemonic        Value          Significance
4186  *
4187  *    DOMAIN            1       argument domain error
4188  *    SING              2       function singularity
4189  *    OVERFLOW          3       overflow range error
4190  *    UNDERFLOW         4       underflow range error
4191  *    TLOSS             5       total loss of precision
4192  *    PLOSS             6       partial loss of precision
4193  *    EDOM             33       Unix domain error code
4194  *    ERANGE           34       Unix range error code
4195  *
4196  * The default version of the file prints the function name,
4197  * passed to it by the pointer fctnam, followed by the
4198  * error condition.  The display is directed to the standard
4199  * output device.  The routine then returns to the calling
4200  * program.  Users may wish to modify the program to abort by
4201  * calling exit under severe error conditions such as domain
4202  * errors.
4203  *
4204  * Since all error conditions pass control to this function,
4205  * the display may be easily changed, eliminated, or directed
4206  * to an error logging device.
4207  *
4208  * SEE ALSO:
4209  *
4210  * mconf.h
4211  *
4212  */
4213 \f
4214 /*
4215 Cephes Math Library Release 2.0:  April, 1987
4216 Copyright 1984, 1987 by Stephen L. Moshier
4217 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4218 */
4219
4220 /* include "mconf.h" */
4221
4222 /* Notice: the order of appearance of the following
4223  * messages is bound to the error codes defined
4224  * in mconf.h.
4225  */
4226 static char *ermsg[7] =
4227 {
4228   "unknown",                    /* error code 0 */
4229   "domain",                     /* error code 1 */
4230   "singularity",                /* et seq.      */
4231   "overflow",
4232   "underflow",
4233   "total loss of precision",
4234   "partial loss of precision"
4235 };
4236
4237 int merror = 0;
4238 extern int merror;
4239
4240 void
4241 mtherr (name, code)
4242      char *name;
4243      int code;
4244 {
4245   char errstr[80];
4246
4247   /* Display string passed by calling program,
4248    * which is supposed to be the name of the
4249    * function in which the error occurred.
4250    */
4251
4252   /* Display error message defined
4253    * by the code argument.
4254    */
4255   if ((code <= 0) || (code >= 6))
4256     code = 0;
4257   sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
4258   pedwarn (errstr);
4259   /* Set global error message word */
4260   merror = code + 1;
4261
4262   /* Return to calling
4263    * program
4264    */
4265 }
4266
4267 /* Here is etodec.c .
4268  *
4269  */
4270
4271 /*
4272 ;       convert DEC double precision to e type
4273 ;       double d;
4274 ;       EMUSHORT e[NE];
4275 ;       dectoe (&d, e);
4276 */
4277 void
4278 dectoe (d, e)
4279      unsigned EMUSHORT *d;
4280      unsigned EMUSHORT *e;
4281 {
4282   unsigned EMUSHORT y[NI];
4283   register unsigned EMUSHORT r, *p;
4284
4285   ecleaz (y);                   /* start with a zero */
4286   p = y;                        /* point to our number */
4287   r = *d;                       /* get DEC exponent word */
4288   if (*d & (unsigned int) 0x8000)
4289     *p = 0xffff;                /* fill in our sign */
4290   ++p;                          /* bump pointer to our exponent word */
4291   r &= 0x7fff;                  /* strip the sign bit */
4292   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
4293     goto done;
4294
4295
4296   r >>= 7;                      /* shift exponent word down 7 bits */
4297   r += EXONE - 0201;            /* subtract DEC exponent offset */
4298   /* add our e type exponent offset */
4299   *p++ = r;                     /* to form our exponent */
4300
4301   r = *d++;                     /* now do the high order mantissa */
4302   r &= 0177;                    /* strip off the DEC exponent and sign bits */
4303   r |= 0200;                    /* the DEC understood high order mantissa bit */
4304   *p++ = r;                     /* put result in our high guard word */
4305
4306   *p++ = *d++;                  /* fill in the rest of our mantissa */
4307   *p++ = *d++;
4308   *p = *d;
4309
4310   eshdn8 (y);                   /* shift our mantissa down 8 bits */
4311  done:
4312   emovo (y, e);
4313 }
4314
4315
4316
4317 /*
4318 ;       convert e type to DEC double precision
4319 ;       double d;
4320 ;       EMUSHORT e[NE];
4321 ;       etodec (e, &d);
4322 */
4323 #if 0
4324 static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4325
4326 void
4327 etodec (x, d)
4328      unsigned EMUSHORT *x, *d;
4329 {
4330   unsigned EMUSHORT xi[NI];
4331   register unsigned EMUSHORT r;
4332   int i, j;
4333
4334   emovi (x, xi);
4335   *d = 0;
4336   if (xi[0] != 0)
4337     *d = 0100000;
4338   r = xi[E];
4339   if (r < (EXONE - 128))
4340     goto zout;
4341   i = xi[M + 4];
4342   if ((i & 0200) != 0)
4343     {
4344       if ((i & 0377) == 0200)
4345         {
4346           if ((i & 0400) != 0)
4347             {
4348               /* check all less significant bits */
4349               for (j = M + 5; j < NI; j++)
4350                 {
4351                   if (xi[j] != 0)
4352                     goto yesrnd;
4353                 }
4354             }
4355           goto nornd;
4356         }
4357     yesrnd:
4358       eaddm (decbit, xi);
4359       r -= enormlz (xi);
4360     }
4361
4362  nornd:
4363
4364   r -= EXONE;
4365   r += 0201;
4366   if (r < 0)
4367     {
4368     zout:
4369       *d++ = 0;
4370       *d++ = 0;
4371       *d++ = 0;
4372       *d++ = 0;
4373       return;
4374     }
4375   if (r >= 0377)
4376     {
4377       *d++ = 077777;
4378       *d++ = -1;
4379       *d++ = -1;
4380       *d++ = -1;
4381       return;
4382     }
4383   r &= 0377;
4384   r <<= 7;
4385   eshup8 (xi);
4386   xi[M] &= 0177;
4387   r |= xi[M];
4388   *d++ |= r;
4389   *d++ = xi[M + 1];
4390   *d++ = xi[M + 2];
4391   *d++ = xi[M + 3];
4392 }
4393
4394 #else
4395
4396 void
4397 etodec (x, d)
4398      unsigned EMUSHORT *x, *d;
4399 {
4400   unsigned EMUSHORT xi[NI];
4401   EMULONG exp;
4402   int rndsav;
4403
4404   emovi (x, xi);
4405   exp = (EMULONG) xi[E] - (EXONE - 0201);       /* adjust exponent for offsets */
4406 /* round off to nearest or even */
4407   rndsav = rndprc;
4408   rndprc = 56;
4409   emdnorm (xi, 0, 0, exp, 64);
4410   rndprc = rndsav;
4411   todec (xi, d);
4412 }
4413
4414 void
4415 todec (x, y)
4416      unsigned EMUSHORT *x, *y;
4417 {
4418   unsigned EMUSHORT i;
4419   unsigned EMUSHORT *p;
4420
4421   p = x;
4422   *y = 0;
4423   if (*p++)
4424     *y = 0100000;
4425   i = *p++;
4426   if (i == 0)
4427     {
4428       *y++ = 0;
4429       *y++ = 0;
4430       *y++ = 0;
4431       *y++ = 0;
4432       return;
4433     }
4434   if (i > 0377)
4435     {
4436       *y++ |= 077777;
4437       *y++ = 0xffff;
4438       *y++ = 0xffff;
4439       *y++ = 0xffff;
4440       return;
4441     }
4442   i &= 0377;
4443   i <<= 7;
4444   eshup8 (x);
4445   x[M] &= 0177;
4446   i |= x[M];
4447   *y++ |= i;
4448   *y++ = x[M + 1];
4449   *y++ = x[M + 2];
4450   *y++ = x[M + 3];
4451 }
4452
4453 #endif /* not 0 */
4454
4455 #endif /* EMU_NON_COMPILE not defined */