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 <ja...@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
