gcc/real.c

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