Add PowerPC simulator from Andrew Cagney <cagney@highland.com.au>

1 files changed, 1307 insertions, 0 deletions

diff --git a/sim/ppc/dp-bit.c b/sim/ppc/dp-bit.c
new file mode 100644
index 0000000..3313132
--- /dev/null
+++ b/sim/ppc/dp-bit.c
@@ -0,0 +1,1307 @@
+/* This is a software floating point library which can be used instead of
+   the floating point routines in libgcc1.c for targets without hardware
+   floating point.  */
+
+/* Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file with other programs, and to distribute
+those programs without any restriction coming from the use of this
+file.  (The General Public License restrictions do apply in other
+respects; for example, they cover modification of the file, and
+distribution when not linked into another program.)
+
+This file is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+/* As a special exception, if you link this library with other files,
+   some of which are compiled with GCC, to produce an executable,
+   this library does not by itself cause the resulting executable
+   to be covered by the GNU General Public License.
+   This exception does not however invalidate any other reasons why
+   the executable file might be covered by the GNU General Public License.  */
+
+/* This implements IEEE 754 format arithmetic, but does not provide a
+   mechanism for setting the rounding mode, or for generating or handling
+   exceptions.
+
+   The original code by Steve Chamberlain, hacked by Mark Eichin and Jim
+   Wilson, all of Cygnus Support.  */
+
+/* The intended way to use this file is to make two copies, add `#define FLOAT'
+   to one copy, then compile both copies and add them to libgcc.a.  */
+
+/* The following macros can be defined to change the behaviour of this file:
+   FLOAT: Implement a `float', aka SFmode, fp library.  If this is not
+     defined, then this file implements a `double', aka DFmode, fp library.
+   FLOAT_ONLY: Used with FLOAT, to implement a `float' only library, i.e.
+     don't include float->double conversion which requires the double library.
+     This is useful only for machines which can't support doubles, e.g. some
+     8-bit processors.
+   CMPtype: Specify the type that floating point compares should return.
+     This defaults to SItype, aka int.
+   US_SOFTWARE_GOFAST: This makes all entry points use the same names as the
+     US Software goFast library.  If this is not defined, the entry points use
+     the same names as libgcc1.c.
+   _DEBUG_BITFLOAT: This makes debugging the code a little easier, by adding
+     two integers to the FLO_union_type.  
+   NO_NANS: Disable nan and infinity handling
+   SMALL_MACHINE: Useful when operations on QIs and HIs are faster
+     than on an SI */
+
+#ifndef SFtype
+typedef SFtype __attribute__ ((mode (SF)));
+#endif
+#ifndef DFtype
+typedef DFtype __attribute__ ((mode (DF)));
+#endif
+
+#ifndef HItype
+typedef int HItype __attribute__ ((mode (HI)));
+#endif
+#ifndef SItype
+typedef int SItype __attribute__ ((mode (SI)));
+#endif
+#ifndef DItype
+typedef int DItype __attribute__ ((mode (DI)));
+#endif
+
+/* The type of the result of a fp compare */
+#ifndef CMPtype
+#define CMPtype SItype
+#endif
+
+#ifndef UHItype
+typedef unsigned int UHItype __attribute__ ((mode (HI)));
+#endif
+#ifndef USItype
+typedef unsigned int USItype __attribute__ ((mode (SI)));
+#endif
+#ifndef UDItype
+typedef unsigned int UDItype __attribute__ ((mode (DI)));
+#endif
+
+#define MAX_SI_INT   ((SItype) ((unsigned) (~0)>>1))
+#define MAX_USI_INT  ((USItype) ~0)
+
+
+#ifdef FLOAT_ONLY
+#define NO_DI_MODE
+#endif
+
+#ifdef FLOAT
+#	define NGARDS    7L
+#	define GARDROUND 0x3f
+#	define GARDMASK  0x7f
+#	define GARDMSB   0x40
+#	define EXPBITS 8
+#	define EXPBIAS 127
+#	define FRACBITS 23
+#	define EXPMAX (0xff)
+#	define QUIET_NAN 0x100000L
+#	define FRAC_NBITS 32
+#	define FRACHIGH  0x80000000L
+#	define FRACHIGH2 0xc0000000L
+	typedef USItype fractype;
+	typedef UHItype halffractype;
+	typedef SFtype FLO_type;
+	typedef SItype intfrac;
+
+#else
+#	define PREFIXFPDP dp
+#	define PREFIXSFDF df
+#	define NGARDS 8L
+#	define GARDROUND 0x7f
+#	define GARDMASK  0xff
+#	define GARDMSB   0x80
+#	define EXPBITS 11
+#	define EXPBIAS 1023
+#	define FRACBITS 52
+#	define EXPMAX (0x7ff)
+#	define QUIET_NAN 0x8000000000000LL
+#	define FRAC_NBITS 64
+#	define FRACHIGH  0x8000000000000000LL
+#	define FRACHIGH2 0xc000000000000000LL
+	typedef UDItype fractype;
+	typedef USItype halffractype;
+	typedef DFtype FLO_type;
+	typedef DItype intfrac;
+#endif
+
+#ifdef US_SOFTWARE_GOFAST
+#	ifdef FLOAT
+#		define add 		fpadd
+#		define sub 		fpsub
+#		define multiply 	fpmul
+#		define divide 		fpdiv
+#		define compare 		fpcmp
+#		define si_to_float 	sitofp
+#		define float_to_si 	fptosi
+#		define float_to_usi 	fptoui
+#		define negate 		__negsf2
+#		define sf_to_df		fptodp
+#		define dptofp 		dptofp
+#else
+#		define add 		dpadd
+#		define sub 		dpsub
+#		define multiply 	dpmul
+#		define divide 		dpdiv
+#		define compare 		dpcmp
+#		define si_to_float 	litodp
+#		define float_to_si 	dptoli
+#		define float_to_usi 	dptoul
+#		define negate 		__negdf2
+#		define df_to_sf 	dptofp
+#endif
+#else
+#	ifdef FLOAT
+#		define add 		__addsf3
+#		define sub 		__subsf3
+#		define multiply 	__mulsf3
+#		define divide 		__divsf3
+#		define compare 		__cmpsf2
+#		define _eq_f2 		__eqsf2
+#		define _ne_f2 		__nesf2
+#		define _gt_f2 		__gtsf2
+#		define _ge_f2 		__gesf2
+#		define _lt_f2 		__ltsf2
+#		define _le_f2 		__lesf2
+#		define si_to_float 	__floatsisf
+#		define float_to_si 	__fixsfsi
+#		define float_to_usi 	__fixunssfsi
+#		define negate 		__negsf2
+#		define sf_to_df		__extendsfdf2
+#else
+#		define add 		__adddf3
+#		define sub 		__subdf3
+#		define multiply 	__muldf3
+#		define divide 		__divdf3
+#		define compare 		__cmpdf2
+#		define _eq_f2 		__eqdf2
+#		define _ne_f2 		__nedf2
+#		define _gt_f2 		__gtdf2
+#		define _ge_f2 		__gedf2
+#		define _lt_f2 		__ltdf2
+#		define _le_f2 		__ledf2
+#		define si_to_float 	__floatsidf
+#		define float_to_si 	__fixdfsi
+#		define float_to_usi 	__fixunsdfsi
+#		define negate 		__negdf2
+#		define df_to_sf		__truncdfsf2
+#	endif
+#endif
+
+
+#ifndef INLINE
+#define INLINE __inline__
+#endif
+
+/* Preserve the sticky-bit when shifting fractions to the right.  */
+#define LSHIFT(a) { a = (a & 1) | (a >> 1); }
+
+/* numeric parameters */
+/* F_D_BITOFF is the number of bits offset between the MSB of the mantissa
+   of a float and of a double. Assumes there are only two float types.
+   (double::FRAC_BITS+double::NGARGS-(float::FRAC_BITS-float::NGARDS))
+ */
+#define F_D_BITOFF (52+8-(23+7))
+
+
+#define NORMAL_EXPMIN (-(EXPBIAS)+1)
+#define IMPLICIT_1 (1LL<<(FRACBITS+NGARDS))
+#define IMPLICIT_2 (1LL<<(FRACBITS+1+NGARDS))
+
+/* common types */
+
+typedef enum
+{
+  CLASS_SNAN,
+  CLASS_QNAN,
+  CLASS_ZERO,
+  CLASS_NUMBER,
+  CLASS_INFINITY
+} fp_class_type;
+
+typedef struct
+{
+#ifdef SMALL_MACHINE
+  char class;
+  unsigned char sign;
+  short normal_exp;
+#else
+  fp_class_type class;
+  unsigned int sign;
+  int normal_exp;
+#endif
+
+  union
+    {
+      fractype ll;
+      halffractype l[2];
+    } fraction;
+} fp_number_type;
+
+typedef union
+{
+  FLO_type value;
+#ifdef _DEBUG_BITFLOAT
+  int l[2];
+#endif
+  struct
+    {
+#ifndef FLOAT_BIT_ORDER_MISMATCH
+      unsigned int sign:1 __attribute__ ((packed));
+      unsigned int exp:EXPBITS __attribute__ ((packed));
+      fractype fraction:FRACBITS __attribute__ ((packed));
+#else
+      fractype fraction:FRACBITS __attribute__ ((packed));
+      unsigned int exp:EXPBITS __attribute__ ((packed));
+      unsigned int sign:1 __attribute__ ((packed));
+#endif
+    }
+  bits;
+}
+FLO_union_type;
+
+
+/* end of header */
+
+/* IEEE "special" number predicates */
+
+#ifdef NO_NANS
+
+#define nan() 0
+#define isnan(x) 0
+#define isinf(x) 0
+#else
+
+INLINE
+static fp_number_type *
+nan ()
+{
+  static fp_number_type thenan;
+
+  return &thenan;
+}
+
+INLINE
+static int
+isnan ( fp_number_type *  x)
+{
+  return x->class == CLASS_SNAN || x->class == CLASS_QNAN;
+}
+
+INLINE
+static int
+isinf ( fp_number_type *  x)
+{
+  return x->class == CLASS_INFINITY;
+}
+
+#endif
+
+INLINE
+static int
+iszero ( fp_number_type *  x)
+{
+  return x->class == CLASS_ZERO;
+}
+
+INLINE 
+static void
+flip_sign ( fp_number_type *  x)
+{
+  x->sign = !x->sign;
+}
+
+static FLO_type
+pack_d ( fp_number_type *  src)
+{
+  FLO_union_type dst;
+  fractype fraction = src->fraction.ll;	/* wasn't unsigned before? */
+
+  dst.bits.sign = src->sign;
+
+  if (isnan (src))
+    {
+      dst.bits.exp = EXPMAX;
+      dst.bits.fraction = src->fraction.ll;
+      if (src->class == CLASS_QNAN || 1)
+	{
+	  dst.bits.fraction |= QUIET_NAN;
+	}
+    }
+  else if (isinf (src))
+    {
+      dst.bits.exp = EXPMAX;
+      dst.bits.fraction = 0;
+    }
+  else if (iszero (src))
+    {
+      dst.bits.exp = 0;
+      dst.bits.fraction = 0;
+    }
+  else if (fraction == 0)
+    {
+      dst.value = 0;
+    }
+  else
+    {
+      if (src->normal_exp < NORMAL_EXPMIN)
+	{
+	  /* This number's exponent is too low to fit into the bits
+	     available in the number, so we'll store 0 in the exponent and
+	     shift the fraction to the right to make up for it.  */
+
+	  int shift = NORMAL_EXPMIN - src->normal_exp;
+
+	  dst.bits.exp = 0;
+
+	  if (shift > FRAC_NBITS - NGARDS)
+	    {
+	      /* No point shifting, since it's more that 64 out.  */
+	      fraction = 0;
+	    }
+	  else
+	    {
+	      /* Shift by the value */
+	      fraction >>= shift;
+	    }
+	  fraction >>= NGARDS;
+	  dst.bits.fraction = fraction;
+	}
+      else if (src->normal_exp > EXPBIAS)
+	{
+	  dst.bits.exp = EXPMAX;
+	  dst.bits.fraction = 0;
+	}
+      else
+	{
+	  dst.bits.exp = src->normal_exp + EXPBIAS;
+	  /* IF the gard bits are the all zero, but the first, then we're
+	     half way between two numbers, choose the one which makes the
+	     lsb of the answer 0.  */
+	  if ((fraction & GARDMASK) == GARDMSB)
+	    {
+	      if (fraction & (1 << NGARDS))
+		fraction += GARDROUND + 1;
+	    }
+	  else
+	    {
+	      /* Add a one to the guards to round up */
+	      fraction += GARDROUND;
+	    }
+	  if (fraction >= IMPLICIT_2)
+	    {
+	      fraction >>= 1;
+	      dst.bits.exp += 1;
+	    }
+	  fraction >>= NGARDS;
+	  dst.bits.fraction = fraction;
+	}
+    }
+  return dst.value;
+}
+
+static void
+unpack_d (FLO_union_type * src, fp_number_type * dst)
+{
+  fractype fraction = src->bits.fraction;
+
+  dst->sign = src->bits.sign;
+  if (src->bits.exp == 0)
+    {
+      /* Hmm.  Looks like 0 */
+      if (fraction == 0)
+	{
+	  /* tastes like zero */
+	  dst->class = CLASS_ZERO;
+	}
+      else
+	{
+	  /* Zero exponent with non zero fraction - it's denormalized,
+	     so there isn't a leading implicit one - we'll shift it so
+	     it gets one.  */
+	  dst->normal_exp = src->bits.exp - EXPBIAS + 1;
+	  fraction <<= NGARDS;
+
+	  dst->class = CLASS_NUMBER;
+#if 1
+	  while (fraction < IMPLICIT_1)
+	    {
+	      fraction <<= 1;
+	      dst->normal_exp--;
+	    }
+#endif
+	  dst->fraction.ll = fraction;
+	}
+    }
+  else if (src->bits.exp == EXPMAX)
+    {
+      /* Huge exponent*/
+      if (fraction == 0)
+	{
+	  /* Attached to a zero fraction - means infinity */
+	  dst->class = CLASS_INFINITY;
+	}
+      else
+	{
+	  /* Non zero fraction, means nan */
+	  if (dst->sign)
+	    {
+	      dst->class = CLASS_SNAN;
+	    }
+	  else
+	    {
+	      dst->class = CLASS_QNAN;
+	    }
+	  /* Keep the fraction part as the nan number */
+	  dst->fraction.ll = fraction;
+	}
+    }
+  else
+    {
+      /* Nothing strange about this number */
+      dst->normal_exp = src->bits.exp - EXPBIAS;
+      dst->class = CLASS_NUMBER;
+      dst->fraction.ll = (fraction << NGARDS) | IMPLICIT_1;
+    }
+}
+
+static fp_number_type *
+_fpadd_parts (fp_number_type * a,
+	      fp_number_type * b,
+	      fp_number_type * tmp)
+{
+  intfrac tfraction;
+
+  /* Put commonly used fields in local variables.  */
+  int a_normal_exp;
+  int b_normal_exp;
+  fractype a_fraction;
+  fractype b_fraction;
+
+  if (isnan (a))
+    {
+      return a;
+    }
+  if (isnan (b))
+    {
+      return b;
+    }
+  if (isinf (a))
+    {
+      /* Adding infinities with opposite signs yields a NaN.  */
+      if (isinf (b) && a->sign != b->sign)
+	return nan ();
+      return a;
+    }
+  if (isinf (b))
+    {
+      return b;
+    }
+  if (iszero (b))
+    {
+      return a;
+    }
+  if (iszero (a))
+    {
+      return b;
+    }
+
+  /* Got two numbers. shift the smaller and increment the exponent till
+     they're the same */
+  {
+    int diff;
+
+    a_normal_exp = a->normal_exp;
+    b_normal_exp = b->normal_exp;
+    a_fraction = a->fraction.ll;
+    b_fraction = b->fraction.ll;
+
+    diff = a_normal_exp - b_normal_exp;
+
+    if (diff < 0)
+      diff = -diff;
+    if (diff < FRAC_NBITS)
+      {
+	/* ??? This does shifts one bit at a time.  Optimize.  */
+	while (a_normal_exp > b_normal_exp)
+	  {
+	    b_normal_exp++;
+	    LSHIFT (b_fraction);
+	  }
+	while (b_normal_exp > a_normal_exp)
+	  {
+	    a_normal_exp++;
+	    LSHIFT (a_fraction);
+	  }
+      }
+    else
+      {
+	/* Somethings's up.. choose the biggest */
+	if (a_normal_exp > b_normal_exp)
+	  {
+	    b_normal_exp = a_normal_exp;
+	    b_fraction = 0;
+	  }
+	else
+	  {
+	    a_normal_exp = b_normal_exp;
+	    a_fraction = 0;
+	  }
+      }
+  }
+
+  if (a->sign != b->sign)
+    {
+      if (a->sign)
+	{
+	  tfraction = -a_fraction + b_fraction;
+	}
+      else
+	{
+	  tfraction = a_fraction - b_fraction;
+	}
+      if (tfraction > 0)
+	{
+	  tmp->sign = 0;
+	  tmp->normal_exp = a_normal_exp;
+	  tmp->fraction.ll = tfraction;
+	}
+      else
+	{
+	  tmp->sign = 1;
+	  tmp->normal_exp = a_normal_exp;
+	  tmp->fraction.ll = -tfraction;
+	}
+      /* and renormalize it */
+
+      while (tmp->fraction.ll < IMPLICIT_1 && tmp->fraction.ll)
+	{
+	  tmp->fraction.ll <<= 1;
+	  tmp->normal_exp--;
+	}
+    }
+  else
+    {
+      tmp->sign = a->sign;
+      tmp->normal_exp = a_normal_exp;
+      tmp->fraction.ll = a_fraction + b_fraction;
+    }
+  tmp->class = CLASS_NUMBER;
+  /* Now the fraction is added, we have to shift down to renormalize the
+     number */
+
+  if (tmp->fraction.ll >= IMPLICIT_2)
+    {
+      LSHIFT (tmp->fraction.ll);
+      tmp->normal_exp++;
+    }
+  return tmp;
+
+}
+
+FLO_type
+add (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+  fp_number_type tmp;
+  fp_number_type *res;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  res = _fpadd_parts (&a, &b, &tmp);
+
+  return pack_d (res);
+}
+
+FLO_type
+sub (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+  fp_number_type tmp;
+  fp_number_type *res;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  b.sign ^= 1;
+
+  res = _fpadd_parts (&a, &b, &tmp);
+
+  return pack_d (res);
+}
+
+static fp_number_type *
+_fpmul_parts ( fp_number_type *  a,
+	       fp_number_type *  b,
+	       fp_number_type * tmp)
+{
+  fractype low = 0;
+  fractype high = 0;
+
+  if (isnan (a))
+    {
+      a->sign = a->sign != b->sign;
+      return a;
+    }
+  if (isnan (b))
+    {
+      b->sign = a->sign != b->sign;
+      return b;
+    }
+  if (isinf (a))
+    {
+      if (iszero (b))
+	return nan ();
+      a->sign = a->sign != b->sign;
+      return a;
+    }
+  if (isinf (b))
+    {
+      if (iszero (a))
+	{
+	  return nan ();
+	}
+      b->sign = a->sign != b->sign;
+      return b;
+    }
+  if (iszero (a))
+    {
+      a->sign = a->sign != b->sign;
+      return a;
+    }
+  if (iszero (b))
+    {
+      b->sign = a->sign != b->sign;
+      return b;
+    }
+
+  /* Calculate the mantissa by multiplying both 64bit numbers to get a
+     128 bit number */
+  {
+    fractype x = a->fraction.ll;
+    fractype ylow = b->fraction.ll;
+    fractype yhigh = 0;
+    int bit;
+
+#if defined(NO_DI_MODE)
+    {
+      /* ??? This does multiplies one bit at a time.  Optimize.  */
+      for (bit = 0; bit < FRAC_NBITS; bit++)
+	{
+	  int carry;
+
+	  if (x & 1)
+	    {
+	      carry = (low += ylow) < ylow;
+	      high += yhigh + carry;
+	    }
+	  yhigh <<= 1;
+	  if (ylow & FRACHIGH)
+	    {
+	      yhigh |= 1;
+	    }
+	  ylow <<= 1;
+	  x >>= 1;
+	}
+    }
+#elif defined(FLOAT) 
+    {
+      /* Multiplying two 32 bit numbers to get a 64 bit number  on 
+        a machine with DI, so we're safe */
+
+      DItype answer = (DItype)(a->fraction.ll) * (DItype)(b->fraction.ll);
+      
+      high = answer >> 32;
+      low = answer;
+    }
+#else
+    /* Doing a 64*64 to 128 */
+    {
+      UDItype nl = a->fraction.ll & 0xffffffff;
+      UDItype nh = a->fraction.ll >> 32;
+      UDItype ml = b->fraction.ll & 0xffffffff;
+      UDItype mh = b->fraction.ll >>32;
+      UDItype pp_ll = ml * nl;
+      UDItype pp_hl = mh * nl;
+      UDItype pp_lh = ml * nh;
+      UDItype pp_hh = mh * nh;
+      UDItype res2 = 0;
+      UDItype res0 = 0;
+      UDItype ps_hh__ = pp_hl + pp_lh;
+      if (ps_hh__ < pp_hl)
+	res2 += 0x100000000LL;
+      pp_hl = (ps_hh__ << 32) & 0xffffffff00000000LL;
+      res0 = pp_ll + pp_hl;
+      if (res0 < pp_ll)
+	res2++;
+      res2 += ((ps_hh__ >> 32) & 0xffffffffL) + pp_hh;
+      high = res2;
+      low = res0;
+    }
+#endif
+  }
+
+  tmp->normal_exp = a->normal_exp + b->normal_exp;
+  tmp->sign = a->sign != b->sign;
+#ifdef FLOAT
+  tmp->normal_exp += 2;		/* ??????????????? */
+#else
+  tmp->normal_exp += 4;		/* ??????????????? */
+#endif
+  while (high >= IMPLICIT_2)
+    {
+      tmp->normal_exp++;
+      if (high & 1)
+	{
+	  low >>= 1;
+	  low |= FRACHIGH;
+	}
+      high >>= 1;
+    }
+  while (high < IMPLICIT_1)
+    {
+      tmp->normal_exp--;
+
+      high <<= 1;
+      if (low & FRACHIGH)
+	high |= 1;
+      low <<= 1;
+    }
+  /* rounding is tricky. if we only round if it won't make us round later. */
+#if 0
+  if (low & FRACHIGH2)
+    {
+      if (((high & GARDMASK) != GARDMSB)
+	  && (((high + 1) & GARDMASK) == GARDMSB))
+	{
+	  /* don't round, it gets done again later. */
+	}
+      else
+	{
+	  high++;
+	}
+    }
+#endif
+  if ((high & GARDMASK) == GARDMSB)
+    {
+      if (high & (1 << NGARDS))
+	{
+	  /* half way, so round to even */
+	  high += GARDROUND + 1;
+	}
+      else if (low)
+	{
+	  /* but we really weren't half way */
+	  high += GARDROUND + 1;
+	}
+    }
+  tmp->fraction.ll = high;
+  tmp->class = CLASS_NUMBER;
+  return tmp;
+}
+
+FLO_type
+multiply (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+  fp_number_type tmp;
+  fp_number_type *res;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  res = _fpmul_parts (&a, &b, &tmp);
+
+  return pack_d (res);
+}
+
+static fp_number_type *
+_fpdiv_parts (fp_number_type * a,
+	      fp_number_type * b,
+	      fp_number_type * tmp)
+{
+  fractype low = 0;
+  fractype high = 0;
+  fractype r0, r1, y0, y1, bit;
+  fractype q;
+  fractype numerator;
+  fractype denominator;
+  fractype quotient;
+  fractype remainder;
+
+  if (isnan (a))
+    {
+      return a;
+    }
+  if (isnan (b))
+    {
+      return b;
+    }
+  if (isinf (a) || iszero (a))
+    {
+      if (a->class == b->class)
+	return nan ();
+      return a;
+    }
+  a->sign = a->sign ^ b->sign;
+
+  if (isinf (b))
+    {
+      a->fraction.ll = 0;
+      a->normal_exp = 0;
+      return a;
+    }
+  if (iszero (b))
+    {
+      a->class = CLASS_INFINITY;
+      return b;
+    }
+
+  /* Calculate the mantissa by multiplying both 64bit numbers to get a
+     128 bit number */
+  {
+    int carry;
+    intfrac d0, d1;		/* weren't unsigned before ??? */
+
+    /* quotient =
+       ( numerator / denominator) * 2^(numerator exponent -  denominator exponent)
+     */
+
+    a->normal_exp = a->normal_exp - b->normal_exp;
+    numerator = a->fraction.ll;
+    denominator = b->fraction.ll;
+
+    if (numerator < denominator)
+      {
+	/* Fraction will be less than 1.0 */
+	numerator *= 2;
+	a->normal_exp--;
+      }
+    bit = IMPLICIT_1;
+    quotient = 0;
+    /* ??? Does divide one bit at a time.  Optimize.  */
+    while (bit)
+      {
+	if (numerator >= denominator)
+	  {
+	    quotient |= bit;
+	    numerator -= denominator;
+	  }
+	bit >>= 1;
+	numerator *= 2;
+      }
+
+    if ((quotient & GARDMASK) == GARDMSB)
+      {
+	if (quotient & (1 << NGARDS))
+	  {
+	    /* half way, so round to even */
+	    quotient += GARDROUND + 1;
+	  }
+	else if (numerator)
+	  {
+	    /* but we really weren't half way, more bits exist */
+	    quotient += GARDROUND + 1;
+	  }
+      }
+
+    a->fraction.ll = quotient;
+    return (a);
+  }
+}
+
+FLO_type
+divide (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+  fp_number_type tmp;
+  fp_number_type *res;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  res = _fpdiv_parts (&a, &b, &tmp);
+
+  return pack_d (res);
+}
+
+/* according to the demo, fpcmp returns a comparison with 0... thus
+   a<b -> -1
+   a==b -> 0
+   a>b -> +1
+ */
+
+static int
+_fpcmp_parts (fp_number_type * a, fp_number_type * b)
+{
+#if 0
+  /* either nan -> unordered. Must be checked outside of this routine. */
+  if (isnan (a) && isnan (b))
+    {
+      return 1;			/* still unordered! */
+    }
+#endif
+
+  if (isnan (a) || isnan (b))
+    {
+      return 1;			/* how to indicate unordered compare? */
+    }
+  if (isinf (a) && isinf (b))
+    {
+      /* +inf > -inf, but +inf != +inf */
+      /* b    \a| +inf(0)| -inf(1)
+       ______\+--------+--------
+       +inf(0)| a==b(0)| a<b(-1)
+       -------+--------+--------
+       -inf(1)| a>b(1) | a==b(0)
+       -------+--------+--------
+       So since unordered must be non zero, just line up the columns...
+       */
+      return b->sign - a->sign;
+    }
+  /* but not both... */
+  if (isinf (a))
+    {
+      return a->sign ? -1 : 1;
+    }
+  if (isinf (b))
+    {
+      return b->sign ? 1 : -1;
+    }
+  if (iszero (a) && iszero (b))
+    {
+      return 0;
+    }
+  if (iszero (a))
+    {
+      return b->sign ? 1 : -1;
+    }
+  if (iszero (b))
+    {
+      return a->sign ? -1 : 1;
+    }
+  /* now both are "normal". */
+  if (a->sign != b->sign)
+    {
+      /* opposite signs */
+      return a->sign ? -1 : 1;
+    }
+  /* same sign; exponents? */
+  if (a->normal_exp > b->normal_exp)
+    {
+      return a->sign ? -1 : 1;
+    }
+  if (a->normal_exp < b->normal_exp)
+    {
+      return a->sign ? 1 : -1;
+    }
+  /* same exponents; check size. */
+  if (a->fraction.ll > b->fraction.ll)
+    {
+      return a->sign ? -1 : 1;
+    }
+  if (a->fraction.ll < b->fraction.ll)
+    {
+      return a->sign ? 1 : -1;
+    }
+  /* after all that, they're equal. */
+  return 0;
+}
+
+CMPtype
+compare (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  return _fpcmp_parts (&a, &b);
+}
+
+#ifndef US_SOFTWARE_GOFAST
+
+/* These should be optimized for their specific tasks someday.  */
+
+CMPtype
+_eq_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return 1;			/* false, truth == 0 */
+
+  return _fpcmp_parts (&a, &b) ;
+}
+
+CMPtype
+_ne_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return 1;			/* true, truth != 0 */
+
+  return  _fpcmp_parts (&a, &b) ;
+}
+
+CMPtype
+_gt_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return -1;			/* false, truth > 0 */
+
+  return _fpcmp_parts (&a, &b);
+}
+
+CMPtype
+_ge_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return -1;			/* false, truth >= 0 */
+  return _fpcmp_parts (&a, &b) ;
+}
+
+CMPtype
+_lt_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return 1;			/* false, truth < 0 */
+
+  return _fpcmp_parts (&a, &b);
+}
+
+CMPtype
+_le_f2 (FLO_type arg_a, FLO_type arg_b)
+{
+  fp_number_type a;
+  fp_number_type b;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  unpack_d ((FLO_union_type *) & arg_b, &b);
+
+  if (isnan (&a) || isnan (&b))
+    return 1;			/* false, truth <= 0 */
+
+  return _fpcmp_parts (&a, &b) ;
+}
+
+#endif /* ! US_SOFTWARE_GOFAST */
+
+FLO_type
+si_to_float (SItype arg_a)
+{
+  fp_number_type in;
+
+  in.class = CLASS_NUMBER;
+  in.sign = arg_a < 0;
+  if (!arg_a)
+    {
+      in.class = CLASS_ZERO;
+    }
+  else
+    {
+      in.normal_exp = FRACBITS + NGARDS;
+      if (in.sign) 
+	{
+	  /* Special case for minint, since there is no +ve integer
+	     representation for it */
+	  if (arg_a == 0x80000000)
+	    {
+	      return -2147483648.0;
+	    }
+	  in.fraction.ll = (-arg_a);
+	}
+      else
+	in.fraction.ll = arg_a;
+
+      while (in.fraction.ll < (1LL << (FRACBITS + NGARDS)))
+	{
+	  in.fraction.ll <<= 1;
+	  in.normal_exp -= 1;
+	}
+    }
+  return pack_d (&in);
+}
+
+SItype
+float_to_si (FLO_type arg_a)
+{
+  fp_number_type a;
+  SItype tmp;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  if (iszero (&a))
+    return 0;
+  if (isnan (&a))
+    return 0;
+  /* get reasonable MAX_SI_INT... */
+  if (isinf (&a))
+    return a.sign ? MAX_SI_INT : (-MAX_SI_INT)-1;
+  /* it is a number, but a small one */
+  if (a.normal_exp < 0)
+    return 0;
+  if (a.normal_exp > 30)
+    return a.sign ? (-MAX_SI_INT)-1 : MAX_SI_INT;
+  tmp = a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
+  return a.sign ? (-tmp) : (tmp);
+}
+
+#ifdef US_SOFTWARE_GOFAST
+/* While libgcc2.c defines its own __fixunssfsi and __fixunsdfsi routines,
+   we also define them for GOFAST because the ones in libgcc2.c have the
+   wrong names and I'd rather define these here and keep GOFAST CYG-LOC's
+   out of libgcc2.c.  We can't define these here if not GOFAST because then
+   there'd be duplicate copies.  */
+
+USItype
+float_to_usi (FLO_type arg_a)
+{
+  fp_number_type a;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  if (iszero (&a))
+    return 0;
+  if (isnan (&a))
+    return 0;
+  /* get reasonable MAX_USI_INT... */
+  if (isinf (&a))
+    return a.sign ? MAX_USI_INT : 0;
+  /* it is a negative number */
+  if (a.sign)
+    return 0;
+  /* it is a number, but a small one */
+  if (a.normal_exp < 0)
+    return 0;
+  if (a.normal_exp > 31)
+    return MAX_USI_INT;
+  else if (a.normal_exp > (FRACBITS + NGARDS))
+    return a.fraction.ll << ((FRACBITS + NGARDS) - a.normal_exp);
+  else
+    return a.fraction.ll >> ((FRACBITS + NGARDS) - a.normal_exp);
+}
+#endif
+
+FLO_type
+negate (FLO_type arg_a)
+{
+  fp_number_type a;
+
+  unpack_d ((FLO_union_type *) & arg_a, &a);
+  flip_sign (&a);
+  return pack_d (&a);
+}
+
+#ifdef FLOAT
+
+SFtype
+__make_fp(fp_class_type class,
+	     unsigned int sign,
+	     int exp, 
+	     USItype frac)
+{
+  fp_number_type in;
+
+  in.class = class;
+  in.sign = sign;
+  in.normal_exp = exp;
+  in.fraction.ll = frac;
+  return pack_d (&in);
+}
+
+#ifndef FLOAT_ONLY
+
+/* This enables one to build an fp library that supports float but not double.
+   Otherwise, we would get an undefined reference to __make_dp.
+   This is needed for some 8-bit ports that can't handle well values that
+   are 8-bytes in size, so we just don't support double for them at all.  */
+
+extern DFtype __make_dp (fp_class_type, unsigned int, int, UDItype frac);
+
+DFtype
+sf_to_df (SFtype arg_a)
+{
+  fp_number_type in;
+
+  unpack_d ((FLO_union_type *) & arg_a, &in);
+  return __make_dp (in.class, in.sign, in.normal_exp,
+		    ((UDItype) in.fraction.ll) << F_D_BITOFF);
+}
+
+#endif
+#endif
+
+#ifndef FLOAT
+
+extern SFtype __make_fp (fp_class_type, unsigned int, int, USItype);
+
+DFtype
+__make_dp (fp_class_type class, unsigned int sign, int exp, UDItype frac)
+{
+  fp_number_type in;
+
+  in.class = class;
+  in.sign = sign;
+  in.normal_exp = exp;
+  in.fraction.ll = frac;
+  return pack_d (&in);
+}
+
+SFtype
+df_to_sf (DFtype arg_a)
+{
+  fp_number_type in;
+
+  unpack_d ((FLO_union_type *) & arg_a, &in);
+  return __make_fp (in.class, in.sign, in.normal_exp,
+		    in.fraction.ll >> F_D_BITOFF);
+}
+
+#endif