gcc/floatlib.c

   1 /*
   2 ** libgcc support for software floating point.
   3 ** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
   4 ** Permission is granted to do *anything* you want with this file,
   5 ** commercial or otherwise, provided this message remains intact.  So there!
   6 ** I would appreciate receiving any updates/patches/changes that anyone
   7 ** makes, and am willing to be the repository for said changes (am I
   8 ** making a big mistake?).
   9
  10 Warning! Only single-precision is actually implemented.  This file
  11 won't really be much use until double-precision is supported.
  12
  13 However, once that is done, this file might eventually become a
  14 replacement for libgcc1.c.  It might also make possible
  15 cross-compilation for an IEEE target machine from a non-IEEE
  16 host such as a VAX.
  17
  18 If you'd like to work on completing this, please talk to rms@gnu.ai.mit.edu.
  19
  20 --> Double precision floating support added by James Carlson on 20 April 1998.
  21
  22 **
  23 ** Pat Wood
  24 ** Pipeline Associates, Inc.
  25 ** pipeline!phw@motown.com or
  26 ** sun!pipeline!phw or
  27 ** uunet!motown!pipeline!phw
  28 **
  29 ** 05/01/91 -- V1.0 -- first release to gcc mailing lists
  30 ** 05/04/91 -- V1.1 -- added float and double prototypes and return values
  31 **                  -- fixed problems with adding and subtracting zero
  32 **                  -- fixed rounding in truncdfsf2
  33 **                  -- fixed SWAP define and tested on 386
  34 */
  35
  36 /*
  37 ** The following are routines that replace the libgcc soft floating point
  38 ** routines that are called automatically when -msoft-float is selected.
  39 ** The support single and double precision IEEE format, with provisions
  40 ** for byte-swapped machines (tested on 386).  Some of the double-precision
  41 ** routines work at full precision, but most of the hard ones simply punt
  42 ** and call the single precision routines, producing a loss of accuracy.
  43 ** long long support is not assumed or included.
  44 ** Overall accuracy is close to IEEE (actually 68882) for single-precision
  45 ** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
  46 ** being rounded the wrong way during a multiply.  I'm not fussy enough to
  47 ** bother with it, but if anyone is, knock yourself out.
  48 **
  49 ** Efficiency has only been addressed where it was obvious that something
  50 ** would make a big difference.  Anyone who wants to do this right for
  51 ** best speed should go in and rewrite in assembler.
  52 **
  53 ** I have tested this only on a 68030 workstation and 386/ix integrated
  54 ** in with -msoft-float.
  55 */
  56
  57 /* the following deal with IEEE single-precision numbers */
  58 #define EXCESS          126
  59 #define SIGNBIT         0x80000000
  60 #define HIDDEN          (1 << 23)
  61 #define SIGN(fp)        ((fp) & SIGNBIT)
  62 #define EXP(fp)         (((fp) >> 23) & 0xFF)
  63 #define MANT(fp)        (((fp) & 0x7FFFFF) | HIDDEN)
  64 #define PACK(s,e,m)     ((s) | ((e) << 23) | (m))
  65
  66 /* the following deal with IEEE double-precision numbers */
  67 #define EXCESSD         1022
  68 #define HIDDEND         (1 << 20)
  69 #define EXPD(fp)        (((fp.l.upper) >> 20) & 0x7FF)
  70 #define SIGND(fp)       ((fp.l.upper) & SIGNBIT)
  71 #define MANTD(fp)       (((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
  72                                 (fp.l.lower >> 22))
  73 #define HIDDEND_LL      ((long long)1 << 52)
  74 #define MANTD_LL(fp)    ((fp.ll & (HIDDEND_LL-1)) | HIDDEND_LL)
  75 #define PACKD_LL(s,e,m) (((long long)((s)+((e)<<20))<<32)|(m))
  76
  77 /* define SWAP for 386/960 reverse-byte-order brain-damaged CPUs */
  78 union double_long {
  79     double d;
  80 #ifdef SWAP
  81     struct {
  82       unsigned long lower;
  83       long upper;
  84     } l;
  85 #else
  86     struct {
  87       long upper;
  88       unsigned long lower;
  89     } l;
  90 #endif
  91     long long ll;
  92 };
  93
  94 union float_long
  95   {
  96     float f;
  97     long l;
  98   };
  99
 100 /* add two floats */
 101 float
 102 __addsf3 (float a1, float a2)
 103 {
 104   register long mant1, mant2;
 105   register union float_long fl1, fl2;
 106   register int exp1, exp2;
 107   int sign = 0;
 108
 109   fl1.f = a1;
 110   fl2.f = a2;
 111
 112   /* check for zero args */
 113   if (!fl1.l) {
 114       fl1.f = fl2.f;
 115       goto test_done;
 116   }
 117   if (!fl2.l)
 118       goto test_done;
 119
 120   exp1 = EXP (fl1.l);
 121   exp2 = EXP (fl2.l);
 122
 123   if (exp1 > exp2 + 25)
 124       goto test_done;
 125   if (exp2 > exp1 + 25) {
 126       fl1.f = fl2.f;
 127       goto test_done;
 128   }
 129
 130   /* do everything in excess precision so's we can round later */
 131   mant1 = MANT (fl1.l) << 6;
 132   mant2 = MANT (fl2.l) << 6;
 133
 134   if (SIGN (fl1.l))
 135     mant1 = -mant1;
 136   if (SIGN (fl2.l))
 137     mant2 = -mant2;
 138
 139   if (exp1 > exp2)
 140     {
 141       mant2 >>= exp1 - exp2;
 142     }
 143   else
 144     {
 145       mant1 >>= exp2 - exp1;
 146       exp1 = exp2;
 147     }
 148   mant1 += mant2;
 149
 150   if (mant1 < 0)
 151     {
 152       mant1 = -mant1;
 153       sign = SIGNBIT;
 154     }
 155   else if (!mant1) {
 156       fl1.f = 0;
 157       goto test_done;
 158   }
 159
 160   /* normalize up */
 161   while (!(mant1 & 0xE0000000))
 162     {
 163       mant1 <<= 1;
 164       exp1--;
 165     }
 166
 167   /* normalize down? */
 168   if (mant1 & (1 << 30))
 169     {
 170       mant1 >>= 1;
 171       exp1++;
 172     }
 173
 174   /* round to even */
 175   mant1 += (mant1 & 0x40) ? 0x20 : 0x1F;
 176
 177   /* normalize down? */
 178   if (mant1 & (1 << 30))
 179     {
 180       mant1 >>= 1;
 181       exp1++;
 182     }
 183
 184   /* lose extra precision */
 185   mant1 >>= 6;
 186
 187   /* turn off hidden bit */
 188   mant1 &= ~HIDDEN;
 189
 190   /* pack up and go home */
 191   fl1.l = PACK (sign, exp1, mant1);
 192 test_done:
 193   return (fl1.f);
 194 }
 195
 196 /* subtract two floats */
 197 float
 198 __subsf3 (float a1, float a2)
 199 {
 200   register union float_long fl1, fl2;
 201
 202   fl1.f = a1;
 203   fl2.f = a2;
 204
 205   /* check for zero args */
 206   if (!fl2.l)
 207     return (fl1.f);
 208   if (!fl1.l)
 209     return (-fl2.f);
 210
 211   /* twiddle sign bit and add */
 212   fl2.l ^= SIGNBIT;
 213   return __addsf3 (a1, fl2.f);
 214 }
 215
 216 /* compare two floats */
 217 long
 218 __cmpsf2 (float a1, float a2)
 219 {
 220   register union float_long fl1, fl2;
 221
 222   fl1.f = a1;
 223   fl2.f = a2;
 224
 225   if (SIGN (fl1.l) && SIGN (fl2.l))
 226     {
 227       fl1.l ^= SIGNBIT;
 228       fl2.l ^= SIGNBIT;
 229     }
 230   if (fl1.l < fl2.l)
 231     return (-1);
 232   if (fl1.l > fl2.l)
 233     return (1);
 234   return (0);
 235 }
 236
 237 /* multiply two floats */
 238 float
 239 __mulsf3 (float a1, float a2)
 240 {
 241   register union float_long fl1, fl2;
 242   register unsigned long result;
 243   register int exp;
 244   int sign;
 245
 246   fl1.f = a1;
 247   fl2.f = a2;
 248
 249   if (!fl1.l || !fl2.l) {
 250       fl1.f = 0;
 251       goto test_done;
 252   }
 253
 254   /* compute sign and exponent */
 255   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 256   exp = EXP (fl1.l) - EXCESS;
 257   exp += EXP (fl2.l);
 258
 259   fl1.l = MANT (fl1.l);
 260   fl2.l = MANT (fl2.l);
 261
 262   /* the multiply is done as one 16x16 multiply and two 16x8 multiples */
 263   result = (fl1.l >> 8) * (fl2.l >> 8);
 264   result += ((fl1.l & 0xFF) * (fl2.l >> 8)) >> 8;
 265   result += ((fl2.l & 0xFF) * (fl1.l >> 8)) >> 8;
 266
 267   result >>= 2;
 268   if (result & 0x20000000)
 269     {
 270       /* round */
 271       result += 0x20;
 272       result >>= 6;
 273     }
 274   else
 275     {
 276       /* round */
 277       result += 0x10;
 278       result >>= 5;
 279       exp--;
 280     }
 281   if (result & (HIDDEN<<1)) {
 282     result >>= 1;
 283     exp++;
 284   }
 285
 286   result &= ~HIDDEN;
 287
 288   /* pack up and go home */
 289   fl1.l = PACK (sign, exp, result);
 290 test_done:
 291   return (fl1.f);
 292 }
 293
 294 /* divide two floats */
 295 float
 296 __divsf3 (float a1, float a2)
 297 {
 298   register union float_long fl1, fl2;
 299   register int result;
 300   register int mask;
 301   register int exp, sign;
 302
 303   fl1.f = a1;
 304   fl2.f = a2;
 305
 306   /* subtract exponents */
 307   exp = EXP (fl1.l) - EXP (fl2.l) + EXCESS;
 308
 309   /* compute sign */
 310   sign = SIGN (fl1.l) ^ SIGN (fl2.l);
 311
 312   /* divide by zero??? */
 313   if (!fl2.l)
 314     /* return NaN or -NaN */
 315     return (sign ? 0xFFFFFFFF : 0x7FFFFFFF);
 316
 317   /* numerator zero??? */
 318   if (!fl1.l)
 319     return (0);
 320
 321   /* now get mantissas */
 322   fl1.l = MANT (fl1.l);
 323   fl2.l = MANT (fl2.l);
 324
 325   /* this assures we have 25 bits of precision in the end */
 326   if (fl1.l < fl2.l)
 327     {
 328       fl1.l <<= 1;
 329       exp--;
 330     }
 331
 332   /* now we perform repeated subtraction of fl2.l from fl1.l */
 333   mask = 0x1000000;
 334   result = 0;
 335   while (mask)
 336     {
 337       if (fl1.l >= fl2.l)
 338         {
 339           result |= mask;
 340           fl1.l -= fl2.l;
 341         }
 342       fl1.l <<= 1;
 343       mask >>= 1;
 344     }
 345
 346   /* round */
 347   result += 1;
 348
 349   /* normalize down */
 350   exp++;
 351   result >>= 1;
 352
 353   result &= ~HIDDEN;
 354
 355   /* pack up and go home */
 356   fl1.l = PACK (sign, exp, result);
 357   return (fl1.f);
 358 }
 359
 360 /* convert int to double */
 361 double
 362 __floatsidf (register long a1)
 363 {
 364   register int sign = 0, exp = 31 + EXCESSD;
 365   union double_long dl;
 366
 367   if (!a1)
 368     {
 369       dl.l.upper = dl.l.lower = 0;
 370       return (dl.d);
 371     }
 372
 373   if (a1 < 0)
 374     {
 375       sign = SIGNBIT;
 376       a1 = -a1;
 377     }
 378
 379   while (a1 < 0x1000000)
 380     {
 381       a1 <<= 4;
 382       exp -= 4;
 383     }
 384
 385   while (a1 < 0x40000000)
 386     {
 387       a1 <<= 1;
 388       exp--;
 389     }
 390
 391   /* pack up and go home */
 392   dl.l.upper = sign;
 393   dl.l.upper |= exp << 20;
 394   dl.l.upper |= (a1 >> 10) & ~HIDDEND;
 395   dl.l.lower = a1 << 22;
 396
 397   return (dl.d);
 398 }
 399
 400 double
 401 __floatdidf (register long long a1)
 402 {
 403     register int exp = 63 + EXCESSD;
 404     union double_long dl;
 405
 406     dl.l.upper = dl.l.lower = 0;
 407     if (a1 == 0)
 408         return (dl.d);
 409
 410     if (a1 < 0) {
 411         dl.l.upper = SIGNBIT;
 412         a1 = -a1;
 413     }
 414
 415     while (a1 < (long long)1<<54) {
 416         a1 <<= 8;
 417         exp -= 8;
 418     }
 419     while (a1 < (long long)1<<62) {
 420         a1 <<= 1;
 421         exp -= 1;
 422     }
 423
 424   /* pack up and go home */
 425     dl.ll |= (a1 >> 10) & ~HIDDEND_LL;
 426     dl.l.upper |= exp << 20;
 427
 428     return (dl.d);
 429 }
 430
 431 float
 432 __floatsisf (register long a1)
 433 {
 434     (float)__floatsidf(a1);
 435 }
 436
 437 float
 438 __floatdisf (register long long a1)
 439 {
 440     (float)__floatdidf(a1);
 441 }
 442
 443 /* negate a float */
 444 float
 445 __negsf2 (float a1)
 446 {
 447   register union float_long fl1;
 448
 449   fl1.f = a1;
 450   if (!fl1.l)
 451     return (0);
 452
 453   fl1.l ^= SIGNBIT;
 454   return (fl1.f);
 455 }
 456
 457 /* negate a double */
 458 double
 459 __negdf2 (double a1)
 460 {
 461   register union double_long dl1;
 462
 463   dl1.d = a1;
 464
 465   if (!dl1.l.upper && !dl1.l.lower)
 466       return (dl1.d);
 467
 468   dl1.l.upper ^= SIGNBIT;
 469   return (dl1.d);
 470 }
 471
 472 /* convert float to double */
 473 double
 474 __extendsfdf2 (float a1)
 475 {
 476   register union float_long fl1;
 477   register union double_long dl;
 478   register int exp;
 479
 480   fl1.f = a1;
 481
 482   if (!fl1.l)
 483     {
 484       dl.l.upper = dl.l.lower = 0;
 485       return (dl.d);
 486     }
 487
 488   dl.l.upper = SIGN (fl1.l);
 489   exp = EXP (fl1.l) - EXCESS + EXCESSD;
 490   dl.l.upper |= exp << 20;
 491   dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
 492   dl.l.lower = MANT (fl1.l) << 29;
 493
 494   return (dl.d);
 495 }
 496
 497 /* convert double to float */
 498 float
 499 __truncdfsf2 (double a1)
 500 {
 501   register int exp;
 502   register long mant;
 503   register union float_long fl;
 504   register union double_long dl1;
 505
 506   dl1.d = a1;
 507
 508   if (!dl1.l.upper && !dl1.l.lower)
 509     return (float)(0);
 510
 511   exp = EXPD (dl1) - EXCESSD + EXCESS;
 512
 513   /* shift double mantissa 6 bits so we can round */
 514   mant = MANTD (dl1) >> 6;
 515
 516   /* now round and shift down */
 517   mant += 1;
 518   mant >>= 1;
 519
 520   /* did the round overflow? */
 521   if (mant & 0xFE000000)
 522     {
 523       mant >>= 1;
 524       exp++;
 525     }
 526
 527   mant &= ~HIDDEN;
 528
 529   /* pack up and go home */
 530   fl.l = PACK (SIGND (dl1), exp, mant);
 531   return (fl.f);
 532 }
 533
 534 /* compare two doubles */
 535 long
 536 __cmpdf2 (double a1, double a2)
 537 {
 538   register union double_long dl1, dl2;
 539
 540   dl1.d = a1;
 541   dl2.d = a2;
 542
 543   if (SIGND (dl1) && SIGND (dl2))
 544     {
 545       dl1.l.upper ^= SIGNBIT;
 546       dl2.l.upper ^= SIGNBIT;
 547     }
 548   if (dl1.l.upper < dl2.l.upper)
 549     return (-1);
 550   if (dl1.l.upper > dl2.l.upper)
 551     return (1);
 552   if (dl1.l.lower < dl2.l.lower)
 553     return (-1);
 554   if (dl1.l.lower > dl2.l.lower)
 555     return (1);
 556   return (0);
 557 }
 558
 559 /* convert double to int */
 560 long
 561 __fixdfsi (double a1)
 562 {
 563   register union double_long dl1;
 564   register int exp;
 565   register long l;
 566
 567   dl1.d = a1;
 568
 569   if (!dl1.l.upper && !dl1.l.lower)
 570     return (0);
 571
 572   exp = EXPD (dl1) - EXCESSD - 31;
 573   l = MANTD (dl1);
 574
 575   if (exp > 0)
 576       return SIGND(dl1) ? (1<<31) : ((1ul<<31)-1);
 577
 578   /* shift down until exp = 0 or l = 0 */
 579   if (exp < 0 && exp > -32 && l)
 580     l >>= -exp;
 581   else
 582     return (0);
 583
 584   return (SIGND (dl1) ? -l : l);
 585 }
 586
 587 /* convert double to int */
 588 long long
 589 __fixdfdi (double a1)
 590 {
 591     register union double_long dl1;
 592     register int exp;
 593     register long long l;
 594
 595     dl1.d = a1;
 596
 597     if (!dl1.l.upper && !dl1.l.lower)
 598         return (0);
 599
 600     exp = EXPD (dl1) - EXCESSD - 64;
 601     l = MANTD_LL(dl1);
 602
 603     if (exp > 0) {
 604         l = (long long)1<<63;
 605         if (!SIGND(dl1))
 606             l--;
 607         return l;
 608     }
 609
 610     /* shift down until exp = 0 or l = 0 */
 611     if (exp < 0 && exp > -64 && l)
 612         l >>= -exp;
 613     else
 614         return (0);
 615
 616     return (SIGND (dl1) ? -l : l);
 617 }
 618
 619 /* convert double to unsigned int */
 620 unsigned long
 621 __fixunsdfsi (double a1)
 622 {
 623   register union double_long dl1;
 624   register int exp;
 625   register unsigned long l;
 626
 627   dl1.d = a1;
 628
 629   if (!dl1.l.upper && !dl1.l.lower)
 630     return (0);
 631
 632   exp = EXPD (dl1) - EXCESSD - 32;
 633   l = (((((dl1.l.upper) & 0xFFFFF) | HIDDEND) << 11) | (dl1.l.lower >> 21));
 634
 635   if (exp > 0)
 636     return (0xFFFFFFFFul);      /* largest integer */
 637
 638   /* shift down until exp = 0 or l = 0 */
 639   if (exp < 0 && exp > -32 && l)
 640     l >>= -exp;
 641   else
 642     return (0);
 643
 644   return (l);
 645 }
 646
 647 /* convert double to unsigned int */
 648 unsigned long long
 649 __fixunsdfdi (double a1)
 650 {
 651     register union double_long dl1;
 652     register int exp;
 653     register unsigned long long l;
 654
 655     dl1.d = a1;
 656
 657     if (dl1.ll == 0)
 658         return (0);
 659
 660     exp = EXPD (dl1) - EXCESSD - 64;
 661
 662     l = dl1.ll;
 663
 664     if (exp > 0)
 665         return (unsigned long long)-1;
 666
 667     /* shift down until exp = 0 or l = 0 */
 668     if (exp < 0 && exp > -64 && l)
 669         l >>= -exp;
 670     else
 671         return (0);
 672
 673     return (l);
 674 }
 675
 676 /* addtwo doubles */
 677 double
 678 __adddf3 (double a1, double a2)
 679 {
 680     register long long mant1, mant2;
 681     register union double_long fl1, fl2;
 682     register int exp1, exp2;
 683     int sign = 0;
 684
 685     fl1.d = a1;
 686     fl2.d = a2;
 687
 688     /* check for zero args */
 689     if (!fl2.ll)
 690         goto test_done;
 691     if (!fl1.ll) {
 692         fl1.d = fl2.d;
 693         goto test_done;
 694     }
 695
 696     exp1 = EXPD(fl1);
 697     exp2 = EXPD(fl2);
 698
 699     if (exp1 > exp2 + 54)
 700         goto test_done;
 701     if (exp2 > exp1 + 54) {
 702         fl1.d = fl2.d;
 703         goto test_done;
 704     }
 705
 706     /* do everything in excess precision so's we can round later */
 707     mant1 = MANTD_LL(fl1) << 9;
 708     mant2 = MANTD_LL(fl2) << 9;
 709
 710     if (SIGND(fl1))
 711         mant1 = -mant1;
 712     if (SIGND(fl2))
 713         mant2 = -mant2;
 714
 715     if (exp1 > exp2)
 716         mant2 >>= exp1 - exp2;
 717     else {
 718         mant1 >>= exp2 - exp1;
 719         exp1 = exp2;
 720     }
 721     mant1 += mant2;
 722
 723     if (mant1 < 0) {
 724         mant1 = -mant1;
 725         sign = SIGNBIT;
 726     } else if (!mant1) {
 727         fl1.d = 0;
 728         goto test_done;
 729     }
 730
 731     /* normalize up */
 732     while (!(mant1 & ((long long)7<<61))) {
 733         mant1 <<= 1;
 734         exp1--;
 735     }
 736
 737     /* normalize down? */
 738     if (mant1 & ((long long)3<<62)) {
 739         mant1 >>= 1;
 740         exp1++;
 741     }
 742
 743     /* round to even */
 744     mant1 += (mant1 & (1<<9)) ? (1<<8) : ((1<<8)-1);
 745
 746     /* normalize down? */
 747     if (mant1 & ((long long)3<<62)) {
 748         mant1 >>= 1;
 749         exp1++;
 750     }
 751
 752     /* lose extra precision */
 753     mant1 >>= 9;
 754
 755     /* turn off hidden bit */
 756     mant1 &= ~HIDDEND_LL;
 757
 758     /* pack up and go home */
 759     fl1.ll = PACKD_LL(sign,exp1,mant1);
 760
 761 test_done:
 762     return (fl1.d);
 763 }
 764
 765 /* subtract two doubles */
 766 double
 767 __subdf3 (double a1, double a2)
 768 {
 769     register union double_long fl1, fl2;
 770
 771     fl1.d = a1;
 772     fl2.d = a2;
 773
 774     /* check for zero args */
 775     if (!fl2.ll)
 776         return (fl1.d);
 777     /* twiddle sign bit and add */
 778     fl2.l.upper ^= SIGNBIT;
 779     if (!fl1.ll)
 780         return (fl2.d);
 781     return __adddf3 (a1, fl2.d);
 782 }
 783
 784 /* multiply two doubles */
 785 double
 786 __muldf3 (double a1, double a2)
 787 {
 788     register union double_long fl1, fl2;
 789     register unsigned long long result;
 790     register int exp;
 791     int sign;
 792
 793     fl1.d = a1;
 794     fl2.d = a2;
 795
 796     if (!fl1.ll || !fl2.ll) {
 797         fl1.d = 0;
 798         goto test_done;
 799     }
 800
 801     /* compute sign and exponent */
 802     sign = SIGND(fl1) ^ SIGND(fl2);
 803     exp = EXPD(fl1) - EXCESSD;
 804     exp += EXPD(fl2);
 805
 806     fl1.ll = MANTD_LL(fl1);
 807     fl2.ll = MANTD_LL(fl2);
 808
 809   /* the multiply is done as one 31x31 multiply and two 31x21 multiples */
 810     result = (fl1.ll >> 21) * (fl2.ll >> 21);
 811     result += ((fl1.ll & 0x1FFFFF) * (fl2.ll >> 21)) >> 21;
 812     result += ((fl2.ll & 0x1FFFFF) * (fl1.ll >> 21)) >> 21;
 813
 814     result >>= 2;
 815     if (result & ((long long)1<<61)) {
 816         /* round */
 817         result += 1<<8;
 818         result >>= 9;
 819     } else {
 820         /* round */
 821         result += 1<<7;
 822         result >>= 8;
 823         exp--;
 824     }
 825     if (result & (HIDDEND_LL<<1)) {
 826         result >>= 1;
 827         exp++;
 828     }
 829
 830     result &= ~HIDDEND_LL;
 831
 832     /* pack up and go home */
 833     fl1.ll = PACKD_LL(sign,exp,result);
 834 test_done:
 835     return (fl1.d);
 836 }
 837
 838 /* divide two doubles */
 839 double
 840 __divdf3 (double a1, double a2)
 841 {
 842     register union double_long fl1, fl2;
 843     register long long mask,result;
 844     register int exp, sign;
 845
 846     fl1.d = a1;
 847     fl2.d = a2;
 848
 849     /* subtract exponents */
 850     exp = EXPD(fl1) - EXPD(fl2) + EXCESSD;
 851
 852     /* compute sign */
 853     sign = SIGND(fl1) ^ SIGND(fl2);
 854
 855     /* numerator zero??? */
 856     if (fl1.ll == 0) {
 857         /* divide by zero??? */
 858         if (fl2.ll == 0)
 859             fl1.ll = ((unsigned long long)1<<63)-1;     /* NaN */
 860         else
 861             fl1.ll = 0;
 862         goto test_done;
 863     }
 864
 865     /* return +Inf or -Inf */
 866     if (fl2.ll == 0) {
 867         fl1.ll = PACKD_LL(SIGND(fl1),2047,0);
 868         goto test_done;
 869     }
 870
 871
 872     /* now get mantissas */
 873     fl1.ll = MANTD_LL(fl1);
 874     fl2.ll = MANTD_LL(fl2);
 875
 876     /* this assures we have 54 bits of precision in the end */
 877     if (fl1.ll < fl2.ll) {
 878         fl1.ll <<= 1;
 879         exp--;
 880     }
 881
 882     /* now we perform repeated subtraction of fl2.ll from fl1.ll */
 883     mask = (long long)1<<53;
 884     result = 0;
 885     while (mask) {
 886         if (fl1.ll >= fl2.ll)
 887         {
 888             result |= mask;
 889             fl1.ll -= fl2.ll;
 890         }
 891         fl1.ll <<= 1;
 892         mask >>= 1;
 893     }
 894
 895     /* round */
 896     result += 1;
 897
 898     /* normalize down */
 899     exp++;
 900     result >>= 1;
 901
 902     result &= ~HIDDEND_LL;
 903
 904     /* pack up and go home */
 905     fl1.ll = PACKD_LL(sign, exp, result);
 906
 907 test_done:
 908     return (fl1.d);
 909 }
 910
 911 int
 912 __gtdf2 (double a1, double a2)
 913 {
 914     return __cmpdf2 ((float) a1, (float) a2) > 0;
 915 }
 916
 917 int
 918 __gedf2 (double a1, double a2)
 919 {
 920     return (__cmpdf2 ((float) a1, (float) a2) >= 0) - 1;
 921 }
 922
 923 int
 924 __ltdf2 (double a1, double a2)
 925 {
 926     return - (__cmpdf2 ((float) a1, (float) a2) < 0);
 927 }
 928
 929 int
 930 __ledf2 (double a1, double a2)
 931 {
 932     return __cmpdf2 ((float) a1, (float) a2) > 0;
 933 }
 934
 935 int
 936 __eqdf2 (double a1, double a2)
 937 {
 938     return *(long long *) &a1 == *(long long *) &a2;
 939 }
 940
 941 int
 942 __nedf2 (double a1, double a2)
 943 {
 944     return *(long long *) &a1 != *(long long *) &a2;
 945 }