[src/root.c] Fixed root(x,k) for large very values of k: MPFR could

crash or take much memory, and the function could return incorrect
values. This is done by changing the algorithm for k > 100. This patch
is rather safe as applications that only use small values of k (≤ 100)
won't see any change, and for applications that may use large values
of k, this change may be needed.
[tests/troot.c] Added tests.
(merged changesets r9921-9928,9930-9931,9965-9970,9973-9974,10036,10038
from the trunk, with mpfr_flags_t replaced by unsigned int; for r9965,
the merged was done only on the src directory)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/3.1@10040 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 194 insertions, 39 deletions

diff --git a/tests/troot.c b/tests/troot.c
index 4887e00c1..ae1af0026 100644
--- a/tests/troot.c
+++ b/tests/troot.c
@@ -25,6 +25,19 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 #include "mpfr-test.h"
 
+#define DEFN(N)                                                         \
+  static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)        \
+  { return mpfr_root (y, x, N, rnd); }                                  \
+  static int pow##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)         \
+  { return mpfr_pow_ui (y, x, N, rnd); }
+
+DEFN(2)
+DEFN(3)
+DEFN(4)
+DEFN(5)
+DEFN(17)
+DEFN(120)
+
 static void
 special (void)
 {
@@ -52,7 +65,7 @@ special (void)
       exit (1);
     }
 
-  /* root(-Inf, 17) =  -Inf */
+  /* root(-Inf, 17) = -Inf */
   mpfr_set_inf (x, -1);
   mpfr_root (y, x, 17, MPFR_RNDN);
   if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0)
@@ -69,7 +82,7 @@ special (void)
       exit (1);
     }
 
-  /* root(+/-0) =  +/-0 */
+  /* root(+/-0, k) = +/-0 for k > 0 */
   mpfr_set_ui (x, 0, MPFR_RNDN);
   mpfr_root (y, x, 17, MPFR_RNDN);
   if (mpfr_cmp_ui (y, 0) || mpfr_sgn (y) < 0)
@@ -190,64 +203,39 @@ special (void)
   i = mpfr_root (y, x, 1, MPFR_RNDN);
   if (mpfr_cmp_ui (x, 17) || i != 0)
     {
-      printf ("Error in root (17^(1/1))\n");
+      printf ("Error in root for 17^(1/1)\n");
       exit (1);
     }
 
-#if 0
-  /* Check for k == 0:
-     For 0 <= x < 1 => +0.
-     For x = 1      => 1.
-     For x > 1,     => +Inf.
-     For x < 0      => NaN.   */
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_INF (y) || !MPFR_IS_POS (y) || i != 0)
-    {
-      printf ("Error in root 17^(1/0)\n");
-      exit (1);
-    }
-  mpfr_set_ui (x, 1, MPFR_RNDN);
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (mpfr_cmp_ui (y, 1) || i != 0)
-    {
-      printf ("Error in root 1^(1/0)\n");
-      exit (1);
-    }
   mpfr_set_ui (x, 0, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
-    {
-      printf ("Error in root 0+^(1/0)\n");
-      exit (1);
-    }
-  MPFR_CHANGE_SIGN (x);
-  i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
+  if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root 0-^(1/0)\n");
+      printf ("Error in root for (+0)^(1/0)\n");
       exit (1);
     }
-  mpfr_set_ui_2exp (x, 17, -5, MPFR_RNDD);
+  mpfr_neg (x, x, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
-  if (!MPFR_IS_ZERO (y) || !MPFR_IS_POS (y) || i != 0)
+  if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root (17/2^5)^(1/0)\n");
+      printf ("Error in root for (-0)^(1/0)\n");
       exit (1);
     }
-#endif
-  mpfr_set_ui (x, 0, MPFR_RNDN);
+
+  mpfr_set_ui (x, 1, MPFR_RNDN);
   i = mpfr_root (y, x, 0, MPFR_RNDN);
   if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root 0+^(1/0)\n");
+      printf ("Error in root for 1^(1/0)\n");
       exit (1);
     }
+
   /* Check for k==2 */
   mpfr_set_si (x, -17, MPFR_RNDD);
   i = mpfr_root (y, x, 2, MPFR_RNDN);
   if (!MPFR_IS_NAN (y) || i != 0)
     {
-      printf ("Error in root (-17)^(1/2)\n");
+      printf ("Error in root for (-17)^(1/2)\n");
       exit (1);
     }
 
@@ -255,11 +243,168 @@ special (void)
   mpfr_clear (y);
 }
 
+/* https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=812779
+ * https://bugzilla.gnome.org/show_bug.cgi?id=756960
+ * is a GNOME Calculator bug (mpfr_root applied on a negative integer,
+ * which is converted to an unsigned integer), but the strange result
+ * is also due to a bug in MPFR.
+ */
+static void
+bigint (void)
+{
+  mpfr_t x, y;
+
+  mpfr_inits2 (64, x, y, (mpfr_ptr) 0);
+
+  mpfr_set_ui (x, 10, MPFR_RNDN);
+  if (sizeof (unsigned long) * CHAR_BIT == 64)
+    {
+      mpfr_root (x, x, ULONG_MAX, MPFR_RNDN);
+      mpfr_set_ui_2exp (y, 1, -63, MPFR_RNDN);
+      mpfr_add_ui (y, y, 1, MPFR_RNDN);
+      if (! mpfr_equal_p (x, y))
+        {
+          printf ("Error in bigint for ULONG_MAX\n");
+          printf ("Expected ");
+          mpfr_dump (y);
+          printf ("Got      ");
+          mpfr_dump (x);
+          exit (1);
+        }
+    }
+
+  mpfr_set_ui (x, 10, MPFR_RNDN);
+  mpfr_root (x, x, 1234567890, MPFR_RNDN);
+  mpfr_set_str_binary (y,
+    "1.00000000000000000000000000001000000000101011000101000110010001");
+  if (! mpfr_equal_p (x, y))
+    {
+      printf ("Error in bigint for 1234567890\n");
+      printf ("Expected ");
+      mpfr_dump (y);
+      printf ("Got      ");
+      mpfr_dump (x);
+      exit (1);
+    }
+
+  mpfr_clears (x, y, (mpfr_ptr) 0);
+}
+
 #define TEST_FUNCTION mpfr_root
 #define INTEGER_TYPE unsigned long
