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