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 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 265 extern 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 /* Constant definitions for math error conditions.  */
 981
 982 #define DOMAIN          1       /* argument domain error */
 983 #define SING            2       /* argument singularity */
 984 #define OVERFLOW        3       /* overflow range error */
 985 #define UNDERFLOW       4       /* underflow range error */
 986 #define TLOSS           5       /* total loss of precision */
 987 #define PLOSS           6       /* partial loss of precision */
 988
 989 #define EDOM            33
 990 #define ERANGE          34
 991
 992 /*  e type constants used by high precision check routines */
 993
 994 /*include "ehead.h"*/
 995 /* 0.0 */
 996 unsigned EMUSHORT ezero[NE] =
 997 {
 998   0, 0000000, 0000000, 0000000, 0000000, 0000000,};
 999 extern unsigned EMUSHORT ezero[];
1000
1001 /* 5.0E-1 */
1002 unsigned EMUSHORT ehalf[NE] =
1003 {
1004   0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1005 extern unsigned EMUSHORT ehalf[];
1006
1007 /* 1.0E0 */
1008 unsigned EMUSHORT eone[NE] =
1009 {
1010   0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1011 extern unsigned EMUSHORT eone[];
1012
1013 /* 2.0E0 */
1014 unsigned EMUSHORT etwo[NE] =
1015 {
1016   0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1017 extern unsigned EMUSHORT etwo[];
1018
1019 /* 3.2E1 */
1020 unsigned EMUSHORT e32[NE] =
1021 {
1022   0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1023 extern unsigned EMUSHORT e32[];
1024
1025 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1026 unsigned EMUSHORT elog2[NE] =
1027 {
1028   0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1029 extern unsigned EMUSHORT elog2[];
1030
1031 /* 1.41421356237309504880168872420969807856967187537695E0 */
1032 unsigned EMUSHORT esqrt2[NE] =
1033 {
1034   0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1035 extern unsigned EMUSHORT esqrt2[];
1036
1037 /* 2/sqrt (PI) =
1038  * 1.12837916709551257389615890312154517168810125865800E0 */
1039 unsigned EMUSHORT eoneopi[NE] =
1040 {
1041   0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
1042 extern unsigned EMUSHORT eoneopi[];
1043
1044 /* 3.14159265358979323846264338327950288419716939937511E0 */
1045 unsigned EMUSHORT epi[NE] =
1046 {
1047   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1048 extern unsigned EMUSHORT epi[];
1049
1050 /* 5.7721566490153286060651209008240243104215933593992E-1 */
1051 unsigned EMUSHORT eeul[NE] =
1052 {
1053   0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
1054 extern unsigned EMUSHORT eeul[];
1055
1056 /*
1057 include "ehead.h"
1058 include "mconf.h"
1059 */
1060
1061
1062
1063 /* Control register for rounding precision.
1064  * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
1065  */
1066 int rndprc = NBITS;
1067 extern int rndprc;
1068
1069 void eaddm (), esubm (), emdnorm (), asctoeg ();
1070 static void toe24 (), toe53 (), toe64 ();
1071 void eremain (), einit (), eiremain ();
1072 int ecmpm (), edivm (), emulm ();
1073 void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
1074 void etodec (), todec (), dectoe ();
1075
1076
1077
1078
1079 void
1080 einit ()
1081 {
1082 }
1083
1084 /*
1085 ; Clear out entire external format number.
1086 ;
1087 ; unsigned EMUSHORT x[];
1088 ; eclear (x);
1089 */
1090
1091 void
1092 eclear (x)
1093      register unsigned EMUSHORT *x;
1094 {
1095   register int i;
1096
1097   for (i = 0; i < NE; i++)
1098     *x++ = 0;
1099 }
1100
1101
1102
1103 /* Move external format number from a to b.
1104  *
1105  * emov (a, b);
1106  */
1107
1108 void
1109 emov (a, b)
1110      register unsigned EMUSHORT *a, *b;
1111 {
1112   register int i;
1113
1114   for (i = 0; i < NE; i++)
1115     *b++ = *a++;
1116 }
1117
1118
1119 /*
1120 ;       Absolute value of external format number
1121 ;
1122 ;       EMUSHORT x[NE];
1123 ;       eabs (x);
1124 */
1125
1126 void
1127 eabs (x)
1128      unsigned EMUSHORT x[];     /* x is the memory address of a short */
1129 {
1130
1131   x[NE - 1] &= 0x7fff;          /* sign is top bit of last word of external format */
1132 }
1133
1134
1135
1136
1137 /*
1138 ;       Negate external format number
1139 ;
1140 ;       unsigned EMUSHORT x[NE];
1141 ;       eneg (x);
1142 */
1143
1144 void
1145 eneg (x)
1146      unsigned EMUSHORT x[];
1147 {
1148
1149   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
1150 }
1151
1152
1153
1154 /* Return 1 if external format number is negative,
1155  * else return zero.
1156  */
1157 int
1158 eisneg (x)
1159      unsigned EMUSHORT x[];
1160 {
1161
1162   if (x[NE - 1] & 0x8000)
1163     return (1);
1164   else
1165     return (0);
1166 }
1167
1168
1169 /* Return 1 if external format number has maximum possible exponent,
1170  * else return zero.
1171  */
1172 int
1173 eisinf (x)
1174      unsigned EMUSHORT x[];
1175 {
1176
1177   if ((x[NE - 1] & 0x7fff) == 0x7fff)
1178     return (1);
1179   else
1180     return (0);
1181 }
1182
1183
1184 /*
1185 ; Fill entire number, including exponent and significand, with
1186 ; largest possible number.  These programs implement a saturation
1187 ; value that is an ordinary, legal number.  A special value
1188 ; "infinity" may also be implemented; this would require tests
1189 ; for that value and implementation of special rules for arithmetic
1190 ; operations involving inifinity.
1191 */
1192
1193 void
1194 einfin (x)
1195      register unsigned EMUSHORT *x;
1196 {
1197   register int i;
1198
1199 #ifdef INFINITY
1200   for (i = 0; i < NE - 1; i++)
1201     *x++ = 0;
1202   *x |= 32767;
1203 #else
1204   for (i = 0; i < NE - 1; i++)
1205     *x++ = 0xffff;
1206   *x |= 32766;
1207   if (rndprc < NBITS)
1208     {
1209       if (rndprc == 64)
1210         {
1211           *(x - 5) = 0;
1212         }
1213       if (rndprc == 53)
1214         {
1215           *(x - 4) = 0xf800;
1216         }
1217       else
1218         {
1219           *(x - 4) = 0;
1220           *(x - 3) = 0;
1221           *(x - 2) = 0xff00;
1222         }
1223     }
1224 #endif
1225 }
1226
1227
1228
1229 /* Move in external format number,
1230  * converting it to internal format.
1231  */
1232 void
1233 emovi (a, b)
1234      unsigned EMUSHORT *a, *b;
1235 {
1236   register unsigned EMUSHORT *p, *q;
1237   int i;
1238
1239   q = b;
1240   p = a + (NE - 1);             /* point to last word of external number */
1241   /* get the sign bit */
1242   if (*p & 0x8000)
1243     *q++ = 0xffff;
1244   else
1245     *q++ = 0;
1246   /* get the exponent */
1247   *q = *p--;
1248   *q++ &= 0x7fff;               /* delete the sign bit */
1249 #ifdef INFINITY
1250   if ((*(q - 1) & 0x7fff) == 0x7fff)
1251     {
1252       for (i = 2; i < NI; i++)
1253         *q++ = 0;
1254       return;
1255     }
1256 #endif
1257   /* clear high guard word */
1258   *q++ = 0;
1259   /* move in the significand */
1260   for (i = 0; i < NE - 1; i++)
1261     *q++ = *p--;
1262   /* clear low guard word */
1263   *q = 0;
1264 }
1265
1266
1267 /* Move internal format number out,
1268  * converting it to external format.
1269  */
1270 void
1271 emovo (a, b)
1272      unsigned EMUSHORT *a, *b;
1273 {
1274   register unsigned EMUSHORT *p, *q;
1275   unsigned EMUSHORT i;
1276
1277   p = a;
1278   q = b + (NE - 1);             /* point to output exponent */
1279   /* combine sign and exponent */
1280   i = *p++;
1281   if (i)
1282     *q-- = *p++ | 0x8000;
1283   else
1284     *q-- = *p++;
1285 #ifdef INFINITY
1286   if (*(p - 1) == 0x7fff)
1287     {
1288       einfin (b);
1289       return;
1290     }
1291 #endif
1292   /* skip over guard word */
1293   ++p;
1294   /* move the significand */
1295   for (i = 0; i < NE - 1; i++)
1296     *q-- = *p++;
1297 }
1298
1299
1300
1301
1302 /* Clear out internal format number.
1303  */
1304
1305 void
1306 ecleaz (xi)
1307      register unsigned EMUSHORT *xi;
1308 {
1309   register int i;
1310
1311   for (i = 0; i < NI; i++)
1312     *xi++ = 0;
1313 }
1314
1315
1316 /* same, but don't touch the sign. */
1317
1318 void
1319 ecleazs (xi)
1320      register unsigned EMUSHORT *xi;
1321 {
1322   register int i;
1323
1324   ++xi;
1325   for (i = 0; i < NI - 1; i++)
1326     *xi++ = 0;
1327 }
1328
1329
1330
1331 /* Move internal format number from a to b.
1332  */
1333 void
1334 emovz (a, b)
1335      register unsigned EMUSHORT *a, *b;
1336 {
1337   register int i;
1338
1339   for (i = 0; i < NI - 1; i++)
1340     *b++ = *a++;
1341   /* clear low guard word */
1342   *b = 0;
1343 }
1344
1345
1346 /*
1347 ;       Compare significands of numbers in internal format.
1348 ;       Guard words are included in the comparison.
1349 ;
1350 ;       unsigned EMUSHORT a[NI], b[NI];
1351 ;       cmpm (a, b);
1352 ;
1353 ;       for the significands:
1354 ;       returns +1 if a > b
1355 ;                0 if a == b
1356 ;               -1 if a < b
1357 */
1358 int
1359 ecmpm (a, b)
1360      register unsigned EMUSHORT *a, *b;
1361 {
1362   int i;
1363
1364   a += M;                       /* skip up to significand area */
1365   b += M;
1366   for (i = M; i < NI; i++)
1367     {
1368       if (*a++ != *b++)
1369         goto difrnt;
1370     }
1371   return (0);
1372
1373  difrnt:
1374   if (*(--a) > *(--b))
1375     return (1);
1376   else
1377     return (-1);
1378 }
1379
1380
1381 /*
1382 ;       Shift significand down by 1 bit
1383 */
1384
1385 void
1386 eshdn1 (x)
1387      register unsigned EMUSHORT *x;
1388 {
1389   register unsigned EMUSHORT bits;
1390   int i;
1391
1392   x += M;                       /* point to significand area */
1393
1394   bits = 0;
1395   for (i = M; i < NI; i++)
1396     {
1397       if (*x & 1)
1398         bits |= 1;
1399       *x >>= 1;
1400       if (bits & 2)
1401         *x |= 0x8000;
1402       bits <<= 1;
1403       ++x;
1404     }
1405 }
1406
1407
1408
1409 /*
1410 ;       Shift significand up by 1 bit
1411 */
1412
1413 void
1414 eshup1 (x)
1415      register unsigned EMUSHORT *x;
1416 {
1417   register unsigned EMUSHORT bits;
1418   int i;
1419
1420   x += NI - 1;
1421   bits = 0;
1422
1423   for (i = M; i < NI; i++)
1424     {
1425       if (*x & 0x8000)
1426         bits |= 1;
1427       *x <<= 1;
1428       if (bits & 2)
1429         *x |= 1;
1430       bits <<= 1;
1431       --x;
1432     }
1433 }
1434
1435
1436
1437 /*
1438 ;       Shift significand down by 8 bits
1439 */
1440
1441 void
1442 eshdn8 (x)
1443      register unsigned EMUSHORT *x;
1444 {
1445   register unsigned EMUSHORT newbyt, oldbyt;
1446   int i;
1447
1448   x += M;
1449   oldbyt = 0;
1450   for (i = M; i < NI; i++)
1451     {
1452       newbyt = *x << 8;
1453       *x >>= 8;
1454       *x |= oldbyt;
1455       oldbyt = newbyt;
1456       ++x;
1457     }
1458 }
1459
1460 /*
1461 ;       Shift significand up by 8 bits
1462 */
1463
1464 void
1465 eshup8 (x)
1466      register unsigned EMUSHORT *x;
1467 {
1468   int i;
1469   register unsigned EMUSHORT newbyt, oldbyt;
1470
1471   x += NI - 1;
1472   oldbyt = 0;
1473
1474   for (i = M; i < NI; i++)
1475     {
1476       newbyt = *x >> 8;
1477       *x <<= 8;
1478       *x |= oldbyt;
1479       oldbyt = newbyt;
1480       --x;
1481     }
1482 }
1483
1484 /*
1485 ;       Shift significand up by 16 bits
1486 */
1487
1488 void
1489 eshup6 (x)
1490      register unsigned EMUSHORT *x;
1491 {
1492   int i;
1493   register unsigned EMUSHORT *p;
1494
1495   p = x + M;
1496   x += M + 1;
1497
1498   for (i = M; i < NI - 1; i++)
1499     *p++ = *x++;
1500
1501   *p = 0;
1502 }
1503
1504 /*
1505 ;       Shift significand down by 16 bits
1506 */
1507
1508 void
1509 eshdn6 (x)
1510      register unsigned EMUSHORT *x;
1511 {
1512   int i;
1513   register unsigned EMUSHORT *p;
1514
1515   x += NI - 1;
1516   p = x + 1;
1517
1518   for (i = M; i < NI - 1; i++)
1519     *(--p) = *(--x);
1520
1521   *(--p) = 0;
1522 }
1523 \f
1524 /*
1525 ;       Add significands
1526 ;       x + y replaces y
1527 */
1528
1529 void
1530 eaddm (x, y)
1531      unsigned EMUSHORT *x, *y;
1532 {
1533   register unsigned EMULONG a;
1534   int i;
1535   unsigned int carry;
1536
1537   x += NI - 1;
1538   y += NI - 1;
1539   carry = 0;
1540   for (i = M; i < NI; i++)
1541     {
1542       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1543       if (a & 0x10000)
1544         carry = 1;
1545       else
1546         carry = 0;
1547       *y = (unsigned EMUSHORT) a;
1548       --x;
1549       --y;
1550     }
1551 }
1552
1553 /*
1554 ;       Subtract significands
1555 ;       y - x replaces y
1556 */
1557
1558 void
1559 esubm (x, y)
1560      unsigned EMUSHORT *x, *y;
1561 {
1562   unsigned EMULONG a;
1563   int i;
1564   unsigned int carry;
1565
1566   x += NI - 1;
1567   y += NI - 1;
1568   carry = 0;
1569   for (i = M; i < NI; i++)
1570     {
1571       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1572       if (a & 0x10000)
1573         carry = 1;
1574       else
1575         carry = 0;
1576       *y = (unsigned EMUSHORT) a;
1577       --x;
1578       --y;
1579     }
1580 }
1581
1582
1583 /* Divide significands */
1584
1585 static unsigned EMUSHORT equot[NI];
1586
1587 int
1588 edivm (den, num)
1589      unsigned EMUSHORT den[], num[];
1590 {
1591   int i;
1592   register unsigned EMUSHORT *p, *q;
1593   unsigned EMUSHORT j;
1594
1595   p = &equot[0];
1596   *p++ = num[0];
1597   *p++ = num[1];
1598
1599   for (i = M; i < NI; i++)
1600     {
1601       *p++ = 0;
1602     }
1603
1604   /* Use faster compare and subtraction if denominator
1605    * has only 15 bits of significance.
1606    */
1607   p = &den[M + 2];
1608   if (*p++ == 0)
1609     {
1610       for (i = M + 3; i < NI; i++)
1611         {
1612           if (*p++ != 0)
1613             goto fulldiv;
1614         }
1615       if ((den[M + 1] & 1) != 0)
1616         goto fulldiv;
1617       eshdn1 (num);
1618       eshdn1 (den);
1619
1620       p = &den[M + 1];
1621       q = &num[M + 1];
1622
1623       for (i = 0; i < NBITS + 2; i++)
1624         {
1625           if (*p <= *q)
1626             {
1627               *q -= *p;
1628               j = 1;
1629             }
1630           else
1631             {
1632               j = 0;
1633             }
1634           eshup1 (equot);
1635           equot[NI - 2] |= j;
1636           eshup1 (num);
1637         }
1638       goto divdon;
1639     }
1640
1641   /* The number of quotient bits to calculate is
1642    * NBITS + 1 scaling guard bit + 1 roundoff bit.
1643    */
1644  fulldiv:
1645
1646   p = &equot[NI - 2];
1647   for (i = 0; i < NBITS + 2; i++)
1648     {
1649       if (ecmpm (den, num) <= 0)
1650         {
1651           esubm (den, num);
1652           j = 1;                /* quotient bit = 1 */
1653         }
1654       else
1655         j = 0;
1656       eshup1 (equot);
1657       *p |= j;
1658       eshup1 (num);
1659     }
1660
1661  divdon:
1662
1663   eshdn1 (equot);
1664   eshdn1 (equot);
1665
1666   /* test for nonzero remainder after roundoff bit */
1667   p = &num[M];
1668   j = 0;
1669   for (i = M; i < NI; i++)
1670     {
1671       j |= *p++;
1672     }
1673   if (j)
1674     j = 1;
1675
1676
1677   for (i = 0; i < NI; i++)
1678     num[i] = equot[i];
1679   return ((int) j);
1680 }
1681
1682
1683 /* Multiply significands */
1684 int
1685 emulm (a, b)
1686      unsigned EMUSHORT a[], b[];
1687 {
1688   unsigned EMUSHORT *p, *q;
1689   int i, j, k;
1690
1691   equot[0] = b[0];
1692   equot[1] = b[1];
1693   for (i = M; i < NI; i++)
1694     equot[i] = 0;
1695
1696   p = &a[NI - 2];
1697   k = NBITS;
1698   while (*p == 0)               /* significand is not supposed to be all zero */
1699     {
1700       eshdn6 (a);
1701       k -= 16;
1702     }
1703   if ((*p & 0xff) == 0)
1704     {
1705       eshdn8 (a);
1706       k -= 8;
1707     }
1708
1709   q = &equot[NI - 1];
1710   j = 0;
1711   for (i = 0; i < k; i++)
1712     {
1713       if (*p & 1)
1714         eaddm (b, equot);
1715       /* remember if there were any nonzero bits shifted out */
1716       if (*q & 1)
1717         j |= 1;
1718       eshdn1 (a);
1719       eshdn1 (equot);
1720     }
1721
1722   for (i = 0; i < NI; i++)
1723     b[i] = equot[i];
1724
1725   /* return flag for lost nonzero bits */
1726   return (j);
1727 }
1728
1729
1730
1731 /*
1732  * Normalize and round off.
1733  *
1734  * The internal format number to be rounded is "s".
1735  * Input "lost" indicates whether or not the number is exact.
1736  * This is the so-called sticky bit.
1737  *
1738  * Input "subflg" indicates whether the number was obtained
1739  * by a subtraction operation.  In that case if lost is nonzero
1740  * then the number is slightly smaller than indicated.
1741  *
1742  * Input "exp" is the biased exponent, which may be negative.
1743  * the exponent field of "s" is ignored but is replaced by
1744  * "exp" as adjusted by normalization and rounding.
1745  *
1746  * Input "rcntrl" is the rounding control.
1747  */
1748
1749 static int rlast = -1;
1750 static int rw = 0;
1751 static unsigned EMUSHORT rmsk = 0;
1752 static unsigned EMUSHORT rmbit = 0;
1753 static unsigned EMUSHORT rebit = 0;
1754 static int re = 0;
1755 static unsigned EMUSHORT rbit[NI];
1756
1757 void
1758 emdnorm (s, lost, subflg, exp, rcntrl)
1759      unsigned EMUSHORT s[];
1760      int lost;
1761      int subflg;
1762      EMULONG exp;
1763      int rcntrl;
1764 {
1765   int i, j;
1766   unsigned EMUSHORT r;
1767
1768   /* Normalize */
1769   j = enormlz (s);
1770
1771   /* a blank significand could mean either zero or infinity. */
1772 #ifndef INFINITY
1773   if (j > NBITS)
1774     {
1775       ecleazs (s);
1776       return;
1777     }
1778 #endif
1779   exp -= j;
1780 #ifndef INFINITY
1781   if (exp >= 32767L)
1782     goto overf;
1783 #else
1784   if ((j > NBITS) && (exp < 32767))
1785     {
1786       ecleazs (s);
1787       return;
1788     }
1789 #endif
1790   if (exp < 0L)
1791     {
1792       if (exp > (EMULONG) (-NBITS - 1))
1793         {
1794           j = (int) exp;
1795           i = eshift (s, j);
1796           if (i)
1797             lost = 1;
1798         }
1799       else
1800         {
1801           ecleazs (s);
1802           return;
1803         }
1804     }
1805   /* Round off, unless told not to by rcntrl. */
1806   if (rcntrl == 0)
1807     goto mdfin;
1808   /* Set up rounding parameters if the control register changed. */
1809   if (rndprc != rlast)
1810     {
1811       ecleaz (rbit);
1812       switch (rndprc)
1813         {
1814         default:
1815         case NBITS:
1816           rw = NI - 1;          /* low guard word */
1817           rmsk = 0xffff;
1818           rmbit = 0x8000;
1819           rbit[rw - 1] = 1;
1820           re = NI - 2;
1821           rebit = 1;
1822           break;
1823         case 64:
1824           rw = 7;
1825           rmsk = 0xffff;
1826           rmbit = 0x8000;
1827           rbit[rw - 1] = 1;
1828           re = rw - 1;
1829           rebit = 1;
1830           break;
1831           /* For DEC arithmetic */
1832         case 56:
1833           rw = 6;
1834           rmsk = 0xff;
1835           rmbit = 0x80;
1836           rbit[rw] = 0x100;
1837           re = rw;
1838           rebit = 0x100;
1839           break;
1840         case 53:
1841           rw = 6;
1842           rmsk = 0x7ff;
1843           rmbit = 0x0400;
1844           rbit[rw] = 0x800;
1845           re = rw;
1846           rebit = 0x800;
1847           break;
1848         case 24:
1849           rw = 4;
1850           rmsk = 0xff;
1851           rmbit = 0x80;
1852           rbit[rw] = 0x100;
1853           re = rw;
1854           rebit = 0x100;
1855           break;
1856         }
1857       rlast = rndprc;
1858     }
1859
1860   if (rndprc >= 64)
1861     {
1862       r = s[rw] & rmsk;
1863       if (rndprc == 64)
1864         {
1865           i = rw + 1;
1866           while (i < NI)
1867             {
1868               if (s[i])
1869                 r |= 1;
1870               s[i] = 0;
1871               ++i;
1872             }
1873         }
1874     }
1875   else
1876     {
1877       if (exp <= 0)
1878         eshdn1 (s);
1879       r = s[rw] & rmsk;
1880       /* These tests assume NI = 8 */
1881       i = rw + 1;
1882       while (i < NI)
1883         {
1884           if (s[i])
1885             r |= 1;
1886           s[i] = 0;
1887           ++i;
1888         }
1889       /*
1890          if (rndprc == 24)
1891          {
1892          if (s[5] || s[6])
1893          r |= 1;
1894          s[5] = 0;
1895          s[6] = 0;
1896          }
1897          */
1898       s[rw] &= ~rmsk;
1899     }
1900   if ((r & rmbit) != 0)
1901     {
1902       if (r == rmbit)
1903         {
1904           if (lost == 0)
1905             {                   /* round to even */
1906               if ((s[re] & rebit) == 0)
1907                 goto mddone;
1908             }
1909           else
1910             {
1911               if (subflg != 0)
1912                 goto mddone;
1913             }
1914         }
1915       eaddm (rbit, s);
1916     }
1917  mddone:
1918   if ((rndprc < 64) && (exp <= 0))
1919     {
1920       eshup1 (s);
1921     }
1922   if (s[2] != 0)
1923     {                           /* overflow on roundoff */
1924       eshdn1 (s);
1925       exp += 1;
1926     }
1927  mdfin:
1928   s[NI - 1] = 0;
1929   if (exp >= 32767L)
1930     {
1931 #ifndef INFINITY
1932     overf:
1933 #endif
1934 #ifdef INFINITY
1935       s[1] = 32767;
1936       for (i = 2; i < NI - 1; i++)
1937         s[i] = 0;
1938       pedwarn ("floating point overflow");
1939 #else
1940       s[1] = 32766;
1941       s[2] = 0;
1942       for (i = M + 1; i < NI - 1; i++)
1943         s[i] = 0xffff;
1944       s[NI - 1] = 0;
1945       if (rndprc < 64)
1946         {
1947           s[rw] &= ~rmsk;
1948           if (rndprc == 24)
1949             {
1950               s[5] = 0;
1951               s[6] = 0;
1952             }
1953         }
1954 #endif
1955       return;
1956     }
1957   if (exp < 0)
1958     s[1] = 0;
1959   else
1960     s[1] = (unsigned EMUSHORT) exp;
1961 }
1962
1963
1964
1965 /*
1966 ;       Subtract external format numbers.
1967 ;
1968 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1969 ;       esub (a, b, c);  c = b - a
1970 */
1971
1972 static int subflg = 0;
1973
1974 void
1975 esub (a, b, c)
1976      unsigned EMUSHORT *a, *b, *c;
1977 {
1978
1979   subflg = 1;
1980   eadd1 (a, b, c);
1981 }
1982
1983
1984 /*
1985 ;       Add.
1986 ;
1987 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
1988 ;       eadd (a, b, c);  c = b + a
1989 */
1990 void
1991 eadd (a, b, c)
1992      unsigned EMUSHORT *a, *b, *c;
1993 {
1994
1995   subflg = 0;
1996   eadd1 (a, b, c);
1997 }
1998
1999 void
2000 eadd1 (a, b, c)
2001      unsigned EMUSHORT *a, *b, *c;
2002 {
2003   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2004   int i, lost, j, k;
2005   EMULONG lt, lta, ltb;
2006
2007 #ifdef INFINITY
2008   if (eisinf (a))
2009     {
2010       emov (a, c);
2011       if (subflg)
2012         eneg (c);
2013       return;
2014     }
2015   if (eisinf (b))
2016     {
2017       emov (b, c);
2018       return;
2019     }
2020 #endif
2021   emovi (a, ai);
2022   emovi (b, bi);
2023   if (subflg)
2024     ai[0] = ~ai[0];
2025
2026   /* compare exponents */
2027   lta = ai[E];
2028   ltb = bi[E];
2029   lt = lta - ltb;
2030   if (lt > 0L)
2031     {                           /* put the larger number in bi */
2032       emovz (bi, ci);
2033       emovz (ai, bi);
2034       emovz (ci, ai);
2035       ltb = bi[E];
2036       lt = -lt;
2037     }
2038   lost = 0;
2039   if (lt != 0L)
2040     {
2041       if (lt < (EMULONG) (-NBITS - 1))
2042         goto done;              /* answer same as larger addend */
2043       k = (int) lt;
2044       lost = eshift (ai, k);    /* shift the smaller number down */
2045     }
2046   else
2047     {
2048       /* exponents were the same, so must compare significands */
2049       i = ecmpm (ai, bi);
2050       if (i == 0)
2051         {                       /* the numbers are identical in magnitude */
2052           /* if different signs, result is zero */
2053           if (ai[0] != bi[0])
2054             {
2055               eclear (c);
2056               return;
2057             }
2058           /* if same sign, result is double */
2059           /* double denomalized tiny number */
2060           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2061             {
2062               eshup1 (bi);
2063               goto done;
2064             }
2065           /* add 1 to exponent unless both are zero! */
2066           for (j = 1; j < NI - 1; j++)
2067             {
2068               if (bi[j] != 0)
2069                 {
2070                   /* This could overflow, but let emovo take care of that. */
2071                   ltb += 1;
2072                   break;
2073                 }
2074             }
2075           bi[E] = (unsigned EMUSHORT) ltb;
2076           goto done;
2077         }
2078       if (i > 0)
2079         {                       /* put the larger number in bi */
2080           emovz (bi, ci);
2081           emovz (ai, bi);
2082           emovz (ci, ai);
2083         }
2084     }
2085   if (ai[0] == bi[0])
2086     {
2087       eaddm (ai, bi);
2088       subflg = 0;
2089     }
2090   else
2091     {
2092       esubm (ai, bi);
2093       subflg = 1;
2094     }
2095   emdnorm (bi, lost, subflg, ltb, 64);
2096
2097  done:
2098   emovo (bi, c);
2099 }
2100
2101
2102
2103 /*
2104 ;       Divide.
2105 ;
2106 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2107 ;       ediv (a, b, c); c = b / a
2108 */
2109 void
2110 ediv (a, b, c)
2111      unsigned EMUSHORT *a, *b, *c;
2112 {
2113   unsigned EMUSHORT ai[NI], bi[NI];
2114   int i;
2115   EMULONG lt, lta, ltb;
2116
2117 #ifdef INFINITY
2118   if (eisinf (b))
2119     {
2120       if (eisneg (a) ^ eisneg (b))
2121         *(c + (NE - 1)) = 0x8000;
2122       else
2123         *(c + (NE - 1)) = 0;
2124       einfin (c);
2125       return;
2126     }
2127   if (eisinf (a))
2128     {
2129       eclear (c);
2130       return;
2131     }
2132 #endif
2133   emovi (a, ai);
2134   emovi (b, bi);
2135   lta = ai[E];
2136   ltb = bi[E];
2137   if (bi[E] == 0)
2138     {                           /* See if numerator is zero. */
2139       for (i = 1; i < NI - 1; i++)
2140         {
2141           if (bi[i] != 0)
2142             {
2143               ltb -= enormlz (bi);
2144               goto dnzro1;
2145             }
2146         }
2147       eclear (c);
2148       return;
2149     }
2150  dnzro1:
2151
2152   if (ai[E] == 0)
2153     {                           /* possible divide by zero */
2154       for (i = 1; i < NI - 1; i++)
2155         {
2156           if (ai[i] != 0)
2157             {
2158               lta -= enormlz (ai);
2159               goto dnzro2;
2160             }
2161         }
2162       if (ai[0] == bi[0])
2163         *(c + (NE - 1)) = 0;
2164       else
2165         *(c + (NE - 1)) = 0x8000;
2166       einfin (c);
2167       mtherr ("ediv", SING);
2168       return;
2169     }
2170  dnzro2:
2171
2172   i = edivm (ai, bi);
2173   /* calculate exponent */
2174   lt = ltb - lta + EXONE;
2175   emdnorm (bi, i, 0, lt, 64);
2176   /* set the sign */
2177   if (ai[0] == bi[0])
2178     bi[0] = 0;
2179   else
2180     bi[0] = 0Xffff;
2181   emovo (bi, c);
2182 }
2183
2184
2185
2186 /*
2187 ;       Multiply.
2188 ;
2189 ;       unsigned EMUSHORT a[NE], b[NE], c[NE];
2190 ;       emul (a, b, c); c = b * a
2191 */
2192 void
2193 emul (a, b, c)
2194      unsigned EMUSHORT *a, *b, *c;
2195 {
2196   unsigned EMUSHORT ai[NI], bi[NI];
2197   int i, j;
2198   EMULONG lt, lta, ltb;
2199
2200 #ifdef INFINITY
2201   if (eisinf (a) || eisinf (b))
2202     {
2203       if (eisneg (a) ^ eisneg (b))
2204         *(c + (NE - 1)) = 0x8000;
2205       else
2206         *(c + (NE - 1)) = 0;
2207       einfin (c);
2208       return;
2209     }
2210 #endif
2211   emovi (a, ai);
2212   emovi (b, bi);
2213   lta = ai[E];
2214   ltb = bi[E];
2215   if (ai[E] == 0)
2216     {
2217       for (i = 1; i < NI - 1; i++)
2218         {
2219           if (ai[i] != 0)
2220             {
2221               lta -= enormlz (ai);
2222               goto mnzer1;
2223             }
2224         }
2225       eclear (c);
2226       return;
2227     }
2228  mnzer1:
2229
2230   if (bi[E] == 0)
2231     {
2232       for (i = 1; i < NI - 1; i++)
2233         {
2234           if (bi[i] != 0)
2235             {
2236               ltb -= enormlz (bi);
2237               goto mnzer2;
2238             }
2239         }
2240       eclear (c);
2241       return;
2242     }
2243  mnzer2:
2244
2245   /* Multiply significands */
2246   j = emulm (ai, bi);
2247   /* calculate exponent */
2248   lt = lta + ltb - (EXONE - 1);
2249   emdnorm (bi, j, 0, lt, 64);
2250   /* calculate sign of product */
2251   if (ai[0] == bi[0])
2252     bi[0] = 0;
2253   else
2254     bi[0] = 0xffff;
2255   emovo (bi, c);
2256 }
2257
2258
2259
2260
2261 /*
2262 ; Convert IEEE double precision to e type
2263 ;       double d;
2264 ;       unsigned EMUSHORT x[N+2];
2265 ;       e53toe (&d, x);
2266 */
2267 void
2268 e53toe (e, y)
2269      unsigned EMUSHORT *e, *y;
2270 {
2271 #ifdef DEC
2272
2273   dectoe (e, y);                /* see etodec.c */
2274
2275 #else
2276
2277   register unsigned EMUSHORT r;
2278   register unsigned EMUSHORT *p;
2279   unsigned EMUSHORT yy[NI];
2280   int denorm, k;
2281
2282   denorm = 0;                   /* flag if denormalized number */
2283   ecleaz (yy);
2284 #ifdef IBMPC
2285   e += 3;
2286 #endif
2287   r = *e;
2288   yy[0] = 0;
2289   if (r & 0x8000)
2290     yy[0] = 0xffff;
2291   yy[M] = (r & 0x0f) | 0x10;
2292   r &= ~0x800f;                 /* strip sign and 4 significand bits */
2293 #ifdef INFINITY
2294   if (r == 0x7ff0)
2295     {
2296       einfin (y);
2297       if (r & 0x8000)
2298         eneg (y);
2299       return;
2300     }
2301 #endif
2302   r >>= 4;
2303   /* If zero exponent, then the significand is denormalized.
2304    * So, take back the understood high significand bit. */
2305   if (r == 0)
2306     {
2307       denorm = 1;
2308       yy[M] &= ~0x10;
2309     }
2310   r += EXONE - 01777;
2311   yy[E] = r;
2312   p = &yy[M + 1];
2313 #ifdef IBMPC
2314   *p++ = *(--e);
2315   *p++ = *(--e);
2316   *p++ = *(--e);
2317 #endif
2318 #ifdef MIEEE
2319   ++e;
2320   *p++ = *e++;
2321   *p++ = *e++;
2322   *p++ = *e++;
2323 #endif
2324   (void) eshift (yy, -5);
2325   if (denorm)
2326     {                           /* if zero exponent, then normalize the significand */
2327       if ((k = enormlz (yy)) > NBITS)
2328         ecleazs (yy);
2329       else
2330         yy[E] -= (unsigned EMUSHORT) (k - 1);
2331     }
2332   emovo (yy, y);
2333 #endif /* not DEC */
2334 }
2335
2336 void
2337 e64toe (e, y)
2338      unsigned EMUSHORT *e, *y;
2339 {
2340   unsigned EMUSHORT yy[NI];
2341   unsigned EMUSHORT *p, *q;
2342   int i;
2343
2344   p = yy;
2345   for (i = 0; i < NE - 5; i++)
2346     *p++ = 0;
2347 #ifdef IBMPC
2348   for (i = 0; i < 5; i++)
2349     *p++ = *e++;
2350 #endif
2351 #ifdef DEC
2352   for (i = 0; i < 5; i++)
2353     *p++ = *e++;
2354 #endif
2355 #ifdef MIEEE
2356   p = &yy[0] + (NE - 1);
2357   *p-- = *e++;
2358   ++e;
2359   for (i = 0; i < 4; i++)
2360     *p-- = *e++;
2361 #endif
2362   p = yy;
2363   q = y;
2364 #ifdef INFINITY
2365   if (*p == 0x7fff)
2366     {
2367       einfin (y);
2368       if (*p & 0x8000)
2369         eneg (y);
2370       return;
2371     }
2372 #endif
2373   for (i = 0; i < NE; i++)
2374     *q++ = *p++;
2375 }
2376
2377
2378 /*
2379 ; Convert IEEE single precision to e type
2380 ;       float d;
2381 ;       unsigned EMUSHORT x[N+2];
2382 ;       dtox (&d, x);
2383 */
2384 void
2385 e24toe (e, y)
2386      unsigned EMUSHORT *e, *y;
2387 {
2388   register unsigned EMUSHORT r;
2389   register unsigned EMUSHORT *p;
2390   unsigned EMUSHORT yy[NI];
2391   int denorm, k;
2392
2393   denorm = 0;                   /* flag if denormalized number */
2394   ecleaz (yy);
2395 #ifdef IBMPC
2396   e += 1;
2397 #endif
2398 #ifdef DEC
2399   e += 1;
2400 #endif
2401   r = *e;
2402   yy[0] = 0;
2403   if (r & 0x8000)
2404     yy[0] = 0xffff;
2405   yy[M] = (r & 0x7f) | 0200;
2406   r &= ~0x807f;                 /* strip sign and 7 significand bits */
2407 #ifdef INFINITY
2408   if (r == 0x7f80)
2409     {
2410       einfin (y);
2411       if (r & 0x8000)
2412         eneg (y);
2413       return;
2414     }
2415 #endif
2416   r >>= 7;
2417   /* If zero exponent, then the significand is denormalized.
2418    * So, take back the understood high significand bit. */
2419   if (r == 0)
2420     {
2421       denorm = 1;
2422       yy[M] &= ~0200;
2423     }
2424   r += EXONE - 0177;
2425   yy[E] = r;
2426   p = &yy[M + 1];
2427 #ifdef IBMPC
2428   *p++ = *(--e);
2429 #endif
2430 #ifdef DEC
2431   *p++ = *(--e);
2432 #endif
2433 #ifdef MIEEE
2434   ++e;
2435   *p++ = *e++;
2436 #endif
2437   (void) eshift (yy, -8);
2438   if (denorm)
2439     {                           /* if zero exponent, then normalize the significand */
2440       if ((k = enormlz (yy)) > NBITS)
2441         ecleazs (yy);
2442       else
2443         yy[E] -= (unsigned EMUSHORT) (k - 1);
2444     }
2445   emovo (yy, y);
2446 }
2447
2448
2449 void
2450 etoe64 (x, e)
2451      unsigned EMUSHORT *x, *e;
2452 {
2453   unsigned EMUSHORT xi[NI];
2454   EMULONG exp;
2455   int rndsav;
2456
2457   emovi (x, xi);
2458   /* adjust exponent for offset */
2459   exp = (EMULONG) xi[E];
2460 #ifdef INFINITY
2461   if (eisinf (x))
2462     goto nonorm;
2463 #endif
2464   /* round off to nearest or even */
2465   rndsav = rndprc;
2466   rndprc = 64;
2467   emdnorm (xi, 0, 0, exp, 64);
2468   rndprc = rndsav;
2469  nonorm:
2470   toe64 (xi, e);
2471 }
2472
2473 /* move out internal format to ieee long double */
2474 static void
2475 toe64 (a, b)
2476      unsigned EMUSHORT *a, *b;
2477 {
2478   register unsigned EMUSHORT *p, *q;
2479   unsigned EMUSHORT i;
2480
2481   p = a;
2482 #ifdef MIEEE
2483   q = b;
2484 #else
2485   q = b + 4;                    /* point to output exponent */
2486 #if LONG_DOUBLE_TYPE_SIZE == 96
2487   /* Clear the last two bytes of 12-byte Intel format */
2488   *(q+1) = 0;
2489 #endif
2490 #endif
2491
2492   /* combine sign and exponent */
2493   i = *p++;
2494 #ifdef MIEEE
2495   if (i)
2496     *q++ = *p++ | 0x8000;
2497   else
2498     *q++ = *p++;
2499   *q++ = 0;
2500 #else
2501   if (i)
2502     *q-- = *p++ | 0x8000;
2503   else
2504     *q-- = *p++;
2505 #endif
2506   /* skip over guard word */
2507   ++p;
2508   /* move the significand */
2509 #ifdef MIEEE
2510   for (i = 0; i < 4; i++)
2511     *q++ = *p++;
2512 #else
2513   for (i = 0; i < 4; i++)
2514     *q-- = *p++;
2515 #endif
2516 }
2517
2518
2519 /*
2520 ; e type to IEEE double precision
2521 ;       double d;
2522 ;       unsigned EMUSHORT x[NE];
2523 ;       etoe53 (x, &d);
2524 */
2525
2526 #ifdef DEC
2527
2528 void
2529 etoe53 (x, e)
2530      unsigned EMUSHORT *x, *e;
2531 {
2532   etodec (x, e);                /* see etodec.c */
2533 }
2534
2535 static void
2536 toe53 (x, y)
2537      unsigned EMUSHORT *x, *y;
2538 {
2539   todec (x, y);
2540 }
2541
2542 #else
2543
2544 void
2545 etoe53 (x, e)
2546      unsigned EMUSHORT *x, *e;
2547 {
2548   unsigned EMUSHORT xi[NI];
2549   EMULONG exp;
2550   int rndsav;
2551
2552   emovi (x, xi);
2553   /* adjust exponent for offsets */
2554   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
2555 #ifdef INFINITY
2556   if (eisinf (x))
2557     goto nonorm;
2558 #endif
2559   /* round off to nearest or even */
2560   rndsav = rndprc;
2561   rndprc = 53;
2562   emdnorm (xi, 0, 0, exp, 64);
2563   rndprc = rndsav;
2564  nonorm:
2565   toe53 (xi, e);
2566 }
2567
2568
2569 static void
2570 toe53 (x, y)
2571      unsigned EMUSHORT *x, *y;
2572 {
2573   unsigned EMUSHORT i;
2574   unsigned EMUSHORT *p;
2575
2576
2577   p = &x[0];
2578 #ifdef IBMPC
2579   y += 3;
2580 #endif
2581   *y = 0;                       /* output high order */
2582   if (*p++)
2583     *y = 0x8000;                /* output sign bit */
2584
2585   i = *p++;
2586   if (i >= (unsigned int) 2047)
2587     {                           /* Saturate at largest number less than infinity. */
2588 #ifdef INFINITY
2589       *y |= 0x7ff0;
2590 #ifdef IBMPC
2591       *(--y) = 0;
2592       *(--y) = 0;
2593       *(--y) = 0;
2594 #endif
2595 #ifdef MIEEE
2596       ++y;
2597       *y++ = 0;
2598       *y++ = 0;
2599       *y++ = 0;
2600 #endif
2601 #else
2602       *y |= (unsigned EMUSHORT) 0x7fef;
2603 #ifdef IBMPC
2604       *(--y) = 0xffff;
2605       *(--y) = 0xffff;
2606       *(--y) = 0xffff;
2607 #endif
2608 #ifdef MIEEE
2609       ++y;
2610       *y++ = 0xffff;
2611       *y++ = 0xffff;
2612       *y++ = 0xffff;
2613 #endif
2614 #endif
2615       return;
2616     }
2617   if (i == 0)
2618     {
2619       (void) eshift (x, 4);
2620     }
2621   else
2622     {
2623       i <<= 4;
2624       (void) eshift (x, 5);
2625     }
2626   i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
2627   *y |= (unsigned EMUSHORT) i;  /* high order output already has sign bit set */
2628 #ifdef IBMPC
2629   *(--y) = *p++;
2630   *(--y) = *p++;
2631   *(--y) = *p;
2632 #endif
2633 #ifdef MIEEE
2634   ++y;
2635   *y++ = *p++;
2636   *y++ = *p++;
2637   *y++ = *p++;
2638 #endif
2639 }
2640
2641 #endif /* not DEC */
2642
2643
2644
2645 /*
2646 ; e type to IEEE single precision
2647 ;       float d;
2648 ;       unsigned EMUSHORT x[N+2];
2649 ;       xtod (x, &d);
2650 */
2651 void
2652 etoe24 (x, e)
2653      unsigned EMUSHORT *x, *e;
2654 {
2655   EMULONG exp;
2656   unsigned EMUSHORT xi[NI];
2657   int rndsav;
2658
2659   emovi (x, xi);
2660   /* adjust exponent for offsets */
2661   exp = (EMULONG) xi[E] - (EXONE - 0177);
2662 #ifdef INFINITY
2663   if (eisinf (x))
2664     goto nonorm;
2665 #endif
2666   /* round off to nearest or even */
2667   rndsav = rndprc;
2668   rndprc = 24;
2669   emdnorm (xi, 0, 0, exp, 64);
2670   rndprc = rndsav;
2671  nonorm:
2672   toe24 (xi, e);
2673 }
2674
2675 static void
2676 toe24 (x, y)
2677      unsigned EMUSHORT *x, *y;
2678 {
2679   unsigned EMUSHORT i;
2680   unsigned EMUSHORT *p;
2681
2682   p = &x[0];
2683 #ifdef IBMPC
2684   y += 1;
2685 #endif
2686 #ifdef DEC
2687   y += 1;
2688 #endif
2689   *y = 0;                       /* output high order */
2690   if (*p++)
2691     *y = 0x8000;                /* output sign bit */
2692
2693   i = *p++;
2694   if (i >= 255)
2695     {                           /* Saturate at largest number less than infinity. */
2696 #ifdef INFINITY
2697       *y |= (unsigned EMUSHORT) 0x7f80;
2698 #ifdef IBMPC
2699       *(--y) = 0;
2700 #endif
2701 #ifdef DEC
2702       *(--y) = 0;
2703 #endif
2704 #ifdef MIEEE
2705       ++y;
2706       *y = 0;
2707 #endif
2708 #else
2709       *y |= (unsigned EMUSHORT) 0x7f7f;
2710 #ifdef IBMPC
2711       *(--y) = 0xffff;
2712 #endif
2713 #ifdef DEC
2714       *(--y) = 0xffff;
2715 #endif
2716 #ifdef MIEEE
2717       ++y;
2718       *y = 0xffff;
2719 #endif
2720 #endif
2721       return;
2722     }
2723   if (i == 0)
2724     {
2725       (void) eshift (x, 7);
2726     }
2727   else
2728     {
2729       i <<= 7;
2730       (void) eshift (x, 8);
2731     }
2732   i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
2733   *y |= i;                      /* high order output already has sign bit set */
2734 #ifdef IBMPC
2735   *(--y) = *p;
2736 #endif
2737 #ifdef DEC
2738   *(--y) = *p;
2739 #endif
2740 #ifdef MIEEE
2741   ++y;
2742   *y = *p;
2743 #endif
2744 }
2745
2746
2747 /* Compare two e type numbers.
2748  *
2749  * unsigned EMUSHORT a[NE], b[NE];
2750  * ecmp (a, b);
2751  *
2752  *  returns +1 if a > b
2753  *           0 if a == b
2754  *          -1 if a < b
2755  */
2756 int
2757 ecmp (a, b)
2758      unsigned EMUSHORT *a, *b;
2759 {
2760   unsigned EMUSHORT ai[NI], bi[NI];
2761   register unsigned EMUSHORT *p, *q;
2762   register int i;
2763   int msign;
2764
2765   emovi (a, ai);
2766   p = ai;
2767   emovi (b, bi);
2768   q = bi;
2769
2770   if (*p != *q)
2771     {                           /* the signs are different */
2772       /* -0 equals + 0 */
2773       for (i = 1; i < NI - 1; i++)
2774         {
2775           if (ai[i] != 0)
2776             goto nzro;
2777           if (bi[i] != 0)
2778             goto nzro;
2779         }
2780       return (0);
2781     nzro:
2782       if (*p == 0)
2783         return (1);
2784       else
2785         return (-1);
2786     }
2787   /* both are the same sign */
2788   if (*p == 0)
2789     msign = 1;
2790   else
2791     msign = -1;
2792   i = NI - 1;
2793   do
2794     {
2795       if (*p++ != *q++)
2796         {
2797           goto diff;
2798         }
2799     }
2800   while (--i > 0);
2801
2802   return (0);                   /* equality */
2803
2804
2805
2806  diff:
2807
2808   if (*(--p) > *(--q))
2809     return (msign);             /* p is bigger */
2810   else
2811     return (-msign);            /* p is littler */
2812 }
2813
2814
2815
2816
2817 /* Find nearest integer to x = floor (x + 0.5)
2818  *
2819  * unsigned EMUSHORT x[NE], y[NE]
2820  * eround (x, y);
2821  */
2822 void
2823 eround (x, y)
2824      unsigned EMUSHORT *x, *y;
2825 {
2826   eadd (ehalf, x, y);
2827   efloor (y, y);
2828 }
2829
2830
2831
2832
2833 /*
2834 ; convert long integer to e type
2835 ;
2836 ;       long l;
2837 ;       unsigned EMUSHORT x[NE];
2838 ;       ltoe (&l, x);
2839 ; note &l is the memory address of l
2840 */
2841 void
2842 ltoe (lp, y)
2843      long *lp;                  /* lp is the memory address of a long integer */
2844      unsigned EMUSHORT *y;              /* y is the address of a short */
2845 {
2846   unsigned EMUSHORT yi[NI];
2847   unsigned long ll;
2848   int k;
2849
2850   ecleaz (yi);
2851   if (*lp < 0)
2852     {
2853       /* make it positive */
2854       ll = (unsigned long) (-(*lp));
2855       yi[0] = 0xffff;           /* put correct sign in the e type number */
2856     }
2857   else
2858     {
2859       ll = (unsigned long) (*lp);
2860     }
2861   /* move the long integer to yi significand area */
2862   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2863   yi[M + 1] = (unsigned EMUSHORT) ll;
2864
2865   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2866   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2867     ecleaz (yi);                /* it was zero */
2868   else
2869     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
2870   emovo (yi, y);                /* output the answer */
2871 }
2872
2873 /*
2874 ; convert unsigned long integer to e type
2875 ;
2876 ;       unsigned long l;
2877 ;       unsigned EMUSHORT x[NE];
2878 ;       ltox (&l, x);
2879 ; note &l is the memory address of l
2880 */
2881 void
2882 ultoe (lp, y)
2883      unsigned long *lp;         /* lp is the memory address of a long integer */
2884      unsigned EMUSHORT *y;              /* y is the address of a short */
2885 {
2886   unsigned EMUSHORT yi[NI];
2887   unsigned long ll;
2888   int k;
2889
2890   ecleaz (yi);
2891   ll = *lp;
2892
2893   /* move the long integer to ayi significand area */
2894   yi[M] = (unsigned EMUSHORT) (ll >> 16);
2895   yi[M + 1] = (unsigned EMUSHORT) ll;
2896
2897   yi[E] = EXONE + 15;           /* exponent if normalize shift count were 0 */
2898   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
2899     ecleaz (yi);                /* it was zero */
2900   else
2901     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
2902   emovo (yi, y);                /* output the answer */
2903 }
2904
2905
2906 /*
2907 ;       Find long integer and fractional parts
2908
2909 ;       long i;
2910 ;       unsigned EMUSHORT x[NE], frac[NE];
2911 ;       xifrac (x, &i, frac);
2912
2913   The integer output has the sign of the input.  The fraction is
2914 the positive fractional part of abs (x).
2915 */
2916 void
2917 eifrac (x, i, frac)
2918      unsigned EMUSHORT *x;
2919      long *i;
2920      unsigned EMUSHORT *frac;
2921 {
2922   unsigned EMUSHORT xi[NI];
2923   int k;
2924
2925   emovi (x, xi);
2926   k = (int) xi[E] - (EXONE - 1);
2927   if (k <= 0)
2928     {
2929       /* if exponent <= 0, integer = 0 and real output is fraction */
2930       *i = 0L;
2931       emovo (xi, frac);
2932       return;
2933     }
2934   if (k > (HOST_BITS_PER_LONG - 1))
2935     {
2936       /*
2937          ;      long integer overflow: output large integer
2938          ;      and correct fraction
2939          */
2940       if (xi[0])
2941         *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
2942       else
2943         *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
2944       (void) eshift (xi, k);
2945       pedwarn ("Overflow on truncation to integer");
2946       goto lab11;
2947     }
2948
2949   if (k > 16)
2950     {
2951       /*
2952          ; shift more than 16 bits: shift up k-16, output the integer,
2953          ; then complete the shift to get the fraction.
2954          */
2955       k -= 16;
2956       (void) eshift (xi, k);
2957
2958       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
2959       eshup6 (xi);
2960       goto lab10;
2961     }
2962
2963   /* shift not more than 16 bits */
2964   (void) eshift (xi, k);
2965   *i = (long) xi[M] & 0xffff;
2966
2967  lab10:
2968
2969   if (xi[0])
2970     *i = -(*i);
2971  lab11:
2972
2973   xi[0] = 0;
2974   xi[E] = EXONE - 1;
2975   xi[M] = 0;
2976   if ((k = enormlz (xi)) > NBITS)
2977     ecleaz (xi);
2978   else
2979     xi[E] -= (unsigned EMUSHORT) k;
2980
2981   emovo (xi, frac);
2982 }
2983
2984
2985 /*
2986 ;       Find unsigned long integer and fractional parts
2987
2988 ;       unsigned long i;
2989 ;       unsigned EMUSHORT x[NE], frac[NE];
2990 ;       xifrac (x, &i, frac);
2991
2992   A negative e type input yields integer output = 0
2993   but correct fraction.
2994 */
2995 void
2996 euifrac (x, i, frac)
2997      unsigned EMUSHORT *x;
2998      long *i;
2999      unsigned EMUSHORT *frac;
3000 {
3001   unsigned EMUSHORT xi[NI];
3002   int k;
3003
3004   emovi (x, xi);
3005   k = (int) xi[E] - (EXONE - 1);
3006   if (k <= 0)
3007     {
3008       /* if exponent <= 0, integer = 0 and argument is fraction */
3009       *i = 0L;
3010       emovo (xi, frac);
3011       return;
3012     }
3013   if (k > 32)
3014     {
3015       /*
3016          ;      long integer overflow: output large integer
3017          ;      and correct fraction
3018          */
3019       *i = ~(0L);
3020       (void) eshift (xi, k);
3021       pedwarn ("Overflow on truncation to unsigned integer");
3022       goto lab10;
3023     }
3024
3025   if (k > 16)
3026     {
3027       /*
3028          ; shift more than 16 bits: shift up k-16, output the integer,
3029          ; then complete the shift to get the fraction.
3030          */
3031       k -= 16;
3032       (void) eshift (xi, k);
3033
3034       *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
3035       eshup6 (xi);
3036       goto lab10;
3037     }
3038
3039   /* shift not more than 16 bits */
3040   (void) eshift (xi, k);
3041   *i = (long) xi[M] & 0xffff;
3042
3043  lab10:
3044
3045   if (xi[0])
3046     *i = 0L;
3047
3048   xi[0] = 0;
3049   xi[E] = EXONE - 1;
3050   xi[M] = 0;
3051   if ((k = enormlz (xi)) > NBITS)
3052     ecleaz (xi);
3053   else
3054     xi[E] -= (unsigned EMUSHORT) k;
3055
3056   emovo (xi, frac);
3057 }
3058
3059
3060
3061 /*
3062 ;       Shift significand
3063 ;
3064 ;       Shifts significand area up or down by the number of bits
3065 ;       given by the variable sc.
3066 */
3067 int
3068 eshift (x, sc)
3069      unsigned EMUSHORT *x;
3070      int sc;
3071 {
3072   unsigned EMUSHORT lost;
3073   unsigned EMUSHORT *p;
3074
3075   if (sc == 0)
3076     return (0);
3077
3078   lost = 0;
3079   p = x + NI - 1;
3080
3081   if (sc < 0)
3082     {
3083       sc = -sc;
3084       while (sc >= 16)
3085         {
3086           lost |= *p;           /* remember lost bits */
3087           eshdn6 (x);
3088           sc -= 16;
3089         }
3090
3091       while (sc >= 8)
3092         {
3093           lost |= *p & 0xff;
3094           eshdn8 (x);
3095           sc -= 8;
3096         }
3097
3098       while (sc > 0)
3099         {
3100           lost |= *p & 1;
3101           eshdn1 (x);
3102           sc -= 1;
3103         }
3104     }
3105   else
3106     {
3107       while (sc >= 16)
3108         {
3109           eshup6 (x);
3110           sc -= 16;
3111         }
3112
3113       while (sc >= 8)
3114         {
3115           eshup8 (x);
3116           sc -= 8;
3117         }
3118
3119       while (sc > 0)
3120         {
3121           eshup1 (x);
3122           sc -= 1;
3123         }
3124     }
3125   if (lost)
3126     lost = 1;
3127   return ((int) lost);
3128 }
3129
3130
3131
3132 /*
3133 ;       normalize
3134 ;
3135 ; Shift normalizes the significand area pointed to by argument
3136 ; shift count (up = positive) is returned.
3137 */
3138 int
3139 enormlz (x)
3140      unsigned EMUSHORT x[];
3141 {
3142   register unsigned EMUSHORT *p;
3143   int sc;
3144
3145   sc = 0;
3146   p = &x[M];
3147   if (*p != 0)
3148     goto normdn;
3149   ++p;
3150   if (*p & 0x8000)
3151     return (0);                 /* already normalized */
3152   while (*p == 0)
3153     {
3154       eshup6 (x);
3155       sc += 16;
3156       /* With guard word, there are NBITS+16 bits available.
3157        * return true if all are zero.
3158        */
3159       if (sc > NBITS)
3160         return (sc);
3161     }
3162   /* see if high byte is zero */
3163   while ((*p & 0xff00) == 0)
3164     {
3165       eshup8 (x);
3166       sc += 8;
3167     }
3168   /* now shift 1 bit at a time */
3169   while ((*p & 0x8000) == 0)
3170     {
3171       eshup1 (x);
3172       sc += 1;
3173       if (sc > NBITS)
3174         {
3175           mtherr ("enormlz", UNDERFLOW);
3176           return (sc);
3177         }
3178     }
3179   return (sc);
3180
3181   /* Normalize by shifting down out of the high guard word
3182      of the significand */
3183  normdn:
3184
3185   if (*p & 0xff00)
3186     {
3187       eshdn8 (x);
3188       sc -= 8;
3189     }
3190   while (*p != 0)
3191     {
3192       eshdn1 (x);
3193       sc -= 1;
3194
3195       if (sc < -NBITS)
3196         {
3197           mtherr ("enormlz", OVERFLOW);
3198           return (sc);
3199         }
3200     }
3201   return (sc);
3202 }
3203
3204
3205
3206
3207 /* Convert e type number to decimal format ASCII string.
3208  * The constants are for 64 bit precision.
3209  */
3210
3211 #define NTEN 12
3212 #define MAXP 4096
3213
3214 static unsigned EMUSHORT etens[NTEN + 1][NE] =
3215 {
3216   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
3217   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
3218   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
3219   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
3220   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
3221   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
3222   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
3223   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
3224   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
3225   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
3226   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
3227   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3228   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
3229 };
3230
3231 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
3232 {
3233   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
3234   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
3235   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3236   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3237   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3238   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3239   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3240   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3241   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3242   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3243   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3244   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3245   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
3246 };
3247
3248 void
3249 e24toasc (x, string, ndigs)
3250      unsigned EMUSHORT x[];
3251      char *string;
3252      int ndigs;
3253 {
3254   unsigned EMUSHORT w[NI];
3255
3256 #ifdef INFINITY
3257 #ifdef IBMPC
3258   if ((x[1] & 0x7f80) == 0x7f80)
3259 #else
3260   if ((x[0] & 0x7f80) == 0x7f80)
3261 #endif
3262     {
3263 #ifdef IBMPC
3264       if (x[1] & 0x8000)
3265 #else
3266       if (x[0] & 0x8000)
3267 #endif
3268         sprintf (string, " -Infinity ");
3269       else
3270         sprintf (string, " Infinity ");
3271       return;
3272     }
3273 #endif
3274   e24toe (x, w);
3275   etoasc (w, string, ndigs);
3276 }
3277
3278
3279 void
3280 e53toasc (x, string, ndigs)
3281      unsigned EMUSHORT x[];
3282      char *string;
3283      int ndigs;
3284 {
3285   unsigned EMUSHORT w[NI];
3286
3287 #ifdef INFINITY
3288 #ifdef IBMPC
3289   if ((x[3] & 0x7ff0) == 0x7ff0)
3290 #else
3291   if ((x[0] & 0x7ff0) == 0x7ff0)
3292 #endif
3293     {
3294 #ifdef IBMPC
3295       if (x[3] & 0x8000)
3296 #else
3297       if (x[0] & 0x8000)
3298 #endif
3299         sprintf (string, " -Infinity ");
3300       else
3301         sprintf (string, " Infinity ");
3302       return;
3303     }
3304 #endif
3305   e53toe (x, w);
3306   etoasc (w, string, ndigs);
3307 }
3308
3309
3310 void
3311 e64toasc (x, string, ndigs)
3312      unsigned EMUSHORT x[];
3313      char *string;
3314      int ndigs;
3315 {
3316   unsigned EMUSHORT w[NI];
3317
3318 #ifdef INFINITY
3319 #ifdef IBMPC
3320   if ((x[4] & 0x7fff) == 0x7fff)
3321 #else
3322   if ((x[0] & 0x7fff) == 0x7fff)
3323 #endif
3324     {
3325 #ifdef IBMPC
3326       if (x[4] & 0x8000)
3327 #else
3328       if (x[0] & 0x8000)
3329 #endif
3330         sprintf (string, " -Infinity ");
3331       else
3332         sprintf (string, " Infinity ");
3333       return;
3334     }
3335 #endif
3336   e64toe (x, w);
3337   etoasc (w, string, ndigs);
3338 }
3339
3340
3341 static char wstring[80];        /* working storage for ASCII output */
3342
3343 void
3344 etoasc (x, string, ndigs)
3345      unsigned EMUSHORT x[];
3346      char *string;
3347      int ndigs;
3348 {
3349   EMUSHORT digit;
3350   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
3351   unsigned EMUSHORT *p, *r, *ten;
3352   unsigned EMUSHORT sign;
3353   int i, j, k, expon, rndsav;
3354   char *s, *ss;
3355   unsigned EMUSHORT m;
3356
3357   ss = string;
3358   s = wstring;
3359   while ((*s++ = *ss++) != '\0')
3360     ;
3361   rndsav = rndprc;
3362   rndprc = NBITS;               /* set to full precision */
3363   emov (x, y);                  /* retain external format */
3364   if (y[NE - 1] & 0x8000)
3365     {
3366       sign = 0xffff;
3367       y[NE - 1] &= 0x7fff;
3368     }
3369   else
3370     {
3371       sign = 0;
3372     }
3373   expon = 0;
3374   ten = &etens[NTEN][0];
3375   emov (eone, t);
3376   /* Test for zero exponent */
3377   if (y[NE - 1] == 0)
3378     {
3379       for (k = 0; k < NE - 1; k++)
3380         {
3381           if (y[k] != 0)
3382             goto tnzro;         /* denormalized number */
3383         }
3384       goto isone;               /* legal all zeros */
3385     }
3386  tnzro:
3387
3388   /* Test for infinity.  Don't bother with illegal infinities.
3389    */
3390   if (y[NE - 1] == 0x7fff)
3391     {
3392       if (sign)
3393         sprintf (wstring, " -Infinity ");
3394       else
3395         sprintf (wstring, " Infinity ");
3396       goto bxit;
3397     }
3398
3399   /* Test for exponent nonzero but significand denormalized.
3400    * This is an error condition.
3401    */
3402   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
3403     {
3404       mtherr ("etoasc", DOMAIN);
3405       sprintf (wstring, "NaN");
3406       goto bxit;
3407     }
3408
3409   /* Compare to 1.0 */
3410   i = ecmp (eone, y);
3411   if (i == 0)
3412     goto isone;
3413
3414   if (i < 0)
3415     {                           /* Number is greater than 1 */
3416       /* Convert significand to an integer and strip trailing decimal zeros. */
3417       emov (y, u);
3418       u[NE - 1] = EXONE + NBITS - 1;
3419
3420       p = &etens[NTEN - 4][0];
3421       m = 16;
3422       do
3423         {
3424           ediv (p, u, t);
3425           efloor (t, w);
3426           for (j = 0; j < NE - 1; j++)
3427             {
3428               if (t[j] != w[j])
3429                 goto noint;
3430             }
3431           emov (t, u);
3432           expon += (int) m;
3433         noint:
3434           p += NE;
3435           m >>= 1;
3436         }
3437       while (m != 0);
3438
3439       /* Rescale from integer significand */
3440       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
3441       emov (u, y);
3442       /* Find power of 10 */
3443       emov (eone, t);
3444       m = MAXP;
3445       p = &etens[0][0];
3446       while (ecmp (ten, u) <= 0)
3447         {
3448           if (ecmp (p, u) <= 0)
3449             {
3450               ediv (p, u, u);
3451               emul (p, t, t);
3452               expon += (int) m;
3453             }
3454           m >>= 1;
3455           if (m == 0)
3456             break;
3457           p += NE;
3458         }
3459     }
3460   else
3461     {                           /* Number is less than 1.0 */
3462       /* Pad significand with trailing decimal zeros. */
3463       if (y[NE - 1] == 0)
3464         {
3465           while ((y[NE - 2] & 0x8000) == 0)
3466             {
3467               emul (ten, y, y);
3468               expon -= 1;
3469             }
3470         }
3471       else
3472         {
3473           emovi (y, w);
3474           for (i = 0; i < NDEC + 1; i++)
3475             {
3476               if ((w[NI - 1] & 0x7) != 0)
3477                 break;
3478               /* multiply by 10 */
3479               emovz (w, u);
3480               eshdn1 (u);
3481               eshdn1 (u);
3482               eaddm (w, u);
3483               u[1] += 3;
3484               while (u[2] != 0)
3485                 {
3486                   eshdn1 (u);
3487                   u[1] += 1;
3488                 }
3489               if (u[NI - 1] != 0)
3490                 break;
3491               if (eone[NE - 1] <= u[1])
3492                 break;
3493               emovz (u, w);
3494               expon -= 1;
3495             }
3496           emovo (w, y);
3497         }
3498       k = -MAXP;
3499       p = &emtens[0][0];
3500       r = &etens[0][0];
3501       emov (y, w);
3502       emov (eone, t);
3503       while (ecmp (eone, w) > 0)
3504         {
3505           if (ecmp (p, w) >= 0)
3506             {
3507               emul (r, w, w);
3508               emul (r, t, t);
3509               expon += k;
3510             }
3511           k /= 2;
3512           if (k == 0)
3513             break;
3514           p += NE;
3515           r += NE;
3516         }
3517       ediv (t, eone, t);
3518     }
3519  isone:
3520   /* Find the first (leading) digit. */
3521   emovi (t, w);
3522   emovz (w, t);
3523   emovi (y, w);
3524   emovz (w, y);
3525   eiremain (t, y);
3526   digit = equot[NI - 1];
3527   while ((digit == 0) && (ecmp (y, ezero) != 0))
3528     {
3529       eshup1 (y);
3530       emovz (y, u);
3531       eshup1 (u);
3532       eshup1 (u);
3533       eaddm (u, y);
3534       eiremain (t, y);
3535       digit = equot[NI - 1];
3536       expon -= 1;
3537     }
3538   s = wstring;
3539   if (sign)
3540     *s++ = '-';
3541   else
3542     *s++ = ' ';
3543   *s++ = (char) digit + '0';
3544   *s++ = '.';
3545   /* Examine number of digits requested by caller. */
3546   if (ndigs < 0)
3547     ndigs = 0;
3548   if (ndigs > NDEC)
3549     ndigs = NDEC;
3550   /* Generate digits after the decimal point. */
3551   for (k = 0; k <= ndigs; k++)
3552     {
3553       /* multiply current number by 10, without normalizing */
3554       eshup1 (y);
3555       emovz (y, u);
3556       eshup1 (u);
3557       eshup1 (u);
3558       eaddm (u, y);
3559       eiremain (t, y);
3560       *s++ = (char) equot[NI - 1] + '0';
3561     }
3562   digit = equot[NI - 1];
3563   --s;
3564   ss = s;
3565   /* round off the ASCII string */
3566   if (digit > 4)
3567     {
3568       /* Test for critical rounding case in ASCII output. */
3569       if (digit == 5)
3570         {
3571           emovo (y, t);
3572           if (ecmp (t, ezero) != 0)
3573             goto roun;          /* round to nearest */
3574           if ((*(s - 1) & 1) == 0)
3575             goto doexp;         /* round to even */
3576         }
3577       /* Round up and propagate carry-outs */
3578     roun:
3579       --s;
3580       k = *s & 0x7f;
3581       /* Carry out to most significant digit? */
3582       if (k == '.')
3583         {
3584           --s;
3585           k = *s;
3586           k += 1;
3587           *s = (char) k;
3588           /* Most significant digit carries to 10? */
3589           if (k > '9')
3590             {
3591               expon += 1;
3592               *s = '1';
3593             }
3594           goto doexp;
3595         }
3596       /* Round up and carry out from less significant digits */
3597       k += 1;
3598       *s = (char) k;
3599       if (k > '9')
3600         {
3601           *s = '0';
3602           goto roun;
3603         }
3604     }
3605  doexp:
3606   /*
3607      if (expon >= 0)
3608      sprintf (ss, "e+%d", expon);
3609      else
3610      sprintf (ss, "e%d", expon);
3611      */
3612   sprintf (ss, "e%d", expon);
3613  bxit:
3614   rndprc = rndsav;
3615   /* copy out the working string */
3616   s = string;
3617   ss = wstring;
3618   while (*ss == ' ')            /* strip possible leading space */
3619     ++ss;
3620   while ((*s++ = *ss++) != '\0')
3621     ;
3622 }
3623
3624
3625
3626
3627 /*
3628 ;                                                               ASCTOQ
3629 ;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
3630 ;                                       SLM, 3 JAN 78
3631 ;
3632 ;       Convert ASCII string to quadruple precision floating point
3633 ;
3634 ;               Numeric input is free field decimal number
3635 ;               with max of 15 digits with or without
3636 ;               decimal point entered as ASCII from teletype.
3637 ;       Entering E after the number followed by a second
3638 ;       number causes the second number to be interpreted
3639 ;       as a power of 10 to be multiplied by the first number
3640 ;       (i.e., "scientific" notation).
3641 ;
3642 ;       Usage:
3643 ;               asctoq (string, q);
3644 */
3645
3646 /* ASCII to single */
3647 void
3648 asctoe24 (s, y)
3649      char *s;
3650      unsigned EMUSHORT *y;
3651 {
3652   asctoeg (s, y, 24);
3653 }
3654
3655
3656 /* ASCII to double */
3657 void
3658 asctoe53 (s, y)
3659      char *s;
3660      unsigned EMUSHORT *y;
3661 {
3662 #ifdef DEC
3663   asctoeg (s, y, 56);
3664 #else
3665   asctoeg (s, y, 53);
3666 #endif
3667 }
3668
3669
3670 /* ASCII to long double */
3671 void
3672 asctoe64 (s, y)
3673      char *s;
3674      unsigned EMUSHORT *y;
3675 {
3676   asctoeg (s, y, 64);
3677 }
3678
3679 /* ASCII to super double */
3680 void
3681 asctoe (s, y)
3682      char *s;
3683      unsigned EMUSHORT *y;
3684 {
3685   asctoeg (s, y, NBITS);
3686 }
3687
3688 /* Space to make a copy of the input string: */
3689 static char lstr[82];
3690
3691 void
3692 asctoeg (ss, y, oprec)
3693      char *ss;
3694      unsigned EMUSHORT *y;
3695      int oprec;
3696 {
3697   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
3698   int esign, decflg, sgnflg, nexp, exp, prec, lost;
3699   int k, trail, c, rndsav;
3700   EMULONG lexp;
3701   unsigned EMUSHORT nsign, *p;
3702   char *sp, *s;
3703
3704   /* Copy the input string. */
3705   s = ss;
3706   while (*s == ' ')             /* skip leading spaces */
3707     ++s;
3708   sp = lstr;
3709   for (k = 0; k < 79; k++)
3710     {
3711       if ((*sp++ = *s++) == '\0')
3712         break;
3713     }
3714   *sp = '\0';
3715   s = lstr;
3716
3717   rndsav = rndprc;
3718   rndprc = NBITS;               /* Set to full precision */
3719   lost = 0;
3720   nsign = 0;
3721   decflg = 0;
3722   sgnflg = 0;
3723   nexp = 0;
3724   exp = 0;
3725   prec = 0;
3726   ecleaz (yy);
3727   trail = 0;
3728
3729  nxtcom:
3730   k = *s - '0';
3731   if ((k >= 0) && (k <= 9))
3732     {
3733       /* Ignore leading zeros */
3734       if ((prec == 0) && (decflg == 0) && (k == 0))
3735         goto donchr;
3736       /* Identify and strip trailing zeros after the decimal point. */
3737       if ((trail == 0) && (decflg != 0))
3738         {
3739           sp = s;
3740           while ((*sp >= '0') && (*sp <= '9'))
3741             ++sp;
3742           /* Check for syntax error */
3743           c = *sp & 0x7f;
3744           if ((c != 'e') && (c != 'E') && (c != '\0')
3745               && (c != '\n') && (c != '\r') && (c != ' ')
3746               && (c != ','))
3747             goto error;
3748           --sp;
3749           while (*sp == '0')
3750             *sp-- = 'z';
3751           trail = 1;
3752           if (*s == 'z')
3753             goto donchr;
3754         }
3755       /* If enough digits were given to more than fill up the yy register,
3756        * continuing until overflow into the high guard word yy[2]
3757        * guarantees that there will be a roundoff bit at the top
3758        * of the low guard word after normalization.
3759        */
3760       if (yy[2] == 0)
3761         {
3762           if (decflg)
3763             nexp += 1;          /* count digits after decimal point */
3764           eshup1 (yy);          /* multiply current number by 10 */
3765           emovz (yy, xt);
3766           eshup1 (xt);
3767           eshup1 (xt);
3768           eaddm (xt, yy);
3769           ecleaz (xt);
3770           xt[NI - 2] = (unsigned EMUSHORT) k;
3771           eaddm (xt, yy);
3772         }
3773       else
3774         {
3775           lost |= k;
3776         }
3777       prec += 1;
3778       goto donchr;
3779     }
3780
3781   switch (*s)
3782     {
3783     case 'z':
3784       break;
3785     case 'E':
3786     case 'e':
3787       goto expnt;
3788     case '.':                   /* decimal point */
3789       if (decflg)
3790         goto error;
3791       ++decflg;
3792       break;
3793     case '-':
3794       nsign = 0xffff;
3795       if (sgnflg)
3796         goto error;
3797       ++sgnflg;
3798       break;
3799     case '+':
3800       if (sgnflg)
3801         goto error;
3802       ++sgnflg;
3803       break;
3804     case ',':
3805     case ' ':
3806     case '\0':
3807     case '\n':
3808     case '\r':
3809       goto daldone;
3810     case 'i':
3811     case 'I':
3812       ecleaz (yy);
3813       yy[E] = 0x7fff;           /* infinity */
3814       goto aexit;
3815     default:
3816     error:
3817       mtherr ("asctoe", DOMAIN);
3818       eclear (y);
3819       goto aexit;
3820     }
3821  donchr:
3822   ++s;
3823   goto nxtcom;
3824
3825   /* Exponent interpretation */
3826  expnt:
3827
3828   esign = 1;
3829   exp = 0;
3830   ++s;
3831   /* check for + or - */
3832   if (*s == '-')
3833     {
3834       esign = -1;
3835       ++s;
3836     }
3837   if (*s == '+')
3838     ++s;
3839   while ((*s >= '0') && (*s <= '9'))
3840     {
3841       exp *= 10;
3842       exp += *s++ - '0';
3843     }
3844   if (esign < 0)
3845     exp = -exp;
3846
3847  daldone:
3848   nexp = exp - nexp;
3849   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3850   while ((nexp > 0) && (yy[2] == 0))
3851     {
3852       emovz (yy, xt);
3853       eshup1 (xt);
3854       eshup1 (xt);
3855       eaddm (yy, xt);
3856       eshup1 (xt);
3857       if (xt[2] != 0)
3858         break;
3859       nexp -= 1;
3860       emovz (xt, yy);
3861     }
3862   if ((k = enormlz (yy)) > NBITS)
3863     {
3864       ecleaz (yy);
3865       goto aexit;
3866     }
3867   lexp = (EXONE - 1 + NBITS) - k;
3868   emdnorm (yy, lost, 0, lexp, 64);
3869   /* convert to external format */
3870
3871
3872   /* Multiply by 10**nexp.  If precision is 64 bits,
3873    * the maximum relative error incurred in forming 10**n
3874    * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3875    * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3876    * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3877    */
3878   lexp = yy[E];
3879   if (nexp == 0)
3880     {
3881       k = 0;
3882       goto expdon;
3883     }
3884   esign = 1;
3885   if (nexp < 0)
3886     {
3887       nexp = -nexp;
3888       esign = -1;
3889       if (nexp > 4096)
3890         {                       /* Punt.  Can't handle this without 2 divides. */
3891           emovi (etens[0], tt);
3892           lexp -= tt[E];
3893           k = edivm (tt, yy);
3894           lexp += EXONE;
3895           nexp -= 4096;
3896         }
3897     }
3898   p = &etens[NTEN][0];
3899   emov (eone, xt);
3900   exp = 1;
3901   do
3902     {
3903       if (exp & nexp)
3904         emul (p, xt, xt);
3905       p -= NE;
3906       exp = exp + exp;
3907     }
3908   while (exp <= MAXP);
3909
3910   emovi (xt, tt);
3911   if (esign < 0)
3912     {
3913       lexp -= tt[E];
3914       k = edivm (tt, yy);
3915       lexp += EXONE;
3916     }
3917   else
3918     {
3919       lexp += tt[E];
3920       k = emulm (tt, yy);
3921       lexp -= EXONE - 1;
3922     }
3923
3924  expdon:
3925
3926   /* Round and convert directly to the destination type */
3927   if (oprec == 53)
3928     lexp -= EXONE - 0x3ff;
3929   else if (oprec == 24)
3930     lexp -= EXONE - 0177;
3931 #ifdef DEC
3932   else if (oprec == 56)
3933     lexp -= EXONE - 0201;
3934 #endif
3935   rndprc = oprec;
3936   emdnorm (yy, k, 0, lexp, 64);
3937
3938  aexit:
3939
3940   rndprc = rndsav;
3941   yy[0] = nsign;
3942   switch (oprec)
3943     {
3944 #ifdef DEC
3945     case 56:
3946       todec (yy, y);            /* see etodec.c */
3947       break;
3948 #endif
3949     case 53:
3950       toe53 (yy, y);
3951       break;
3952     case 24:
3953       toe24 (yy, y);
3954       break;
3955     case 64:
3956       toe64 (yy, y);
3957       break;
3958     case NBITS:
3959       emovo (yy, y);
3960       break;
3961     }
3962 }
3963
3964
3965
3966 /* y = largest integer not greater than x
3967  * (truncated toward minus infinity)
3968  *
3969  * unsigned EMUSHORT x[NE], y[NE]
3970  *
3971  * efloor (x, y);
3972  */
3973 static unsigned EMUSHORT bmask[] =
3974 {
3975   0xffff,
3976   0xfffe,
3977   0xfffc,
3978   0xfff8,
3979   0xfff0,
3980   0xffe0,
3981   0xffc0,
3982   0xff80,
3983   0xff00,
3984   0xfe00,
3985   0xfc00,
3986   0xf800,
3987   0xf000,
3988   0xe000,
3989   0xc000,
3990   0x8000,
3991   0x0000,
3992 };
3993
3994 void
3995 efloor (x, y)
3996      unsigned EMUSHORT x[], y[];
3997 {
3998   register unsigned EMUSHORT *p;
3999   int e, expon, i;
4000   unsigned EMUSHORT f[NE];
4001
4002   emov (x, f);                  /* leave in external format */
4003   expon = (int) f[NE - 1];
4004   e = (expon & 0x7fff) - (EXONE - 1);
4005   if (e <= 0)
4006     {
4007       eclear (y);
4008       goto isitneg;
4009     }
4010   /* number of bits to clear out */
4011   e = NBITS - e;
4012   emov (f, y);
4013   if (e <= 0)
4014     return;
4015
4016   p = &y[0];
4017   while (e >= 16)
4018     {
4019       *p++ = 0;
4020       e -= 16;
4021     }
4022   /* clear the remaining bits */
4023   *p &= bmask[e];
4024   /* truncate negatives toward minus infinity */
4025  isitneg:
4026
4027   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
4028     {
4029       for (i = 0; i < NE - 1; i++)
4030         {
4031           if (f[i] != y[i])
4032             {
4033               esub (eone, y, y);
4034               break;
4035             }
4036         }
4037     }
4038 }
4039
4040
4041 /* unsigned EMUSHORT x[], s[];
4042  * int *exp;
4043  *
4044  * efrexp (x, exp, s);
4045  *
4046  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
4047  * For example, 1.1 = 0.55 * 2**1
4048  * Handles denormalized numbers properly using long integer exp.
4049  */
4050 void
4051 efrexp (x, exp, s)
4052      unsigned EMUSHORT x[];
4053      int *exp;
4054      unsigned EMUSHORT s[];
4055 {
4056   unsigned EMUSHORT xi[NI];
4057   EMULONG li;
4058
4059   emovi (x, xi);
4060   li = (EMULONG) ((EMUSHORT) xi[1]);
4061
4062   if (li == 0)
4063     {
4064       li -= enormlz (xi);
4065     }
4066   xi[1] = 0x3ffe;
4067   emovo (xi, s);
4068   *exp = (int) (li - 0x3ffe);
4069 }
4070
4071
4072
4073 /* unsigned EMUSHORT x[], y[];
4074  * long pwr2;
4075  *
4076  * eldexp (x, pwr2, y);
4077  *
4078  * Returns y = x * 2**pwr2.
4079  */
4080 void
4081 eldexp (x, pwr2, y)
4082      unsigned EMUSHORT x[];
4083      int pwr2;
4084      unsigned EMUSHORT y[];
4085 {
4086   unsigned EMUSHORT xi[NI];
4087   EMULONG li;
4088   int i;
4089
4090   emovi (x, xi);
4091   li = xi[1];
4092   li += pwr2;
4093   i = 0;
4094   emdnorm (xi, i, i, li, 64);
4095   emovo (xi, y);
4096 }
4097
4098
4099 /* c = remainder after dividing b by a
4100  * Least significant integer quotient bits left in equot[].
4101  */
4102 void
4103 eremain (a, b, c)
4104      unsigned EMUSHORT a[], b[], c[];
4105 {
4106   unsigned EMUSHORT den[NI], num[NI];
4107
4108   if (ecmp (a, ezero) == 0)
4109     {
4110       mtherr ("eremain", SING);
4111       eclear (c);
4112       return;
4113     }
4114   emovi (a, den);
4115   emovi (b, num);
4116   eiremain (den, num);
4117   /* Sign of remainder = sign of quotient */
4118   if (a[0] == b[0])
4119     num[0] = 0;
4120   else
4121     num[0] = 0xffff;
4122   emovo (num, c);
4123 }
4124
4125 void
4126 eiremain (den, num)
4127      unsigned EMUSHORT den[], num[];
4128 {
4129   EMULONG ld, ln;
4130   unsigned EMUSHORT j;
4131
4132   ld = den[E];
4133   ld -= enormlz (den);
4134   ln = num[E];
4135   ln -= enormlz (num);
4136   ecleaz (equot);
4137   while (ln >= ld)
4138     {
4139       if (ecmpm (den, num) <= 0)
4140         {
4141           esubm (den, num);
4142           j = 1;
4143         }
4144       else
4145         {
4146           j = 0;
4147         }
4148       eshup1 (equot);
4149       equot[NI - 1] |= j;
4150       eshup1 (num);
4151       ln -= 1;
4152     }
4153   emdnorm (num, 0, 0, ln, 0);
4154 }
4155
4156 /*                                                      mtherr.c
4157  *
4158  *      Library common error handling routine
4159  *
4160  *
4161  *
4162  * SYNOPSIS:
4163  *
4164  * char *fctnam;
4165  * int code;
4166  * void mtherr ();
4167  *
4168  * mtherr (fctnam, code);
4169  *
4170  *
4171  *
4172  * DESCRIPTION:
4173  *
4174  * This routine may be called to report one of the following
4175  * error conditions (in the include file mconf.h).
4176  *
4177  *   Mnemonic        Value          Significance
4178  *
4179  *    DOMAIN            1       argument domain error
4180  *    SING              2       function singularity
4181  *    OVERFLOW          3       overflow range error
4182  *    UNDERFLOW         4       underflow range error
4183  *    TLOSS             5       total loss of precision
4184  *    PLOSS             6       partial loss of precision
4185  *    EDOM             33       Unix domain error code
4186  *    ERANGE           34       Unix range error code
4187  *
4188  * The default version of the file prints the function name,
4189  * passed to it by the pointer fctnam, followed by the
4190  * error condition.  The display is directed to the standard
4191  * output device.  The routine then returns to the calling
4192  * program.  Users may wish to modify the program to abort by
4193  * calling exit under severe error conditions such as domain
4194  * errors.
4195  *
4196  * Since all error conditions pass control to this function,
4197  * the display may be easily changed, eliminated, or directed
4198  * to an error logging device.
4199  *
4200  * SEE ALSO:
4201  *
4202  * mconf.h
4203  *
4204  */
4205 \f
4206 /*
4207 Cephes Math Library Release 2.0:  April, 1987
4208 Copyright 1984, 1987 by Stephen L. Moshier
4209 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
4210 */
4211
4212 /* include "mconf.h" */
4213
4214 /* Notice: the order of appearance of the following
4215  * messages is bound to the error codes defined
4216  * in mconf.h.
4217  */
4218 static char *ermsg[7] =
4219 {
4220   "unknown",                    /* error code 0 */
4221   "domain",                     /* error code 1 */
4222   "singularity",                /* et seq.      */
4223   "overflow",
4224   "underflow",
4225   "total loss of precision",
4226   "partial loss of precision"
4227 };
4228
4229 int merror = 0;
4230 extern int merror;
4231
4232 void
4233 mtherr (name, code)
4234      char *name;
4235      int code;
4236 {
4237   char errstr[80];
4238
4239   /* Display string passed by calling program,
4240    * which is supposed to be the name of the
4241    * function in which the error occurred.
4242    */
4243
4244   /* Display error message defined
4245    * by the code argument.
4246    */
4247   if ((code <= 0) || (code >= 6))
4248     code = 0;
4249   sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
4250   pedwarn (errstr);
4251   /* Set global error message word */
4252   merror = code + 1;
4253
4254   /* Return to calling
4255    * program
4256    */
4257 }
4258
4259 /* Here is etodec.c .
4260  *
4261  */
4262
4263 /*
4264 ;       convert DEC double precision to e type
4265 ;       double d;
4266 ;       EMUSHORT e[NE];
4267 ;       dectoe (&d, e);
4268 */
4269 void
4270 dectoe (d, e)
4271      unsigned EMUSHORT *d;
4272      unsigned EMUSHORT *e;
4273 {
4274   unsigned EMUSHORT y[NI];
4275   register unsigned EMUSHORT r, *p;
4276
4277   ecleaz (y);                   /* start with a zero */
4278   p = y;                        /* point to our number */
4279   r = *d;                       /* get DEC exponent word */
4280   if (*d & (unsigned int) 0x8000)
4281     *p = 0xffff;                /* fill in our sign */
4282   ++p;                          /* bump pointer to our exponent word */
4283   r &= 0x7fff;                  /* strip the sign bit */
4284   if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
4285     goto done;
4286
4287
4288   r >>= 7;                      /* shift exponent word down 7 bits */
4289   r += EXONE - 0201;            /* subtract DEC exponent offset */
4290   /* add our e type exponent offset */
4291   *p++ = r;                     /* to form our exponent */
4292
4293   r = *d++;                     /* now do the high order mantissa */
4294   r &= 0177;                    /* strip off the DEC exponent and sign bits */
4295   r |= 0200;                    /* the DEC understood high order mantissa bit */
4296   *p++ = r;                     /* put result in our high guard word */
4297
4298   *p++ = *d++;                  /* fill in the rest of our mantissa */
4299   *p++ = *d++;
4300   *p = *d;
4301
4302   eshdn8 (y);                   /* shift our mantissa down 8 bits */
4303  done:
4304   emovo (y, e);
4305 }
4306
4307
4308
4309 /*
4310 ;       convert e type to DEC double precision
4311 ;       double d;
4312 ;       EMUSHORT e[NE];
4313 ;       etodec (e, &d);
4314 */
4315 #if 0
4316 static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
4317
4318 void
4319 etodec (x, d)
4320      unsigned EMUSHORT *x, *d;
4321 {
4322   unsigned EMUSHORT xi[NI];
4323   register unsigned EMUSHORT r;
4324   int i, j;
4325
4326   emovi (x, xi);
4327   *d = 0;
4328   if (xi[0] != 0)
4329     *d = 0100000;
4330   r = xi[E];
4331   if (r < (EXONE - 128))
4332     goto zout;
4333   i = xi[M + 4];
4334   if ((i & 0200) != 0)
4335     {
4336       if ((i & 0377) == 0200)
4337         {
4338           if ((i & 0400) != 0)
4339             {
4340               /* check all less significant bits */
4341               for (j = M + 5; j < NI; j++)
4342                 {
4343                   if (xi[j] != 0)
4344                     goto yesrnd;
4345                 }
4346             }
4347           goto nornd;
4348         }
4349     yesrnd:
4350       eaddm (decbit, xi);
4351       r -= enormlz (xi);
4352     }
4353
4354  nornd:
4355
4356   r -= EXONE;
4357   r += 0201;
4358   if (r < 0)
4359     {
4360     zout:
4361       *d++ = 0;
4362       *d++ = 0;
4363       *d++ = 0;
4364       *d++ = 0;
4365       return;
4366     }
4367   if (r >= 0377)
4368     {
4369       *d++ = 077777;
4370       *d++ = -1;
4371       *d++ = -1;
4372       *d++ = -1;
4373       return;
4374     }
4375   r &= 0377;
4376   r <<= 7;
4377   eshup8 (xi);
4378   xi[M] &= 0177;
4379   r |= xi[M];
4380   *d++ |= r;
4381   *d++ = xi[M + 1];
4382   *d++ = xi[M + 2];
4383   *d++ = xi[M + 3];
4384 }
4385
4386 #else
4387
4388 void
4389 etodec (x, d)
4390      unsigned EMUSHORT *x, *d;
4391 {
4392   unsigned EMUSHORT xi[NI];
4393   EMULONG exp;
4394   int rndsav;
4395
4396   emovi (x, xi);
4397   exp = (EMULONG) xi[E] - (EXONE - 0201);       /* adjust exponent for offsets */
4398 /* round off to nearest or even */
4399   rndsav = rndprc;
4400   rndprc = 56;
4401   emdnorm (xi, 0, 0, exp, 64);
4402   rndprc = rndsav;
4403   todec (xi, d);
4404 }
4405
4406 void
4407 todec (x, y)
4408      unsigned EMUSHORT *x, *y;
4409 {
4410   unsigned EMUSHORT i;
4411   unsigned EMUSHORT *p;
4412
4413   p = x;
4414   *y = 0;
4415   if (*p++)
4416     *y = 0100000;
4417   i = *p++;
4418   if (i == 0)
4419     {
4420       *y++ = 0;
4421       *y++ = 0;
4422       *y++ = 0;
4423       *y++ = 0;
4424       return;
4425     }
4426   if (i > 0377)
4427     {
4428       *y++ |= 077777;
4429       *y++ = 0xffff;
4430       *y++ = 0xffff;
4431       *y++ = 0xffff;
4432       return;
4433     }
4434   i &= 0377;
4435   i <<= 7;
4436   eshup8 (x);
4437   x[M] &= 0177;
4438   i |= x[M];
4439   *y++ |= i;
4440   *y++ = x[M + 1];
4441   *y++ = x[M + 2];
4442   *y++ = x[M + 3];
4443 }
4444
4445 #endif /* not 0 */
4446
4447 #endif /* EMU_NON_COMPILE not defined */