[src/urandom.c] Generate the exception flags for mpfr_urandom almost as

  expected:
    * Set the underflow flag if the drawn exponent is less than emin.
      This corresponds to "underflow before rounding" while the normal
      rule in MPFR is "underflow after rounding". This is not fixable
      in the 3.1 branch since the significand is not drawn in this
      case, and drawing the significand to determine whether there is
      an underflow would change the state of the PRNG at the end of
      the function, breaking the ABI and the users' expectations.
    * Set the inexact flag on underflow (this was not done before).
[tests/turandom.c] Added underflow and overflow tests.
[doc/mpfr.texi] Documented the above issue with the underflow flag for
  mpfr_urandom.
(merged changesets r11220,11635,11637-11641,11643-11647 from the trunk;
instead of merging the conflicting r11636, did all the changes manually;
replaced mpfr_flags_t by unsigned int)

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

3 files changed, 220 insertions, 28 deletions

diff --git a/doc/mpfr.texi b/doc/mpfr.texi
index 8de7cad68..167fb6f13 100644
--- a/doc/mpfr.texi
+++ b/doc/mpfr.texi
@@ -2740,6 +2740,15 @@ created using the GMP @code{gmp_randinit} function (see the GMP manual).
 
 Note: the note for @code{mpfr_urandomb} holds too. In addition, the exponent
 range and the rounding mode might have a side effect on the next random state.
+
+The rule for the underflow flag is here ``underflow before rounding''
+instead of the usual ``underflow after rounding''. The reason is that
+the exponent is drawn first, and if it is smaller than the minimum
+exponent, the significand is not drawn. To fix the behavior on the
+underflow flag, one would have to draw the significand in some cases,
+meaning that the behavior of the random generator would change, thus
+it would break the ABI for the MPFR 3.1 branch. However, the observed
+behavior always corresponds to an existing number.
 @end deftypefun
 
 @deftypefun int mpfr_grandom (mpfr_t @var{rop1}, mpfr_t @var{rop2}, gmp_randstate_t @var{state}, mpfr_rnd_t @var{rnd})
diff --git a/src/urandom.c b/src/urandom.c
index 0145be52c..0283bdccb 100644
--- a/src/urandom.c
+++ b/src/urandom.c
@@ -37,6 +37,25 @@ random_rounding_bit (gmp_randstate_t rstate)
   return r & MPFR_LIMB_ONE;
 }
 
+/* NOTE: The current behavior is to consider "underflow before rounding"
+   (the significand does not need to be drawn), while the rule in MPFR
+   is "underflow after rounding". This is unfixable in this 3.1 branch
+   without changing the behavior of the PRNG (thus breaking the ABI). */
+
+/* The mpfr_urandom() function is implemented in the following way for
+   the generic case.
+   1. One determines the exponent exp: 0 with probability 1/2, -1 with
+      probability 1/4, -2 with probability 1/8, etc.
+   2. One draws a 1-ulp interval ]a,b[ containing the exact result (the
+      interval can be regarded as open since it has the same measure as
+      the closed interval).
+   3. Rounding is done. For the directed rounding modes, the rounded value
+      is uniquely determined. For rounding to nearest, ]a,m[ and ]m,b[,
+      where m = (a+b)/2, have the same measure, so that one gets a or b
+      with equal probabilities.
+   Note: Only low-level functions are used (except just before a "return"),
+   so that we do not need MPFR_SAVE_EXPO_*.
+*/
 
 int
 mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode)
@@ -46,33 +65,39 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode)
   mp_size_t nlimbs;
   mp_size_t n;
   mpfr_exp_t exp;
-  mpfr_exp_t emin;
   int cnt;
   int inex;
 
   rp = MPFR_MANT (rop);
   nbits = MPFR_PREC (rop);
-  nlimbs = MPFR_LIMB_SIZE (rop);
   MPFR_SET_POS (rop);
