[Bug math/32411] THREEp96 seems wrong

[jakub at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=32411

--- Comment #4 from Jakub Jelinek <jakub at redhat dot com> ---
Tested glibc with it and no regressions there.
And additionally testing in the comparefloats repo on top of the above patch
with
diff --git a/expq.h b/expq.h
index 3994f40..72ea9c0 100644
--- a/expq.h
+++ b/expq.h
@@ -165,6 +165,24 @@ expq (__float128 x, int use_actual_pow_2_96)
       xl -= __expq_table[T_EXPL_ARG1+2*tval1+1];

       /* Calculate t/32768.  */
+      if (use_actual_pow_2_96 == -1) {
+        __float128 t2 = x + threetimes2pow94;
+        t2 -= threetimes2pow94;
+        int tval3 = (int) (t2 * TWO15);
+        t = x + THREEp96;
+        t -= THREEp96;
+        tval2 = (int) (t * TWO15);
+        if (tval2 == tval2)
+          return 0;
+        __float128 x2 = x, xl2 = xl;
+        x -= __expq_table[T_EXPL_ARG2+2*tval2];
+        xl -= __expq_table[T_EXPL_ARG2+2*tval2+1];
+        x2 -= __expq_table[T_EXPL_ARG2+2*tval3];
+        xl2 -= __expq_table[T_EXPL_ARG2+2*tval3+1];
+        if (__builtin_fabsf128 (x + xl) > __builtin_fabsf128 (x2 + xl2))
+          __builtin_abort ();
+        return 0;
+      }
       if (use_actual_pow_2_96) {
         t = x + threetimes2pow94;
         t -= threetimes2pow94;
@@ -242,6 +260,8 @@ expq (__float128 x, int use_actual_pow_2_96)
          return result;
        }
     }
+  else if (use_actual_pow_2_96 == -1)
+    return 0;
   /* Exceptional cases:  */
   else if (__builtin_isless (x, himark))
     {
diff --git a/main.c b/main.c
index 1ce856a..81b04de 100644
--- a/main.c
+++ b/main.c
@@ -104,6 +104,8 @@ int main(void) {
           }
         u.i |= (r & (unsigned __int128) 0x7fff) << 112;
         __float128 in = u.f;
+        volatile __float128 nv = expq(in, -1);
+        continue;
     //for (int i = /*-114330000*/1; i <= 1135000000/* */; i++) {
     //    __float128 in = ((__float128)i)*1e-5Q;
         __float128 trunk_val = expq(in, 0);
which so far passed 2 billion iterations, so unless I've screwed up something
that means with the fix for the THREEp96 constant x + xl (i.e. the remainder
value for which we use the polynom) is never larger than previously.  Also note
that threetimes2pow94 is misnamed in the repo,
because that one is actually the true threep96.

-- 
You are receiving this mail because:
You are on the CC list for the bug.