urandomb.c: added comments and cleaned up code.

mpfr.texi: improved description of mpfr_urandomb.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@5698 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 22 insertions, 10 deletions

diff --git a/mpfr.texi b/mpfr.texi
index fd2404626..4a9499bb9 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -2243,9 +2243,13 @@ is set to +0.
 
 @deftypefun int mpfr_urandomb (mpfr_t @var{rop}, gmp_randstate_t @var{state})
 Generate a uniformly distributed random float in the interval
-@math{0 @le{} @var{rop} < 1}.
+@math{0 @le{} @var{rop} < 1}. More precisely, the number can be seen as a
+float with a random non-normalized significand and exponent 0, which is then
+normalized (thus if @var{e} denotes the exponent after normalization, then
+the least @math{-@var{e}} significant bits of the significand are always 0).
 Return 0, unless the exponent is not in the current exponent range, in
-which case @var{rop} is set to NaN and a non-zero value is returned. The
+which case @var{rop} is set to NaN and a non-zero value is returned (this
+should never happen in practice, except in very specific cases). The
 second argument is a @code{gmp_randstate_t} structure which should be
 created using the GMP @code{gmp_randinit} function, see the GMP manual.
 @end deftypefun
diff --git a/urandomb.c b/urandomb.c
index 5666e2e66..1d2aa56fd 100644
--- a/urandomb.c
+++ b/urandomb.c
@@ -37,20 +37,22 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
   mp_exp_t exp;
   int cnt;
 
-  MPFR_CLEAR_FLAGS(rop);
+  MPFR_CLEAR_FLAGS (rop);
 
-  rp = MPFR_MANT(rop);
-  nbits = MPFR_PREC(rop);
-  nlimbs = MPFR_LIMB_SIZE(rop);
+  rp = MPFR_MANT (rop);
+  nbits = MPFR_PREC (rop);
+  nlimbs = MPFR_LIMB_SIZE (rop);
   MPFR_SET_POS (rop);
 
+  /* Uniform non-normalized significand */
   _gmp_rand (rp, rstate, nlimbs * BITS_PER_MP_LIMB);
 
   /* If nbits isn't a multiple of BITS_PER_MP_LIMB, mask the low bits */
   cnt = nlimbs * BITS_PER_MP_LIMB - nbits;
-  if (MPFR_LIKELY(cnt != 0))
+  if (MPFR_LIKELY (cnt != 0))
     rp[0] &= ~MPFR_LIMB_MASK (cnt);
 
+  /* Count the null significant limbs and remaining limbs */
   exp = 0;
   k = 0;
   while (nlimbs != 0 && rp[nlimbs - 1] == 0)
@@ -60,22 +62,28 @@ mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
       exp -= BITS_PER_MP_LIMB;
     }
 
-  if (MPFR_LIKELY(nlimbs != 0)) /* otherwise value is zero */
+  if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
     {
       count_leading_zeros (cnt, rp[nlimbs - 1]);
+      /* Normalization */
       if (mpfr_set_exp (rop, exp - cnt))
         {
+          /* If the exponent is not in the current exponent range, we
+             choose to return a NaN as this is probably a user error.
+             Indeed this can happen only if the exponent range has been
+             reduced to a very small interval and/or the precision is
+             huge (very unlikely). */
           MPFR_SET_NAN (rop);
           __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
           return 1;
         }
       if (cnt != 0)
         mpn_lshift (rp + k, rp, nlimbs, cnt);
-      if (k)
+      if (k != 0)
         MPN_ZERO (rp, k);
     }
   else
-    MPFR_SET_ZERO(rop);
+    MPFR_SET_ZERO (rop);
 
   return 0;
 }