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