![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to real.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [GCC 2.x] / gcc|
Return to real.c CVS log | Up to [GCC 2.x] / gcc |
1.1 ! root 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 ([email protected]). ! 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: ! 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: ! 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: ! 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 */
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.