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