-  exp = 0;
-  emin = mpfr_get_emin ();
-  if (MPFR_UNLIKELY (emin > 0))
+
+  if (MPFR_UNLIKELY (__gmpfr_emin > 0))
     {
+      /* The minimum positive representable number 2^(emin-1) is >= 1,
+         so that we need to round to +0 or 2^(emin-1). For the directed
+         rounding modes, the rounded value is uniquely determined. For
+         rounding to nearest: if emin = 1, one has probability 1/2 for
+         each; otherwise (i.e. if emin > 1), the rounded value is 0. */
+      __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW;
       if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
-          || (emin == 1 && rnd_mode == MPFR_RNDN
+          || (__gmpfr_emin == 1 && rnd_mode == MPFR_RNDN
               && random_rounding_bit (rstate)))
         {
-          mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode);
-          return +1;
+          mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode);
+          MPFR_RET (+1);
         }
       else
         {
           MPFR_SET_ZERO (rop);
-          return -1;
+          MPFR_RET (-1);
         }
     }
 
-  /* Exponent */
+  exp = 0;
+  MPFR_ASSERTD (exp >= __gmpfr_emin);
+
+  /* Step 1 (exponent). */
 #define DRAW_BITS 8 /* we draw DRAW_BITS at a time */
   cnt = DRAW_BITS;
   MPFR_ASSERTN(DRAW_BITS <= GMP_NUMB_BITS);
@@ -89,7 +114,7 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode)
         }
       exp -= cnt;  /* no integer overflow */
 
-      if (MPFR_UNLIKELY (exp < emin))
+      if (MPFR_UNLIKELY (exp < __gmpfr_emin))
         {
           /* To get here, we have been drawing more than -emin zeros
              in a row, then return 0 or the smallest representable
@@ -99,46 +124,46 @@ mpfr_urandom (mpfr_ptr rop, gmp_randstate_t rstate, mpfr_rnd_t rnd_mode)
              If exp == emin - 1, the rounding bit is set, except
              if cnt == DRAW_BITS in which case the rounding bit is
              outside rp[0] and must be generated. */
+          __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW;
           if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
-              || (rnd_mode == MPFR_RNDN && exp == emin - 1
+              || (rnd_mode == MPFR_RNDN && exp == __gmpfr_emin - 1
                   && (cnt != DRAW_BITS || random_rounding_bit (rstate))))
             {
-              mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode);
-              return +1;
+              mpfr_set_ui_2exp (rop, 1, __gmpfr_emin - 1, rnd_mode);
+              MPFR_RET (+1);
             }
           else
             {
               MPFR_SET_ZERO (rop);
-              return -1;
+              MPFR_RET (-1);
             }
         }
+      MPFR_ASSERTD (exp >= __gmpfr_emin);
     }
-  MPFR_EXP (rop) = exp; /* Warning: may be outside the current
-                           exponent range */
 
+  MPFR_ASSERTD (exp >= __gmpfr_emin);
+  MPFR_EXP (rop) = exp; /* Warning: may be larger than emax */
 
-  /* Significand: we need generate only nbits-1 bits, since the most
-     significant is 1 */
+  /* Step 2 (significand): we need generate only nbits-1 bits, since the
+     most significant bit is 1. */
   mpfr_rand_raw (rp, rstate, nbits - 1);
+  nlimbs = MPFR_LIMB_SIZE (rop);
   n = nlimbs * GMP_NUMB_BITS - nbits;
   if (MPFR_LIKELY (n != 0)) /* this will put the low bits to zero */
     mpn_lshift (rp, rp, nlimbs, n);
