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