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