Re: Irix6 long doubles implemented wrong? (27_io/ostream_inserter_arith)

[Alexandre Oliva <aoliva at redhat dot com>]

On Jan 19, 2003, Alexandre Oliva <aoliva@redhat.com> wrote:

> the minimum normal exponent should be set such that all 107
> significant bits are available (and it's 107, not 106, as I
> implemented it) for normals.

It can't possibly be 107.  IRIX must be doing something wrong in
setting LDBL_MANT_DIG to 107, since there are only 53 bits of
precision in each double, even if you count the implicit 1s.  Besides,
the emulation I implemented worked correctly down to the least
significant bits with IRIX's printf, so I'm checking these bits in
too, per rth's approval.

Index: gcc/ChangeLog
from  Alexandre Oliva  <aoliva@redhat.com>
	* fp-bit.h: Define macros for TFmode floating-point constants
	in IBM-extended TFmode types.
	(TMODES): Define if __LDBL_MANT_DIG__ has the newly-supported
	widths.
	* config/fp-bit.c (pack_d, unpack_d): Support IBM-extended
	TFmode type.

Index: gcc/config/fp-bit.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/config/fp-bit.c,v
retrieving revision 1.36.4.1
diff -u -p -r1.36.4.1 fp-bit.c
--- gcc/config/fp-bit.c 26 Jan 2003 09:30:38 -0000 1.36.4.1
+++ gcc/config/fp-bit.c 26 Jan 2003 10:01:20 -0000
@@ -322,9 +322,66 @@ pack_d ( fp_number_type *  src)
   dst.bits.exp = exp;
   dst.bits.sign = sign;
 #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 non-zero 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);
   dst.value_raw |= ((fractype) (exp & ((1 << EXPBITS) - 1))) << FRACBITS;
   dst.value_raw |= ((fractype) (sign & 1)) << (FRACBITS | EXPBITS);
+# endif
 #endif
 
 #if defined(FLOAT_WORD_ORDER_MISMATCH) && !defined(FLOAT)
@@ -381,9 +438,42 @@ unpack_d (FLO_union_type * src, fp_numbe
   exp = src->bits.exp;
   sign = src->bits.sign;
 #else
+# if defined TFLOAT && defined HALFFRACBITS
+ {
+   halffractype high, low;
+   
+   high = src->value_raw >> HALFSHIFT;
+   low = src->value_raw & (((fractype)1 << HALFSHIFT) - 1);
+
+   fraction = high & ((((fractype)1) << HALFFRACBITS) - 1);
+   fraction <<= FRACBITS - HALFFRACBITS;
+   exp = ((int)(high >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
+   sign = ((int)(high >> (((HALFFRACBITS + EXPBITS))))) & 1;
+
+   if (exp != EXPMAX && exp != 0 && low != 0)
+     {
+       int lowexp = ((int)(low >> HALFFRACBITS)) & ((1 << EXPBITS) - 1);
+       int shift;
+       fractype xlow;
+
+       xlow = low & ((((fractype)1) << HALFFRACBITS) - 1);
+       if (lowexp)
+	 xlow |= (((halffractype)1) << HALFFRACBITS);
+       else
+	 lowexp = 1;
+       shift = (FRACBITS - HALFFRACBITS) - (exp - lowexp);
+       if (shift > 0)
+	 xlow <<= shift;
+       else if (shift < 0)
+	 xlow >>= -shift;
+       fraction += xlow;
+     }
+ }
+# else
   fraction = src->value_raw & ((((fractype)1) << FRACBITS) - 1);
   exp = ((int)(src->value_raw >> FRACBITS)) & ((1 << EXPBITS) - 1);
   sign = ((int)(src->value_raw >> (FRACBITS + EXPBITS))) & 1;
+# endif
 #endif
 
   dst->sign = sign;
Index: gcc/config/fp-bit.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/config/fp-bit.h,v
retrieving revision 1.7.18.1
diff -u -p -r1.7.18.1 fp-bit.h
--- gcc/config/fp-bit.h 26 Jan 2003 09:30:38 -0000 1.7.18.1
+++ gcc/config/fp-bit.h 26 Jan 2003 10:01:20 -0000
@@ -87,7 +87,7 @@ Boston, MA 02111-1307, USA.  */
 #endif
 #endif /* ! FINE_GRAINED_LIBRARIES */
 
-#if __LDBL_MANT_DIG__ == 113
+#if __LDBL_MANT_DIG__ == 113 || __LDBL_MANT_DIG__ == 106
 # define TMODES
 #endif
 
@@ -150,6 +150,18 @@ typedef unsigned int UTItype __attribute
 #	define FRACHIGH  ((TItype)0x8 << 124)
 #	define FRACHIGH2 ((TItype)0xc << 124)
 #	define FRACBITS 112
+# endif
+
+# if __LDBL_MANT_DIG__ == 106 /* IBM extended (double+double) */
+#	define EXPBITS 11
+#	define EXPBIAS 1023
+#	define EXPMAX (0x7ff)
+#	define QUIET_NAN ((TItype)0x8 << (48 + 64))
+#	define FRACHIGH  ((TItype)0x8 << 124)
+#	define FRACHIGH2 ((TItype)0xc << 124)
+#	define FRACBITS 105
+#	define HALFFRACBITS 52
+#	define HALFSHIFT 64
 # endif
 
 #	define pack_d __pack_t
-- 
Alexandre Oliva   Enjoy Guarana', see http://www.ic.unicamp.br/~oliva/
Red Hat GCC Developer                 aoliva@{redhat.com, gcc.gnu.org}
CS PhD student at IC-Unicamp        oliva@{lsd.ic.unicamp.br, gnu.org}
Free Software Evangelist                Professional serial bug killer