Fix fp-bit.c's paired double handling

[Richard Sandiford]

Geoff's recent real.c patch meant the gcc and fp-bit.c no longer agreed
on the canonical representation of paired doubles.  This was causing
conversion.c to fail on IRIX.

I agree that Geoff's patch, which causes the high double to be rounded
like any other double, is right for IRIX too.  However, fp-bit is still
following the old behaviour, in which the high part used a truncation of
the full significand.  This patch makes it round the high part instead.

The unpack code also mishandled the case where a nonzero low part has
the opposite sign to a power-of-two high part: it would subtract the low
part fraction from zero to form the significand (so lots of spurious
bits would be set) and it would also use the wrong exponent.

Tested on mips-sgi-irix6.5, fixes the conversion failures.
I also tried some other tests by hand.  OK to install?

Richard


	* config/fp-bit.c (pack_d): When using paired doubles to implement
	a long double, round the high part separately.
	(unpack_d): Fix the case in which the high part is a power of two
	and the low part is a nonzero value of the opposite sign.

Index: config/fp-bit.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/config/fp-bit.c,v
retrieving revision 1.43
diff -c -p -F^\([(a-zA-Z0-9_]\|#define\) -r1.43 fp-bit.c
*** config/fp-bit.c	22 Aug 2003 09:33:25 -0000	1.43
--- config/fp-bit.c	20 Jan 2004 08:14:05 -0000
*************** pack_d ( fp_number_type *  src)
*** 330,387 ****
  #else
  # if defined TFLOAT && defined HALFFRACBITS
   {
!    halffractype high, low;
  
!    high = (fraction >> (FRACBITS - HALFFRACBITS));
!    high &= (((fractype)1) << HALFFRACBITS) - 1;
!    high |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << HALFFRACBITS;
!    high |= ((fractype) (sign & 1)) << (HALFFRACBITS | EXPBITS);
  
!    low = (halffractype)fraction &
!      ((((halffractype)1) << (FRACBITS - HALFFRACBITS)) - 1);
  
     if (exp == EXPMAX || exp == 0 || low == 0)
       low = 0;
     else
       {
!        exp -= HALFFRACBITS + 1;
! 
!        while (exp > 0
! 	      && low < ((halffractype)1 << HALFFRACBITS))
  	 {
  	   low <<= 1;
! 	   exp--;
  	 }
  
!        if (exp <= 0)
  	 {
  	   halffractype roundmsb, round;
  
! 	   exp = -exp + 1;
! 
! 	   roundmsb = (1 << (exp - 1));
  	   round = low & ((roundmsb << 1) - 1);
  
! 	   low >>= exp;
! 	   exp = 0;
  
! 	   if (round > roundmsb || (round == roundmsb && (low & 1)))
  	     {
  	       low++;
! 	       if (low >= ((halffractype)1 << HALFFRACBITS))
! 		 /* We don't shift left, since it has just become the
! 		    smallest normal number, whose implicit 1 bit is
! 		    now indicated by the nonzero exponent.  */
! 		 exp++;
  	     }
  	 }
  
!        low &= ((halffractype)1 << HALFFRACBITS) - 1;
!        low |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << HALFFRACBITS;
!        low |= ((fractype) (sign & 1)) << (HALFFRACBITS | EXPBITS);
       }
! 
!    dst.value_raw = (((fractype) high) << HALFSHIFT) | low;
   }
  # else
    dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
--- 330,405 ----
  #else
  # if defined TFLOAT && defined HALFFRACBITS
   {
!    halffractype high, low, unity;
!    int lowsign, lowexp;
! 
!    unity = (halffractype) 1 << HALFFRACBITS;
  
!    /* Set HIGH to the high double's significand, masking out the implicit 1.
!       Set LOW to the low double's full significand.  */
!    high = (fraction >> (FRACBITS - HALFFRACBITS)) & (unity - 1);
!    low = fraction & (unity * 2 - 1);
! 
!    /* Get the initial sign and exponent of the low double.  */
!    lowexp = exp - HALFFRACBITS - 1;
!    lowsign = sign;
! 
!    /* HIGH should be rounded like a normal double, making |LOW| <=
!       0.5 ULP of HIGH.  Assume round-to-nearest.  */
!    if (exp < EXPMAX)
!      if (low > unity || (low == unity && (high & 1) == 1))
!        {
! 	 /* Round HIGH up and adjust LOW to match.  */
! 	 high++;
! 	 if (high == unity)
! 	   {
! 	     /* May make it infinite, but that's OK.  */
! 	     high = 0;
! 	     exp++;
! 	   }
! 	 low = unity * 2 - low;
! 	 lowsign ^= 1;
!        }
  
!    high |= (halffractype) exp << HALFFRACBITS;
!    high |= (halffractype) sign << (HALFFRACBITS + EXPBITS);
  
     if (exp == EXPMAX || exp == 0 || low == 0)
       low = 0;
     else
       {
!        while (lowexp > 0 && low < unity)
  	 {
  	   low <<= 1;
! 	   lowexp--;
  	 }
  
!        if (lowexp <= 0)
  	 {
  	   halffractype roundmsb, round;
+ 	   int shift;
  
! 	   shift = 1 - lowexp;
! 	   roundmsb = (1 << (shift - 1));
  	   round = low & ((roundmsb << 1) - 1);
  
! 	   low >>= shift;
! 	   lowexp = 0;
  
! 	   if (round > roundmsb || (round == roundmsb && (low & 1) == 1))
  	     {
  	       low++;
! 	       if (low == unity)
! 		 /* LOW rounds up to the smallest normal number.  */
! 		 lowexp++;
  	     }
  	 }
  
!        low &= unity - 1;
!        low |= (halffractype) lowexp << HALFFRACBITS;
!        low |= (halffractype) lowsign << (HALFFRACBITS + EXPBITS);
       }
!    dst.value_raw = ((fractype) high << HALFSHIFT) | low;
   }
  # else
    dst.value_raw = fraction & ((((fractype)1) << FRACBITS) - (fractype)1);
*************** unpack_d (FLO_union_type * src, fp_numbe
*** 475,482 ****
  	 xlow >>= -shift;
         if (sign == lowsign)
  	 fraction += xlow;
!        else
  	 fraction -= xlow;
       }
   }
  # else
--- 493,508 ----
  	 xlow >>= -shift;
         if (sign == lowsign)
  	 fraction += xlow;
!        else if (fraction >= xlow)
  	 fraction -= xlow;
+        else
+ 	 {
+ 	   /* The high part is a power of two but the full number is lower.
+ 	      This code will leave the implicit 1 in FRACTION, but we'd
+ 	      have added that below anyway.  */
+ 	   fraction = (((fractype) 1 << FRACBITS) - xlow) << 1;
+ 	   exp--;
+ 	 }
       }
   }
  # else