-
-  /* Set the msb to 1 since it was fixed by the exponent choice */
   rp[nlimbs - 1] |= MPFR_LIMB_HIGHBIT;
 
   /* Rounding */
   if (rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDA
       || (rnd_mode == MPFR_RNDN && random_rounding_bit (rstate)))
     {
-      /* Take care of the exponent range: it may have been reduced */
-      if (exp < emin)
-        mpfr_set_ui_2exp (rop, 1, emin - 1, rnd_mode);
-      else if (exp > mpfr_get_emax ())
-        mpfr_set_inf (rop, +1); /* overflow, flag set by mpfr_check_range */
+      if (MPFR_UNLIKELY (exp > __gmpfr_emax))
+        mpfr_set_inf (rop, +1);  /* overflow */
       else
         mpfr_nextabove (rop);
       inex = +1;
+      /* There is an overflow in the first case and possibly in the second
+         case. If this occurs, the flag will be set by mpfr_check_range. */
     }
   else
     inex = -1;
diff --git a/tests/turandom.c b/tests/turandom.c
index ec5ff8ff9..bedd98ca4 100644
--- a/tests/turandom.c
+++ b/tests/turandom.c
@@ -38,6 +38,7 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
   long count = 0;
   int i;
   int inex = 1;
+  unsigned int ex_flags, flags;
 
   size_tab = (nbtests >= 1000 ? nbtests / 50 : 20);
   tab = (int *) calloc (size_tab, sizeof(int));
@@ -60,7 +61,10 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
 
   for (k = 0; k < nbtests; k++)
     {
+      mpfr_clear_flags ();
+      ex_flags = MPFR_FLAGS_INEXACT;
       i = mpfr_urandom (x, RANDS, rnd);
+      flags = __gmpfr_flags;
       inex = (i != 0) && inex;
       /* check that lower bits are zero */
       if (MPFR_MANT(x)[0] & MPFR_LIMB_MASK(sh) && !MPFR_IS_ZERO (x))
@@ -76,6 +80,17 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
           mpfr_print_binary (x); puts ("");
           exit (1);
         }
+      /* check the flags (an underflow is theoretically possible, but
+         impossible in practice due to the huge exponent range) */
+      if (flags != ex_flags)
+        {
+          printf ("Error: mpfr_urandom() returns incorrect flags.\n");
+          printf ("Expected ");
+          flags_out (ex_flags);
+          printf ("Got      ");
+          flags_out (flags);
+          exit (1);
+        }
 
       d = mpfr_get_d1 (x); av += d; var += d*d;
       i = (int)(size_tab * d);
