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
__addsf3(float a1,float a2)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
__subsf3(float a1,float a2)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
__cmpsf2(float a1,float a2)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
__mulsf3(float a1,float a2)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
__divsf3(float a1,float a2)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
__floatsidf(long a1)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
__floatdidf(long long a1)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
__floatsisf(long a1)431 __floatsisf (long a1)
432 {
433 (float)__floatsidf(a1);
434 }
435
436 float
__floatdisf(long long a1)437 __floatdisf (long long a1)
438 {
439 (float)__floatdidf(a1);
440 }
441
442 /* negate a float */
443 float
__negsf2(float a1)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
__negdf2(double a1)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
__extendsfdf2(float a1)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
__truncdfsf2(double a1)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
__cmpdf2(double a1,double a2)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
__fixdfsi(double a1)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
__fixdfdi(double a1)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
__fixunsdfsi(double a1)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
__fixunsdfdi(double a1)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
__adddf3(double a1,double a2)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
__subdf3(double a1,double a2)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
__muldf3(double a1,double a2)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
__divdf3(double a1,double a2)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
__gtdf2(double a1,double a2)911 __gtdf2 (double a1, double a2)
912 {
913 return __cmpdf2 ((float) a1, (float) a2) > 0;
914 }
915
916 int
__gedf2(double a1,double a2)917 __gedf2 (double a1, double a2)
918 {
919 return (__cmpdf2 ((float) a1, (float) a2) >= 0) - 1;
920 }
921
922 int
__ltdf2(double a1,double a2)923 __ltdf2 (double a1, double a2)
924 {
925 return - (__cmpdf2 ((float) a1, (float) a2) < 0);
926 }
927
928 int
__ledf2(double a1,double a2)929 __ledf2 (double a1, double a2)
930 {
931 return __cmpdf2 ((float) a1, (float) a2) > 0;
932 }
933
934 int
__eqdf2(double a1,double a2)935 __eqdf2 (double a1, double a2)
936 {
937 return *(long long *) &a1 == *(long long *) &a2;
938 }
939
940 int
__nedf2(double a1,double a2)941 __nedf2 (double a1, double a2)
942 {
943 return *(long long *) &a1 != *(long long *) &a2;
944 }
945