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