![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fp-bit.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [GCC 2.x] / gcc / config|
Return to fp-bit.c CVS log | Up to [GCC 2.x] / gcc / config |
1.1 ! root 1: /* This is a software floating point library which can be used instead of ! 2: the floating point routines in libgcc1.c for targets without hardware ! 3: floating point. */ ! 4: ! 5: /* Copyright (C) 1994 Free Software Foundation, Inc. ! 6: ! 7: This file is free software; you can redistribute it and/or modify it ! 8: under the terms of the GNU General Public License as published by the ! 9: Free Software Foundation; either version 2, or (at your option) any ! 10: later version. ! 11: ! 12: In addition to the permissions in the GNU General Public License, the ! 13: Free Software Foundation gives you unlimited permission to link the ! 14: compiled version of this file with other programs, and to distribute ! 15: those programs without any restriction coming from the use of this ! 16: file. (The General Public License restrictions do apply in other ! 17: respects; for example, they cover modification of the file, and ! 18: distribution when not linked into another program.) ! 19: ! 20: This file is distributed in the hope that it will be useful, but ! 21: WITHOUT ANY WARRANTY; without even the implied warranty of ! 22: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 23: General Public License for more details. ! 24: ! 25: You should have received a copy of the GNU General Public License ! 26: along with this program; see the file COPYING. If not, write to ! 27: the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ ! 28: ! 29: /* As a special exception, if you link this library with other files, ! 30: some of which are compiled with GCC, to produce an executable, ! 31: this library does not by itself cause the resulting executable ! 32: to be covered by the GNU General Public License. ! 33: This exception does not however invalidate any other reasons why ! 34: the executable file might be covered by the GNU General Public License. */ ! 35: ! 36: /* This implements IEEE 754 format arithmetic, but does not provide a ! 37: mechanism for setting the rounding mode, or for generating or handling ! 38: exceptions. ! 39: ! 40: The original code by Steve Chamberlain, hacked by Mark Eichin and Jim ! 41: Wilson, all of Cygnus Support. */ ! 42: ! 43: /* The intended way to use this file is to make two copies, add `#define FLOAT' ! 44: to one copy, then compile both copies and add them to libgcc.a. */ ! 45: ! 46: /* The following macros can be defined to change the behaviour of this file: ! 47: FLOAT: Implement a `float', aka SFmode, fp library. If this is not ! 48: defined, then this file implements a `double', aka DFmode, fp library. ! 49: FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e. ! 50: don't include float->double conversion which requires the double library. ! 51: This is useful only for machines which can't support doubles, e.g. some ! 52: 8-bit processors. ! 53: CMPtype: Specify the type that floating point compares should return. ! 54: This defaults to SItype, aka int. ! 55: US_SOFTWARE_GOFAST: This makes all entry points use the same names as the ! 56: US Software goFast library. If this is not defined, the entry points use ! 57: the same names as libgcc1.c. ! 58: _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding ! 59: two integers to the FLO_union_type. ! 60: NO_NANS: Disable nan and infinity handling ! 61: SMALL_MACHINE: Useful when operations on QIs and HIs are faster ! 62: than on an SI */ ! 63: ! 64: typedef SFtype __attribute__ ((mode (SF))); ! 65: typedef DFtype __attribute__ ((mode (DF))); ! 66: ! 67: typedef int HItype __attribute__ ((mode (HI))); ! 68: typedef int SItype __attribute__ ((mode (SI))); ! 69: typedef int DItype __attribute__ ((mode (DI))); ! 70: ! 71: /* The type of the result of a fp compare */ ! 72: #ifndef CMPtype ! 73: #define CMPtype SItype ! 74: #endif ! 75: ! 76: typedef unsigned int UHItype __attribute__ ((mode (HI))); ! 77: typedef unsigned int USItype __attribute__ ((mode (SI))); ! 78: typedef unsigned int UDItype __attribute__ ((mode (DI))); ! 79: ! 80: #define MAX_SI_INT ((SItype) ((unsigned) (~0)>>1)) ! 81: #define MAX_USI_INT ((USItype) ~0) ! 82: ! 83: ! 84: #ifdef FLOAT_ONLY ! 85: #define NO_DI_MODE ! 86: #endif ! 87: ! 88: #ifdef FLOAT ! 89: # define NGARDS 7L ! 90: # define GARDROUND 0x3f ! 91: # define GARDMASK 0x7f ! 92: # define GARDMSB 0x40 ! 93: # define EXPBITS 8 ! 94: # define EXPBIAS 127 ! 95: # define FRACBITS 23 ! 96: # define EXPMAX (0xff) ! 97: # define QUIET_NAN 0x100000L ! 98: # define FRAC_NBITS 32 ! 99: # define FRACHIGH 0x80000000L ! 100: # define FRACHIGH2 0xc0000000L ! 101: typedef USItype fractype; ! 102: typedef UHItype halffractype; ! 103: typedef SFtype FLO_type; ! 104: typedef SItype intfrac; ! 105: ! 106: #else ! 107: # define PREFIXFPDP dp ! 108: # define PREFIXSFDF df ! 109: # define NGARDS 8L ! 110: # define GARDROUND 0x7f ! 111: # define GARDMASK 0xff ! 112: # define GARDMSB 0x80 ! 113: # define EXPBITS 11 ! 114: # define EXPBIAS 1023 ! 115: # define FRACBITS 52 ! 116: # define EXPMAX (0x7ff) ! 117: # define QUIET_NAN 0x8000000000000LL ! 118: # define FRAC_NBITS 64 ! 119: # define FRACHIGH 0x8000000000000000LL ! 120: # define FRACHIGH2 0xc000000000000000LL ! 121: typedef UDItype fractype; ! 122: typedef USItype halffractype; ! 123: typedef DFtype FLO_type; ! 124: typedef DItype intfrac; ! 125: #endif ! 126: ! 127: #ifdef US_SOFTWARE_GOFAST ! 128: # ifdef FLOAT ! 129: # define add fpadd ! 130: # define sub fpsub ! 131: # define multiply fpmul ! 132: # define divide fpdiv ! 133: # define compare fpcmp ! 134: # define si_to_float sitofp ! 135: # define float_to_si fptosi ! 136: # define float_to_usi fptoui ! 137: # define negate __negsf2 ! 138: # define sf_to_df fptodp ! 139: # define dptofp dptofp ! 140: #else ! 141: # define add dpadd ! 142: # define sub dpsub ! 143: # define multiply dpmul ! 144: # define divide dpdiv ! 145: # define compare dpcmp ! 146: # define si_to_float litodp ! 147: # define float_to_si dptoli ! 148: # define float_to_usi dptoul ! 149: # define negate __negdf2 ! 150: # define df_to_sf dptofp ! 151: #endif ! 152: #else ! 153: # ifdef FLOAT ! 154: # define add __addsf3 ! 155: # define sub __subsf3 ! 156: # define multiply __mulsf3 ! 157: # define divide __divsf3 ! 158: # define compare __cmpsf2 ! 159: # define _eq_f2 __eqsf2 ! 160: # define _ne_f2 __nesf2 ! 161: # define _gt_f2 __gtsf2 ! 162: # define _ge_f2 __gesf2 ! 163: # define _lt_f2 __ltsf2 ! 164: # define _le_f2 __lesf2 ! 165: # define si_to_float __floatsisf ! 166: # define float_to_si __fixsfsi ! 167: # define float_to_usi __fixunssfsi ! 168: # define negate __negsf2 ! 169: # define sf_to_df __extendsfdf2 ! 170: #else ! 171: # define add __adddf3 ! 172: # define sub __subdf3 ! 173: # define multiply __muldf3 ! 174: # define divide __divdf3 ! 175: # define compare __cmpdf2 ! 176: # define _eq_f2 __eqdf2 ! 177: # define _ne_f2 __nedf2 ! 178: # define _gt_f2 __gtdf2 ! 179: # define _ge_f2 __gedf2 ! 180: # define _lt_f2 __ltdf2 ! 181: # define _le_f2 __ledf2 ! 182: # define si_to_float __floatsidf ! 183: # define float_to_si __fixdfsi ! 184: # define float_to_usi __fixunsdfsi ! 185: # define negate __negdf2 ! 186: # define df_to_sf __truncdfsf2 ! 187: # endif ! 188: #endif ! 189: ! 190: ! 191: #define INLINE __inline__ ! 192: ! 193: /* Preserve the sticky-bit when shifting fractions to the right. */ ! 194: #define LSHIFT(a) { a = (a & 1) | (a >> 1); } ! 195: ! 196: /* numeric parameters */ ! 197: /* F_D_BITOFF is the number of bits offset between the MSB of the mantissa ! 198: of a float and of a double. Assumes there are only two float types. ! 199: (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS)) ! 200: */ ! 201: #define F_D_BITOFF (52+8-(23+7)) ! 202: ! 203: ! 204: #define NORMAL_EXPMIN (-(EXPBIAS)+1) ! 205: #define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS)) ! 206: #define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS)) ! 207: ! 208: /* common types */ ! 209: ! 210: typedef enum ! 211: { ! 212: CLASS_SNAN, ! 213: CLASS_QNAN, ! 214: CLASS_ZERO, ! 215: CLASS_NUMBER, ! 216: CLASS_INFINITY ! 217: } fp_class_type; ! 218: ! 219: typedef struct ! 220: { ! 221: #ifdef SMALL_MACHINE ! 222: char class; ! 223: unsigned char sign; ! 224: short normal_exp; ! 225: #else ! 226: fp_class_type class; ! 227: unsigned int sign; ! 228: int normal_exp; ! 229: #endif ! 230: ! 231: union ! 232: { ! 233: fractype ll; ! 234: halffractype l[2]; ! 235: } fraction; ! 236: } fp_number_type; ! 237: ! 238: typedef union ! 239: { ! 240: FLO_type value; ! 241: #ifdef _DEBUG_BITFLOAT ! 242: int l[2]; ! 243: #endif ! 244: struct ! 245: { ! 246: #ifndef FLOAT_BIT_ORDER_MISMATCH ! 247: unsigned int sign:1 __attribute__ ((packed)); ! 248: unsigned int exp:EXPBITS __attribute__ ((packed)); ! 249: fractype fraction:FRACBITS __attribute__ ((packed)); ! 250: #else ! 251: fractype fraction:FRACBITS __attribute__ ((packed)); ! 252: unsigned int exp:EXPBITS __attribute__ ((packed)); ! 253: unsigned int sign:1 __attribute__ ((packed)); ! 254: #endif ! 255: } ! 256: bits; ! 257: } ! 258: FLO_union_type; ! 259: ! 260: ! 261: /* end of header */ ! 262: ! 263: /* IEEE "special" number predicates */ ! 264: ! 265: #ifdef NO_NANS ! 266: ! 267: #define nan() 0 ! 268: #define isnan(x) 0 ! 269: #define isinf(x) 0 ! 270: #else ! 271: ! 272: INLINE ! 273: static fp_number_type * ! 274: nan () ! 275: { ! 276: static fp_number_type thenan; ! 277: ! 278: return &thenan; ! 279: } ! 280: ! 281: INLINE ! 282: static int ! 283: isnan ( fp_number_type * x) ! 284: { ! 285: return x->class == CLASS_SNAN || x->class == CLASS_QNAN; ! 286: } ! 287: ! 288: INLINE ! 289: static int ! 290: isinf ( fp_number_type * x) ! 291: { ! 292: return x->class == CLASS_INFINITY; ! 293: } ! 294: ! 295: #endif ! 296: ! 297: INLINE ! 298: static int ! 299: iszero ( fp_number_type * x) ! 300: { ! 301: return x->class == CLASS_ZERO; ! 302: } ! 303: ! 304: INLINE ! 305: static void ! 306: flip_sign ( fp_number_type * x) ! 307: { ! 308: x->sign = !x->sign; ! 309: } ! 310: ! 311: static FLO_type ! 312: pack_d ( fp_number_type * src) ! 313: { ! 314: FLO_union_type dst; ! 315: fractype fraction = src->fraction.ll; /* wasn't unsigned before? */ ! 316: ! 317: dst.bits.sign = src->sign; ! 318: ! 319: if (isnan (src)) ! 320: { ! 321: dst.bits.exp = EXPMAX; ! 322: dst.bits.fraction = src->fraction.ll; ! 323: if (src->class == CLASS_QNAN || 1) ! 324: { ! 325: dst.bits.fraction |= QUIET_NAN; ! 326: } ! 327: } ! 328: else if (isinf (src)) ! 329: { ! 330: dst.bits.exp = EXPMAX; ! 331: dst.bits.fraction = 0; ! 332: } ! 333: else if (iszero (src)) ! 334: { ! 335: dst.bits.exp = 0; ! 336: dst.bits.fraction = 0; ! 337: } ! 338: else if (fraction == 0) ! 339: { ! 340: dst.value = 0; ! 341: } ! 342: else ! 343: { ! 344: if (src->normal_exp < NORMAL_EXPMIN) ! 345: { ! 346: /* This number's exponent is too low to fit into the bits ! 347: available in the number, so we'll store 0 in the exponent and ! 348: shift the fraction to the right to make up for it. */ ! 349: ! 350: int shift = NORMAL_EXPMIN - src->normal_exp; ! 351: ! 352: dst.bits.exp = 0; ! 353: ! 354: if (shift > FRAC_NBITS - NGARDS) ! 355: { ! 356: /* No point shifting, since it's more that 64 out. */ ! 357: fraction = 0; ! 358: } ! 359: else ! 360: { ! 361: /* Shift by the value */ ! 362: fraction >>= shift; ! 363: } ! 364: fraction >>= NGARDS; ! 365: dst.bits.fraction = fraction; ! 366: } ! 367: else if (src->normal_exp > EXPBIAS) ! 368: { ! 369: dst.bits.exp = EXPMAX; ! 370: dst.bits.fraction = 0; ! 371: } ! 372: else ! 373: { ! 374: dst.bits.exp = src->normal_exp + EXPBIAS; ! 375: /* IF the gard bits are the all zero, but the first, then we're ! 376: half way between two numbers, choose the one which makes the ! 377: lsb of the answer 0. */ ! 378: if ((fraction & GARDMASK) == GARDMSB) ! 379: { ! 380: if (fraction & (1 << NGARDS)) ! 381: fraction += GARDROUND + 1; ! 382: } ! 383: else ! 384: { ! 385: /* Add a one to the guards to round up */ ! 386: fraction += GARDROUND; ! 387: } ! 388: if (fraction >= IMPLICIT_2) ! 389: { ! 390: fraction >>= 1; ! 391: dst.bits.exp += 1; ! 392: } ! 393: fraction >>= NGARDS; ! 394: dst.bits.fraction = fraction; ! 395: } ! 396: } ! 397: return dst.value; ! 398: } ! 399: ! 400: static void ! 401: unpack_d (FLO_union_type * src, fp_number_type * dst) ! 402: { ! 403: fractype fraction = src->bits.fraction; ! 404: ! 405: dst->sign = src->bits.sign; ! 406: if (src->bits.exp == 0) ! 407: { ! 408: /* Hmm. Looks like 0 */ ! 409: if (fraction == 0) ! 410: { ! 411: /* tastes like zero */ ! 412: dst->class = CLASS_ZERO; ! 413: } ! 414: else ! 415: { ! 416: /* Zero exponent with non zero fraction - it's denormalized, ! 417: so there isn't a leading implicit one - we'll shift it so ! 418: it gets one. */ ! 419: dst->normal_exp = src->bits.exp - EXPBIAS + 1; ! 420: fraction <<= NGARDS; ! 421: ! 422: dst->class = CLASS_NUMBER; ! 423: #if 1 ! 424: while (fraction < IMPLICIT_1) ! 425: { ! 426: fraction <<= 1; ! 427: dst->normal_exp--; ! 428: } ! 429: #endif ! 430: dst->fraction.ll = fraction; ! 431: } ! 432: } ! 433: else if (src->bits.exp == EXPMAX) ! 434: { ! 435: /* Huge exponent*/ ! 436: if (fraction == 0) ! 437: { ! 438: /* Attatched to a zero fraction - means infinity */ ! 439: dst->class = CLASS_INFINITY; ! 440: } ! 441: else ! 442: { ! 443: /* Non zero fraction, means nan */ ! 444: if (dst->sign) ! 445: { ! 446: dst->class = CLASS_SNAN; ! 447: } ! 448: else ! 449: { ! 450: dst->class = CLASS_QNAN; ! 451: } ! 452: /* Keep the fraction part as the nan number */ ! 453: dst->fraction.ll = fraction; ! 454: } ! 455: } ! 456: else ! 457: { ! 458: /* Nothing strange about this number */ ! 459: dst->normal_exp = src->bits.exp - EXPBIAS; ! 460: dst->class = CLASS_NUMBER; ! 461: dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1; ! 462: } ! 463: } ! 464: ! 465: static fp_number_type * ! 466: _fpadd_parts (fp_number_type * a, ! 467: fp_number_type * b, ! 468: fp_number_type * tmp) ! 469: { ! 470: intfrac tfraction; ! 471: ! 472: /* Put commonly used fields in local variables. */ ! 473: int a_normal_exp; ! 474: int b_normal_exp; ! 475: fractype a_fraction; ! 476: fractype b_fraction; ! 477: ! 478: if (isnan (a)) ! 479: { ! 480: return a; ! 481: } ! 482: if (isnan (b)) ! 483: { ! 484: return b; ! 485: } ! 486: if (isinf (a)) ! 487: { ! 488: return a; ! 489: } ! 490: if (isinf (b)) ! 491: { ! 492: return b; ! 493: } ! 494: if (iszero (b)) ! 495: { ! 496: return a; ! 497: } ! 498: if (iszero (a)) ! 499: { ! 500: return b; ! 501: } ! 502: ! 503: /* Got two numbers. shift the smaller and increment the exponent till ! 504: they're the same */ ! 505: { ! 506: int diff; ! 507: ! 508: a_normal_exp = a->normal_exp; ! 509: b_normal_exp = b->normal_exp; ! 510: a_fraction = a->fraction.ll; ! 511: b_fraction = b->fraction.ll; ! 512: ! 513: diff = a_normal_exp - b_normal_exp; ! 514: ! 515: if (diff < 0) ! 516: diff = -diff; ! 517: if (diff < FRAC_NBITS) ! 518: { ! 519: /* ??? This does shifts one bit at a time. Optimize. */ ! 520: while (a_normal_exp > b_normal_exp) ! 521: { ! 522: b_normal_exp++; ! 523: LSHIFT (b_fraction); ! 524: } ! 525: while (b_normal_exp > a_normal_exp) ! 526: { ! 527: a_normal_exp++; ! 528: LSHIFT (a_fraction); ! 529: } ! 530: } ! 531: else ! 532: { ! 533: /* Somethings's up.. choose the biggest */ ! 534: if (a_normal_exp > b_normal_exp) ! 535: { ! 536: b_normal_exp = a_normal_exp; ! 537: b_fraction = 0; ! 538: } ! 539: else ! 540: { ! 541: a_normal_exp = b_normal_exp; ! 542: a_fraction = 0; ! 543: } ! 544: } ! 545: } ! 546: ! 547: if (a->sign != b->sign) ! 548: { ! 549: if (a->sign) ! 550: { ! 551: tfraction = -a_fraction + b_fraction; ! 552: } ! 553: else ! 554: { ! 555: tfraction = a_fraction - b_fraction; ! 556: } ! 557: if (tfraction > 0) ! 558: { ! 559: tmp->sign = 0; ! 560: tmp->normal_exp = a_normal_exp; ! 561: tmp->fraction.ll = tfraction; ! 562: } ! 563: else ! 564: { ! 565: tmp->sign = 1; ! 566: tmp->normal_exp = a_normal_exp; ! 567: tmp->fraction.ll = -tfraction; ! 568: } ! 569: /* and renomalize it */ ! 570: ! 571: while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll) ! 572: { ! 573: tmp->fraction.ll <<= 1; ! 574: tmp->normal_exp--; ! 575: } ! 576: } ! 577: else ! 578: { ! 579: tmp->sign = a->sign; ! 580: tmp->normal_exp = a_normal_exp; ! 581: tmp->fraction.ll = a_fraction + b_fraction; ! 582: } ! 583: tmp->class = CLASS_NUMBER; ! 584: /* Now the fraction is added, we have to shift down to renormalize the ! 585: number */ ! 586: ! 587: if (tmp->fraction.ll >= IMPLICIT_2) ! 588: { ! 589: LSHIFT (tmp->fraction.ll); ! 590: tmp->normal_exp++; ! 591: } ! 592: return tmp; ! 593: ! 594: } ! 595: ! 596: FLO_type ! 597: add (FLO_type arg_a, FLO_type arg_b) ! 598: { ! 599: fp_number_type a; ! 600: fp_number_type b; ! 601: fp_number_type tmp; ! 602: fp_number_type *res; ! 603: ! 604: unpack_d ((FLO_union_type *) & arg_a, &a); ! 605: unpack_d ((FLO_union_type *) & arg_b, &b); ! 606: ! 607: res = _fpadd_parts (&a, &b, &tmp); ! 608: ! 609: return pack_d (res); ! 610: } ! 611: ! 612: FLO_type ! 613: sub (FLO_type arg_a, FLO_type arg_b) ! 614: { ! 615: fp_number_type a; ! 616: fp_number_type b; ! 617: fp_number_type tmp; ! 618: fp_number_type *res; ! 619: ! 620: unpack_d ((FLO_union_type *) & arg_a, &a); ! 621: unpack_d ((FLO_union_type *) & arg_b, &b); ! 622: ! 623: b.sign ^= 1; ! 624: ! 625: res = _fpadd_parts (&a, &b, &tmp); ! 626: ! 627: return pack_d (res); ! 628: } ! 629: ! 630: static fp_number_type * ! 631: _fpmul_parts ( fp_number_type * a, ! 632: fp_number_type * b, ! 633: fp_number_type * tmp) ! 634: { ! 635: fractype low = 0; ! 636: fractype high = 0; ! 637: ! 638: if (isnan (a)) ! 639: { ! 640: a->sign = a->sign != b->sign; ! 641: return a; ! 642: } ! 643: if (isnan (b)) ! 644: { ! 645: b->sign = a->sign != b->sign; ! 646: return b; ! 647: } ! 648: if (isinf (a)) ! 649: { ! 650: if (iszero (b)) ! 651: return nan (); ! 652: a->sign = a->sign != b->sign; ! 653: return a; ! 654: } ! 655: if (isinf (b)) ! 656: { ! 657: if (iszero (a)) ! 658: { ! 659: return nan (); ! 660: } ! 661: b->sign = a->sign != b->sign; ! 662: return b; ! 663: } ! 664: if (iszero (a)) ! 665: { ! 666: a->sign = a->sign != b->sign; ! 667: return a; ! 668: } ! 669: if (iszero (b)) ! 670: { ! 671: b->sign = a->sign != b->sign; ! 672: return b; ! 673: } ! 674: ! 675: /* Calculate the mantissa by multiplying both 64bit numbers to get a ! 676: 128 bit number */ ! 677: { ! 678: fractype x = a->fraction.ll; ! 679: fractype ylow = b->fraction.ll; ! 680: fractype yhigh = 0; ! 681: int bit; ! 682: ! 683: #if defined(NO_DI_MODE) ! 684: { ! 685: /* ??? This does multiplies one bit at a time. Optimize. */ ! 686: for (bit = 0; bit < FRAC_NBITS; bit++) ! 687: { ! 688: int carry; ! 689: ! 690: if (x & 1) ! 691: { ! 692: carry = (low += ylow) < ylow; ! 693: high += yhigh + carry; ! 694: } ! 695: yhigh <<= 1; ! 696: if (ylow & FRACHIGH) ! 697: { ! 698: yhigh |= 1; ! 699: } ! 700: ylow <<= 1; ! 701: x >>= 1; ! 702: } ! 703: } ! 704: #elif defined(FLOAT) ! 705: { ! 706: /* Multiplying two 32 bit numbers to get a 64 bit number on ! 707: a machine with DI, so we're safe */ ! 708: ! 709: DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll); ! 710: ! 711: high = answer >> 32; ! 712: low = answer; ! 713: } ! 714: #else ! 715: /* Doing a 64*64 to 128 */ ! 716: { ! 717: UDItype nl = a->fraction.ll & 0xffffffff; ! 718: UDItype nh = a->fraction.ll >> 32; ! 719: UDItype ml = b->fraction.ll & 0xffffffff; ! 720: UDItype mh = b->fraction.ll >>32; ! 721: UDItype pp_ll = ml * nl; ! 722: UDItype pp_hl = mh * nl; ! 723: UDItype pp_lh = ml * nh; ! 724: UDItype pp_hh = mh * nh; ! 725: UDItype res2 = 0; ! 726: UDItype res0 = 0; ! 727: UDItype ps_hh__ = pp_hl + pp_lh; ! 728: if (ps_hh__ < pp_hl) ! 729: res2 += 0x100000000LL; ! 730: pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL; ! 731: res0 = pp_ll + pp_hl; ! 732: if (res0 < pp_ll) ! 733: res2++; ! 734: res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh; ! 735: high = res2; ! 736: low = res0; ! 737: } ! 738: #endif ! 739: } ! 740: ! 741: tmp->normal_exp = a->normal_exp + b->normal_exp; ! 742: tmp->sign = a->sign != b->sign; ! 743: #ifdef FLOAT ! 744: tmp->normal_exp += 2; /* ??????????????? */ ! 745: #else ! 746: tmp->normal_exp += 4; /* ??????????????? */ ! 747: #endif ! 748: while (high >= IMPLICIT_2) ! 749: { ! 750: tmp->normal_exp++; ! 751: if (high & 1) ! 752: { ! 753: low >>= 1; ! 754: low |= FRACHIGH; ! 755: } ! 756: high >>= 1; ! 757: } ! 758: while (high < IMPLICIT_1) ! 759: { ! 760: tmp->normal_exp--; ! 761: ! 762: high <<= 1; ! 763: if (low & FRACHIGH) ! 764: high |= 1; ! 765: low <<= 1; ! 766: } ! 767: /* rounding is tricky. if we only round if it won't make us round later. */ ! 768: #if 0 ! 769: if (low & FRACHIGH2) ! 770: { ! 771: if (((high & GARDMASK) != GARDMSB) ! 772: && (((high + 1) & GARDMASK) == GARDMSB)) ! 773: { ! 774: /* don't round, it gets done again later. */ ! 775: } ! 776: else ! 777: { ! 778: high++; ! 779: } ! 780: } ! 781: #endif ! 782: if ((high & GARDMASK) == GARDMSB) ! 783: { ! 784: if (high & (1 << NGARDS)) ! 785: { ! 786: /* half way, so round to even */ ! 787: high += GARDROUND + 1; ! 788: } ! 789: else if (low) ! 790: { ! 791: /* but we really weren't half way */ ! 792: high += GARDROUND + 1; ! 793: } ! 794: } ! 795: tmp->fraction.ll = high; ! 796: tmp->class = CLASS_NUMBER; ! 797: return tmp; ! 798: } ! 799: ! 800: FLO_type ! 801: multiply (FLO_type arg_a, FLO_type arg_b) ! 802: { ! 803: fp_number_type a; ! 804: fp_number_type b; ! 805: fp_number_type tmp; ! 806: fp_number_type *res; ! 807: ! 808: unpack_d ((FLO_union_type *) & arg_a, &a); ! 809: unpack_d ((FLO_union_type *) & arg_b, &b); ! 810: ! 811: res = _fpmul_parts (&a, &b, &tmp); ! 812: ! 813: return pack_d (res); ! 814: } ! 815: ! 816: static fp_number_type * ! 817: _fpdiv_parts (fp_number_type * a, ! 818: fp_number_type * b, ! 819: fp_number_type * tmp) ! 820: { ! 821: fractype low = 0; ! 822: fractype high = 0; ! 823: fractype r0, r1, y0, y1, bit; ! 824: fractype q; ! 825: fractype numerator; ! 826: fractype denominator; ! 827: fractype quotient; ! 828: fractype remainder; ! 829: ! 830: if (isnan (a)) ! 831: { ! 832: return a; ! 833: } ! 834: if (isnan (b)) ! 835: { ! 836: return b; ! 837: } ! 838: if (isinf (a) || iszero (a)) ! 839: { ! 840: if (a->class == b->class) ! 841: return nan (); ! 842: return a; ! 843: } ! 844: a->sign = a->sign ^ b->sign; ! 845: ! 846: if (isinf (b)) ! 847: { ! 848: a->fraction.ll = 0; ! 849: a->normal_exp = 0; ! 850: return a; ! 851: } ! 852: if (iszero (b)) ! 853: { ! 854: a->class = CLASS_INFINITY; ! 855: return b; ! 856: } ! 857: ! 858: /* Calculate the mantissa by multiplying both 64bit numbers to get a ! 859: 128 bit number */ ! 860: { ! 861: int carry; ! 862: intfrac d0, d1; /* weren't unsigned before ??? */ ! 863: ! 864: /* quotient = ! 865: ( numerator / denominator) * 2^(numerator exponent - denominator exponent) ! 866: */ ! 867: ! 868: a->normal_exp = a->normal_exp - b->normal_exp; ! 869: numerator = a->fraction.ll; ! 870: denominator = b->fraction.ll; ! 871: ! 872: if (numerator < denominator) ! 873: { ! 874: /* Fraction will be less than 1.0 */ ! 875: numerator *= 2; ! 876: a->normal_exp--; ! 877: } ! 878: bit = IMPLICIT_1; ! 879: quotient = 0; ! 880: /* ??? Does divide one bit at a time. Optimize. */ ! 881: while (bit) ! 882: { ! 883: if (numerator >= denominator) ! 884: { ! 885: quotient |= bit; ! 886: numerator -= denominator; ! 887: } ! 888: bit >>= 1; ! 889: numerator *= 2; ! 890: } ! 891: ! 892: if ((quotient & GARDMASK) == GARDMSB) ! 893: { ! 894: if (quotient & (1 << NGARDS)) ! 895: { ! 896: /* half way, so round to even */ ! 897: quotient += GARDROUND + 1; ! 898: } ! 899: else if (numerator) ! 900: { ! 901: /* but we really weren't half way, more bits exist */ ! 902: quotient += GARDROUND + 1; ! 903: } ! 904: } ! 905: ! 906: a->fraction.ll = quotient; ! 907: return (a); ! 908: } ! 909: } ! 910: ! 911: FLO_type ! 912: divide (FLO_type arg_a, FLO_type arg_b) ! 913: { ! 914: fp_number_type a; ! 915: fp_number_type b; ! 916: fp_number_type tmp; ! 917: fp_number_type *res; ! 918: ! 919: unpack_d ((FLO_union_type *) & arg_a, &a); ! 920: unpack_d ((FLO_union_type *) & arg_b, &b); ! 921: ! 922: res = _fpdiv_parts (&a, &b, &tmp); ! 923: ! 924: return pack_d (res); ! 925: } ! 926: ! 927: /* according to the demo, fpcmp returns a comparison with 0... thus ! 928: a<b -> -1 ! 929: a==b -> 0 ! 930: a>b -> +1 ! 931: */ ! 932: ! 933: static int ! 934: _fpcmp_parts (fp_number_type * a, fp_number_type * b) ! 935: { ! 936: #if 0 ! 937: /* either nan -> unordered. Must be checked outside of this routine. */ ! 938: if (isnan (a) && isnan (b)) ! 939: { ! 940: return 1; /* still unordered! */ ! 941: } ! 942: #endif ! 943: ! 944: if (isnan (a) || isnan (b)) ! 945: { ! 946: return 1; /* how to indicate unordered compare? */ ! 947: } ! 948: if (isinf (a) && isinf (b)) ! 949: { ! 950: /* +inf > -inf, but +inf != +inf */ ! 951: /* b \a| +inf(0)| -inf(1) ! 952: ______\+--------+-------- ! 953: +inf(0)| a==b(0)| a<b(-1) ! 954: -------+--------+-------- ! 955: -inf(1)| a>b(1) | a==b(0) ! 956: -------+--------+-------- ! 957: So since unordered must be non zero, just line up the columns... ! 958: */ ! 959: return b->sign - a->sign; ! 960: } ! 961: /* but not both... */ ! 962: if (isinf (a)) ! 963: { ! 964: return a->sign ? -1 : 1; ! 965: } ! 966: if (isinf (b)) ! 967: { ! 968: return b->sign ? 1 : -1; ! 969: } ! 970: if (iszero (a) && iszero (b)) ! 971: { ! 972: return 0; ! 973: } ! 974: if (iszero (a)) ! 975: { ! 976: return b->sign ? 1 : -1; ! 977: } ! 978: if (iszero (b)) ! 979: { ! 980: return a->sign ? -1 : 1; ! 981: } ! 982: /* now both are "normal". */ ! 983: if (a->sign != b->sign) ! 984: { ! 985: /* opposite signs */ ! 986: return a->sign ? -1 : 1; ! 987: } ! 988: /* same sign; exponents? */ ! 989: if (a->normal_exp > b->normal_exp) ! 990: { ! 991: return a->sign ? -1 : 1; ! 992: } ! 993: if (a->normal_exp < b->normal_exp) ! 994: { ! 995: return a->sign ? 1 : -1; ! 996: } ! 997: /* same exponents; check size. */ ! 998: if (a->fraction.ll > b->fraction.ll) ! 999: { ! 1000: return a->sign ? -1 : 1; ! 1001: } ! 1002: if (a->fraction.ll < b->fraction.ll) ! 1003: { ! 1004: return a->sign ? 1 : -1; ! 1005: } ! 1006: /* after all that, they're equal. */ ! 1007: return 0; ! 1008: } ! 1009: ! 1010: CMPtype ! 1011: compare (FLO_type arg_a, FLO_type arg_b) ! 1012: { ! 1013: fp_number_type a; ! 1014: fp_number_type b; ! 1015: ! 1016: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1017: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1018: ! 1019: return _fpcmp_parts (&a, &b); ! 1020: } ! 1021: ! 1022: #ifndef US_SOFTWARE_GOFAST ! 1023: ! 1024: /* These should be optimized for their specific tasks someday. */ ! 1025: ! 1026: CMPtype ! 1027: _eq_f2 (FLO_type arg_a, FLO_type arg_b) ! 1028: { ! 1029: fp_number_type a; ! 1030: fp_number_type b; ! 1031: ! 1032: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1033: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1034: ! 1035: if (isnan (&a) || isnan (&b)) ! 1036: return 1; /* false, truth == 0 */ ! 1037: ! 1038: return _fpcmp_parts (&a, &b) ; ! 1039: } ! 1040: ! 1041: CMPtype ! 1042: _ne_f2 (FLO_type arg_a, FLO_type arg_b) ! 1043: { ! 1044: fp_number_type a; ! 1045: fp_number_type b; ! 1046: ! 1047: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1048: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1049: ! 1050: if (isnan (&a) || isnan (&b)) ! 1051: return 1; /* true, truth != 0 */ ! 1052: ! 1053: return _fpcmp_parts (&a, &b) ; ! 1054: } ! 1055: ! 1056: CMPtype ! 1057: _gt_f2 (FLO_type arg_a, FLO_type arg_b) ! 1058: { ! 1059: fp_number_type a; ! 1060: fp_number_type b; ! 1061: ! 1062: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1063: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1064: ! 1065: if (isnan (&a) || isnan (&b)) ! 1066: return -1; /* false, truth > 0 */ ! 1067: ! 1068: return _fpcmp_parts (&a, &b); ! 1069: } ! 1070: ! 1071: CMPtype ! 1072: _ge_f2 (FLO_type arg_a, FLO_type arg_b) ! 1073: { ! 1074: fp_number_type a; ! 1075: fp_number_type b; ! 1076: ! 1077: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1078: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1079: ! 1080: if (isnan (&a) || isnan (&b)) ! 1081: return -1; /* false, truth >= 0 */ ! 1082: return _fpcmp_parts (&a, &b) ; ! 1083: } ! 1084: ! 1085: CMPtype ! 1086: _lt_f2 (FLO_type arg_a, FLO_type arg_b) ! 1087: { ! 1088: fp_number_type a; ! 1089: fp_number_type b; ! 1090: ! 1091: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1092: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1093: ! 1094: if (isnan (&a) || isnan (&b)) ! 1095: return 1; /* false, truth < 0 */ ! 1096: ! 1097: return _fpcmp_parts (&a, &b); ! 1098: } ! 1099: ! 1100: CMPtype ! 1101: _le_f2 (FLO_type arg_a, FLO_type arg_b) ! 1102: { ! 1103: fp_number_type a; ! 1104: fp_number_type b; ! 1105: ! 1106: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1107: unpack_d ((FLO_union_type *) & arg_b, &b); ! 1108: ! 1109: if (isnan (&a) || isnan (&b)) ! 1110: return 1; /* false, truth <= 0 */ ! 1111: ! 1112: return _fpcmp_parts (&a, &b) ; ! 1113: } ! 1114: ! 1115: #endif /* ! US_SOFTWARE_GOFAST */ ! 1116: ! 1117: FLO_type ! 1118: si_to_float (SItype arg_a) ! 1119: { ! 1120: fp_number_type in; ! 1121: ! 1122: in.class = CLASS_NUMBER; ! 1123: in.sign = arg_a < 0; ! 1124: if (!arg_a) ! 1125: { ! 1126: in.class = CLASS_ZERO; ! 1127: } ! 1128: else ! 1129: { ! 1130: in.normal_exp = FRACBITS + NGARDS; ! 1131: if (in.sign) ! 1132: { ! 1133: /* Special case for minint, since there is no +ve integer ! 1134: representation for it */ ! 1135: if (arg_a == 0x80000000) ! 1136: { ! 1137: return -2147483648.0; ! 1138: } ! 1139: in.fraction.ll = (-arg_a); ! 1140: } ! 1141: else ! 1142: in.fraction.ll = arg_a; ! 1143: ! 1144: while (in.fraction.ll < (1LL << (FRACBITS + NGARDS))) ! 1145: { ! 1146: in.fraction.ll <<= 1; ! 1147: in.normal_exp -= 1; ! 1148: } ! 1149: } ! 1150: return pack_d (&in); ! 1151: } ! 1152: ! 1153: SItype ! 1154: float_to_si (FLO_type arg_a) ! 1155: { ! 1156: fp_number_type a; ! 1157: SItype tmp; ! 1158: ! 1159: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1160: if (iszero (&a)) ! 1161: return 0; ! 1162: if (isnan (&a)) ! 1163: return 0; ! 1164: /* get reasonable MAX_SI_INT... */ ! 1165: if (isinf (&a)) ! 1166: return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1; ! 1167: /* it is a number, but a small one */ ! 1168: if (a.normal_exp < 0) ! 1169: return 0; ! 1170: if (a.normal_exp > 30) ! 1171: return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT; ! 1172: tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp); ! 1173: return a.sign ? (-tmp) : (tmp); ! 1174: } ! 1175: ! 1176: #ifdef US_SOFTWARE_GOFAST ! 1177: /* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines, ! 1178: we also define them for GOFAST because the ones in libgcc2.c have the ! 1179: wrong names and I'd rather define these here and keep GOFAST CYG-LOC's ! 1180: out of libgcc2.c. We can't define these here if not GOFAST because then ! 1181: there'd be duplicate copies. */ ! 1182: ! 1183: USItype ! 1184: float_to_usi (FLO_type arg_a) ! 1185: { ! 1186: fp_number_type a; ! 1187: ! 1188: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1189: if (iszero (&a)) ! 1190: return 0; ! 1191: if (isnan (&a)) ! 1192: return 0; ! 1193: /* get reasonable MAX_USI_INT... */ ! 1194: if (isinf (&a)) ! 1195: return a.sign ? MAX_USI_INT : 0; ! 1196: /* it is a negative number */ ! 1197: if (a.sign) ! 1198: return 0; ! 1199: /* it is a number, but a small one */ ! 1200: if (a.normal_exp < 0) ! 1201: return 0; ! 1202: if (a.normal_exp > 31) ! 1203: return MAX_USI_INT; ! 1204: else if (a.normal_exp > (FRACBITS + NGARDS)) ! 1205: return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp); ! 1206: else ! 1207: return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp); ! 1208: } ! 1209: #endif ! 1210: ! 1211: FLO_type ! 1212: negate (FLO_type arg_a) ! 1213: { ! 1214: fp_number_type a; ! 1215: ! 1216: unpack_d ((FLO_union_type *) & arg_a, &a); ! 1217: flip_sign (&a); ! 1218: return pack_d (&a); ! 1219: } ! 1220: ! 1221: #ifdef FLOAT ! 1222: ! 1223: SFtype ! 1224: __make_fp(fp_class_type class, ! 1225: unsigned int sign, ! 1226: int exp, ! 1227: USItype frac) ! 1228: { ! 1229: fp_number_type in; ! 1230: ! 1231: in.class = class; ! 1232: in.sign = sign; ! 1233: in.normal_exp = exp; ! 1234: in.fraction.ll = frac; ! 1235: return pack_d (&in); ! 1236: } ! 1237: ! 1238: #ifndef FLOAT_ONLY ! 1239: ! 1240: /* This enables one to build an fp library that supports float but not double. ! 1241: Otherwise, we would get an undefined reference to __make_dp. ! 1242: This is needed for some 8-bit ports that can't handle well values that ! 1243: are 8-bytes in size, so we just don't support double for them at all. */ ! 1244: ! 1245: extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac); ! 1246: ! 1247: DFtype ! 1248: sf_to_df (SFtype arg_a) ! 1249: { ! 1250: fp_number_type in; ! 1251: ! 1252: unpack_d ((FLO_union_type *) & arg_a, &in); ! 1253: return __make_dp (in.class, in.sign, in.normal_exp, ! 1254: ((UDItype) in.fraction.ll) << F_D_BITOFF); ! 1255: } ! 1256: ! 1257: #endif ! 1258: #endif ! 1259: ! 1260: #ifndef FLOAT ! 1261: ! 1262: extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype); ! 1263: ! 1264: DFtype ! 1265: __make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac) ! 1266: { ! 1267: fp_number_type in; ! 1268: ! 1269: in.class = class; ! 1270: in.sign = sign; ! 1271: in.normal_exp = exp; ! 1272: in.fraction.ll = frac; ! 1273: return pack_d (&in); ! 1274: } ! 1275: ! 1276: SFtype ! 1277: df_to_sf (DFtype arg_a) ! 1278: { ! 1279: fp_number_type in; ! 1280: ! 1281: unpack_d ((FLO_union_type *) & arg_a, &in); ! 1282: return __make_fp (in.class, in.sign, in.normal_exp, ! 1283: in.fraction.ll >> F_D_BITOFF); ! 1284: } ! 1285: ! 1286: #endif
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.