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