[PATCH] Fix libquadmath on FLT_EVAL_METHOD != 0 targets

[PATCH] Fix libquadmath on FLT_EVAL_METHOD != 0 targets

[Jakub Jelinek]

Hi!

As
#include <quadmath.h>
#include <stdio.h>

void printq(__float128 x)
{
  char buf[100];
  quadmath_snprintf(buf, 100, "%40.30Qa", x);
  printf("%s\n", buf);
}

int
main()
{
  __float128 twopi = 6.28318530717958647692528676655900576839433879875021164194989Q;
  __float128 three = 3.0Q;
  __float128 two = 2.0Q;
  __float128 ang = two * twopi / three;
  __float128 c = cosq(ang);
  __float128 correctc = -.5;
  __float128 s = sinq(ang);
  __float128 corrects = 0.866025403784438646763723170752936183471402626905190314027903Q;
  printq(twopi);
  printq(three);
  printq(two);
  printq(ang);
  printq(c);
  printq(correctc);
  printq(s);
  printq(corrects);
  return 0;
}
shows, the following two loops (at least the first one matters a lot)
in __quadmath_kernel_rem_pio2 rely on FLT_EVAL_METHOD 0, so sinq
is giving a result within 1ulp or so on x86_64, but a number that is far off
on i686.  If libquadmath is recompiled with -O0, it works even on i686.
Fixed by forced rounding into double precision (don't want to rely on
-fexcess-precision=* setting used to compile the library, so added a
- volatile), committed to trunk/4.6.

2011-08-01  Jakub Jelinek  <jakub@redhat.com>

	* math/rem_pio2q.c (__quadmath_kernel_rem_pio2): Fix up fq to y
	conversion for prec 3 and __FLT_EVAL_METHOD__ != 0.

--- libquadmath/math/rem_pio2q.c.jj	2010-12-14 08:11:24.000000000 +0100
+++ libquadmath/math/rem_pio2q.c	2011-08-01 10:45:27.000000000 +0200
@@ -282,14 +282,20 @@ recompute:
 		break;
 	    case 3:	/* painful */
 		for (i=jz;i>0;i--) {
-		    fw      = fq[i-1]+fq[i];
-		    fq[i]  += fq[i-1]-fw;
-		    fq[i-1] = fw;
+#if __FLT_EVAL_METHOD__ != 0
+		    volatile
+#endif
+		    double fv = (double)(fq[i-1]+fq[i]);
+		    fq[i]  += fq[i-1]-fv;
+		    fq[i-1] = fv;
 		}
 		for (i=jz;i>1;i--) {
-		    fw      = fq[i-1]+fq[i];
-		    fq[i]  += fq[i-1]-fw;
-		    fq[i-1] = fw;
+#if __FLT_EVAL_METHOD__ != 0
+		    volatile
+#endif
+		    double fv = (double)(fq[i-1]+fq[i]);
+		    fq[i]  += fq[i-1]-fv;
+		    fq[i-1] = fv;
 		}
 		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 		if(ih==0) {

	Jakub