-#define INT_RAND_FUNCTION() (INTEGER_TYPE) (randlimb () % 3 +2)
+#define INT_RAND_FUNCTION() \
+  (INTEGER_TYPE) (randlimb () & 1 ? randlimb () : randlimb () % 3 + 2)
 #include "tgeneric_ui.c"
 
+static void
+exact_powers (unsigned long bmax, unsigned long kmax)
+{
+  long b, k;
+  mpz_t z;
+  mpfr_t x, y;
+  int inex, neg;
+
+  mpz_init (z);
+  for (b = 2; b <= bmax; b++)
+    for (k = 1; k <= kmax; k++)
+      {
+        mpz_ui_pow_ui (z, b, k);
+        mpfr_init2 (x, mpz_sizeinbase (z, 2));
+        mpfr_set_ui (x, b, MPFR_RNDN);
+        mpfr_pow_ui (x, x, k, MPFR_RNDN);
+        mpz_set_ui (z, b);
+        mpfr_init2 (y, mpz_sizeinbase (z, 2));
+        for (neg = 0; neg <= 1; neg++)
+          {
+            inex = mpfr_root (y, x, k, MPFR_RNDN);
+            if (inex != 0)
+              {
+                printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                printf ("Expected inex=0, got %d\n", inex);
+                exit (1);
+              }
+            if (neg && (k & 1) == 0)
+              {
+                if (!MPFR_IS_NAN (y))
+                  {
+                    printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                    printf ("Expected y=NaN\n");
+                    printf ("Got      ");
+                    mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
+                    printf ("\n");
+                    exit (1);
+                  }
+              }
+            else if (MPFR_IS_NAN (y) || mpfr_cmp_si (y, b) != 0)
+              {
+                printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
+                printf ("Expected y=%ld\n", b);
+                printf ("Got      ");
+                mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
+                printf ("\n");
+                exit (1);
+              }
+            mpfr_neg (x, x, MPFR_RNDN);
+            b = -b;
+          }
+        mpfr_clear (x);
+        mpfr_clear (y);
+      }
+  mpz_clear (z);
+}
+
+/* Compare root(x,2^h) with pow(x,2^(-h)). */
+static void
+cmp_pow (void)
+{
+  mpfr_t x, y1, y2;
+  int h;
+
+  mpfr_inits2 (128, x, y1, y2, (mpfr_ptr) 0);
+
+  for (h = 1; h < sizeof (unsigned long) * CHAR_BIT; h++)
+    {
+      unsigned long k = (unsigned long) 1 << h;
+      int i;
+
+      for (i = 0; i < 10; i++)
+        {
+          mpfr_rnd_t rnd;
+          unsigned int flags1, flags2;
+          int inex1, inex2;
+
+          tests_default_random (x, 0, __gmpfr_emin, __gmpfr_emax, 1);
+          rnd = RND_RAND ();
+          mpfr_set_ui_2exp (y1, 1, -h, MPFR_RNDN);
+          mpfr_clear_flags ();
+          inex1 = mpfr_pow (y1, x, y1, rnd);
+          flags1 = __gmpfr_flags;
+          mpfr_clear_flags ();
+          inex2 = mpfr_root (y2, x, k, rnd);
+          flags2 = __gmpfr_flags;
+          if (!(mpfr_equal_p (y1, y2) && SAME_SIGN (inex1, inex2) &&
+                flags1 == flags2))
+            {
+              printf ("Error in cmp_pow on h=%d, i=%d, rnd=%s\n",
+                      h, i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
+              printf ("x = ");
+              mpfr_dump (x);
+              printf ("pow  = ");
+              mpfr_dump (y1);
+              printf ("with inex = %d, flags =", inex1);
+              flags_out (flags1);
+              printf ("root = ");
+              mpfr_dump (y2);
+              printf ("with inex = %d, flags =", inex2);
+              flags_out (flags2);
+              exit (1);
+            }
+        }
+    }
+
+  mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
+}
+
 int
 main (void)
 {
@@ -270,7 +415,10 @@ main (void)
 
   tests_start_mpfr ();
 
+  exact_powers (3, 1000);
   special ();
+  bigint ();
+  cmp_pow ();
 
   mpfr_init (x);
 
@@ -329,6 +477,13 @@ main (void)
 
   test_generic_ui (2, 200, 30);
 
+  bad_cases (root2, pow2, "mpfr_root[2]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root3, pow3, "mpfr_root[3]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root4, pow4, "mpfr_root[4]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root5, pow5, "mpfr_root[5]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root17, pow17, "mpfr_root[17]", 8, -256, 255, 4, 128, 800, 40);
+  bad_cases (root120, pow120, "mpfr_root[120]", 8, -256, 255, 4, 128, 800, 40);
+
   tests_end_mpfr ();
   return 0;
 }