@@ -98,7 +113,20 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
   for (k = 0; k < 5; k++)
     {
       set_emin (k+1);
+      mpfr_clear_flags ();
+      ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
       inex = mpfr_urandom (x, RANDS, rnd);
+      flags = __gmpfr_flags;
+      if (flags != ex_flags)
+        {
+          printf ("Error: mpfr_urandom() returns incorrect flags"
+                  " for emin = %d.\n", k+1);
+          printf ("Expected ");
+          flags_out (ex_flags);
+          printf ("Got      ");
+          flags_out (flags);
+          exit (1);
+        }
       if ((   (rnd == MPFR_RNDZ || rnd == MPFR_RNDD)
               && (!MPFR_IS_ZERO (x) || inex != -1))
           || ((rnd == MPFR_RNDU || rnd == MPFR_RNDA)
@@ -130,7 +158,7 @@ test_urandom (long nbtests, mpfr_prec_t prec, mpfr_rnd_t rnd, long bit_index,
   printf ("Average = %.5f\nVariance = %.5f\n", av, var);
   printf ("Repartition for urandom with rounding mode %s. "
           "Each integer should be close to %d.\n",
-         mpfr_print_rnd_mode (rnd), (int)th);
+          mpfr_print_rnd_mode (rnd), (int) th);
 
   for (k = 0; k < size_tab; k++)
     {
@@ -209,6 +237,134 @@ bug20170123 (void)
   mpfr_set_emin (emin);
 }
 
+static void
+underflow_tests (void)
+{
+  mpfr_t x;
+  mpfr_exp_t emin;
+  int i, k;
+  int inex;
+  int rnd;
+  unsigned int ex_flags, flags;
+
+  emin = mpfr_get_emin ();
+  mpfr_init2 (x, 4);
+  ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; /* if underflow */
+  for (i = 2; i >= -4; i--)
+    {
+      mpfr_set_emin (i);
+      RND_LOOP (rnd)
+        for (k = 0; k < 100; k++)
+          {
+            mpfr_clear_flags ();
+            inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd);
+            flags = __gmpfr_flags;
+            MPFR_ASSERTN (mpfr_inexflag_p ());
+            if (MPFR_IS_NEG (x))
+              {
+                printf ("Error in underflow_tests: got a negative sign"
+                        " for i=%d rnd=%s k=%d.\n",
+                        i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                exit (1);
+              }
+            if (MPFR_IS_ZERO (x))
+              {
+                if (rnd == MPFR_RNDU || rnd == MPFR_RNDA)
+                  {
+                    printf ("Error in underflow_tests: the value cannot"
+                            " be 0 for i=%d rnd=%s k=%d.\n",
+                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                    exit (1);
+                  }
+                if (flags != ex_flags)
+                  {
+                    printf ("Error in underflow_tests: incorrect flags"
+                            " for i=%d rnd=%s k=%d.\n",
+                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                    printf ("Expected ");
+                    flags_out (ex_flags);
+                    printf ("Got      ");
+                    flags_out (flags);
+                    exit (1);
+                  }
+              }
+            if (inex == 0 || (MPFR_IS_ZERO (x) && inex > 0))
+              {
+                printf ("Error in underflow_tests: incorrect inex (%d)"
+                        " for i=%d rnd=%s k=%d.\n", inex,
+                        i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                exit (1);
+              }
+          }
+    }
+  mpfr_clear (x);
+  mpfr_set_emin (emin);
+}
+
+static void
+overflow_tests (void)
+{
+  mpfr_t x;
+  mpfr_exp_t emax;
+  int i, k;
+  int inex;
+  int rnd;
+  unsigned int ex_flags, flags;
+
+  emax = mpfr_get_emax ();
+  mpfr_init2 (x, 4);
+  ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; /* if overflow */
+  for (i = -4; i <= 0; i++)
+    {
+      mpfr_set_emax (i);
+      RND_LOOP (rnd)
+        for (k = 0; k < 100; k++)
+          {
+            mpfr_clear_flags ();
+            inex = mpfr_urandom (x, mpfr_rands, (mpfr_rnd_t) rnd);
+            flags = __gmpfr_flags;
+            MPFR_ASSERTN (mpfr_inexflag_p ());
+            if (MPFR_IS_NEG (x))
+              {
+                printf ("Error in overflow_tests: got a negative sign"
+                        " for i=%d rnd=%s k=%d.\n",
+                        i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                exit (1);
+              }
+            if (MPFR_IS_INF (x))
+              {
+                if (rnd == MPFR_RNDD || rnd == MPFR_RNDZ)
+                  {
+                    printf ("Error in overflow_tests: the value cannot"
+                            " be +inf for i=%d rnd=%s k=%d.\n",
+                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                    exit (1);
+                  }
+                if (flags != ex_flags)
+                  {
+                    printf ("Error in overflow_tests: incorrect flags"
+                            " for i=%d rnd=%s k=%d.\n",
+                            i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                    printf ("Expected ");
+                    flags_out (ex_flags);
+                    printf ("Got      ");
+                    flags_out (flags);
+                    exit (1);
+                  }
+              }
+            if (inex == 0 || (MPFR_IS_INF (x) && inex < 0))
+              {
+                printf ("Error in overflow_tests: incorrect inex (%d)"
+                        " for i=%d rnd=%s k=%d.\n", inex,
+                        i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), k);
+                exit (1);
+              }
+          }
+    }
+  mpfr_clear (x);
+  mpfr_set_emax (emax);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -260,8 +416,10 @@ main (int argc, char *argv[])
         }
     }
 
-  bug20100914 ();
+  underflow_tests ();
+  overflow_tests ();
 
+  bug20100914 ();
   bug20170123 ();
 
   tests_end_mpfr ();