[Bug math/32411] THREEp96 seems wrong

[jakub at redhat dot com]

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

--- Comment #3 from Jakub Jelinek <jakub at redhat dot com> ---
That doesn't seem to be random numbers but instead all multiplies of 1e-5Q in
certain range.

I've modified it slightly:
diff --git a/main.c b/main.c
index 1ce856a..375642c 100644
--- a/main.c
+++ b/main.c
@@ -5,6 +5,7 @@
 #include <mpfr.h>

 #include <quadmath.h>
+#include <stdlib.h>

 void serialize_qfloat(__float128 f, char *buffer_atleast100) {
        quadmath_snprintf(buffer_atleast100, 100, "%.80Qe", f);
@@ -60,11 +61,11 @@ int compare_values(__float128 exponent, __float128 earlier,
__float128 later) {
 }

 int main(void) {
-    int equal = 0;
-    int not_equal = 0;
+    long long int equal = 0;
+    long long int not_equal = 0;

-    int earlier_is_better = 0;
-    int later_is_better = 0;
+    long long int earlier_is_better = 0;
+    long long int later_is_better = 0;


     // Set the precision to 16384 bits (~4932 decimal digits)
@@ -77,8 +78,34 @@ int main(void) {
             print_qfloat(expq(NUM, 1));
         print_qfloat(expq(NUM, 0));*/

-    for (int i = /*-114330000*/1; i <= 1135000000/* */; i++) {
-        __float128 in = ((__float128)i)*1e-5Q;
+    for (long long i = 0; i < 10000000000LL; i++) {
+        union U { __float128 f; unsigned __int128 i; } u;
+        u.i = rand () & 0x7fffffff;
+        u.i |= (rand () & (unsigned __int128) 0x7fffffff) << 31;
+        u.i |= (rand () & (unsigned __int128) 0x7fffffff) << 62;
+        unsigned r = rand ();
+        u.i |= (r & (unsigned __int128) 0x7ffff) << 93;
+        r >>= 19;
+        u.i |= (r & (unsigned __int128) 1) << 127;
+        if ((r & 14) == 0)
+          {
+            r = rand () & 0x7fffffff;
+            r = (unsigned) r % 0x400d;
+          }
+        else if ((r & 14) == 2)
+          {
+            r = rand () & 0x7fffffff;
+            r = ((unsigned) r % 50) + 0x400d - 50;
+          }
+        else
+          {
+            r = rand () & 0x7fffffff;
+            r = ((unsigned) r % 26) + 0x400d - 26;
+          }
+        u.i |= (r & (unsigned __int128) 0x7fff) << 112;
+        __float128 in = u.f;
+    //for (int i = /*-114330000*/1; i <= 1135000000/* */; i++) {
+    //    __float128 in = ((__float128)i)*1e-5Q;
         __float128 trunk_val = expq(in, 0);
         __float128 new_val = expq(in, 1);

@@ -92,7 +119,7 @@ int main(void) {
             }
         }
     }
-    __builtin_printf("total: %d not_equal: %d earlier_score: %d later_score:
%d\n", equal + not_equal, not_equal, earlier_is_better, later_is_better);
+    __builtin_printf("total: %lld not_equal: %lld earlier_score: %lld
later_score: %lld\n", equal + not_equal, not_equal, earlier_is_better,
later_is_better);

     // Free the memory allocated for MPFR variables
     mpfr_clear(result);

so that it tests pseudo random mantissas and exponents but with strong
preference of the +-0.0002 to +-16383.9999 numbers so that tval1/tval2 aren't
actually 0 all the time.
With this I got
total: 10000000000 not_equal: 29861 earlier_score: 4606 later_score: 25255
and so ~85% of the cases with differences being better and ~15% worse.

It seems that using 3.p94 rather than 3.p96 results in tval2 sometimes one
bigger and sometimes one smaller.
And the code actually uses both delta and result table entries using the worse
tval2 choice, so it isn't completely wrong result, but just one which can have
worse errors.

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