master 2ef037c0dd 1/2: Improve random bignum generation

[Paul Eggert]

branch: master
commit 2ef037c0dd3510a51ad73fdead1ded09848166f4
Author: Paul Eggert <eggert@cs.ucla.edu>
Commit: Paul Eggert <eggert@cs.ucla.edu>

    Improve random bignum generation
    
    * src/bignum.c (get_random_limb, get_random_limb_lim)
    (get_random_bignum): New functions, for more-efficient
    generation of random bignums without using Frem etc.
    * src/fns.c (get_random_fixnum): New function.
    (Frandom): Use it, and get_random_bignum.
    Be consistent about signalling nonpositive integer arguments;
    since zero is invalid, Qnatnump is not quite right here.
    * src/sysdep.c (get_random_ulong): New function.
---
 src/bignum.c | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 src/bignum.h |  1 +
 src/fns.c    | 77 ++++++++++++++++++-------------------------------
 src/lisp.h   |  1 +
 src/sysdep.c | 10 +++++++
 5 files changed, 132 insertions(+), 50 deletions(-)

diff --git a/src/bignum.c b/src/bignum.c
index cb5322f291..e4e4d45d68 100644
--- a/src/bignum.c
+++ b/src/bignum.c
@@ -476,3 +476,96 @@ check_int_nonnegative (Lisp_Object x)
   CHECK_INTEGER (x);
   return NILP (Fnatnump (x)) ? 0 : check_integer_range (x, 0, INT_MAX);
 }
+
+/* Return a random mp_limb_t.  */
+
+static mp_limb_t
+get_random_limb (void)
+{
+  if (GMP_NUMB_BITS <= ULONG_WIDTH)
+    return get_random_ulong ();
+
+  /* Work around GCC -Wshift-count-overflow false alarm.  */
+  int shift = GMP_NUMB_BITS <= ULONG_WIDTH ? 0 : ULONG_WIDTH;
+
+  /* This is in case someone builds GMP with unusual definitions for
+     MINI_GMP_LIMB_TYPE or _LONG_LONG_LIMB.  */
+  mp_limb_t r = 0;
+  for (int i = 0; i < GMP_NUMB_BITS; i += ULONG_WIDTH)
+    r = (r << shift) | get_random_ulong ();
+  return r;
+}
+
+/* Return a random mp_limb_t I in the range 0 <= I < LIM.
+   If LIM is zero, simply return a random mp_limb_t.  */
+
+static mp_limb_t
+get_random_limb_lim (mp_limb_t lim)
+{
+  /* Return the remainder of a random mp_limb_t R divided by LIM,
+     except reject the rare case where R is so close to the maximum
+     mp_limb_t that the remainder isn't random.  */
+  mp_limb_t difflim = - lim, diff, remainder;
+  do
+    {
+      mp_limb_t r = get_random_limb ();
+      if (lim == 0)
+       return r;
+      remainder = r % lim;
+      diff = r - remainder;
+    }
+  while (difflim < diff);
+
+  return remainder;
+}
+
+/* Return a random Lisp integer I in the range 0 <= I < LIMIT,
+   where LIMIT is a positive bignum.  */
+
+Lisp_Object
+get_random_bignum (struct Lisp_Bignum const *limit)
+{
+  mpz_t const *lim = bignum_val (limit);
+  mp_size_t nlimbs = mpz_size (*lim);
+  eassume (0 < nlimbs);
+  mp_limb_t *r_limb = mpz_limbs_write (mpz[0], nlimbs);
+  mp_limb_t const *lim_limb = mpz_limbs_read (*lim);
+  mp_limb_t limhi = lim_limb[nlimbs - 1];
+  eassert (limhi);
+  bool edgy;
+
+  do
+    {
+      /* Generate the result one limb at a time, most significant first.
+        Choose the most significant limb RHI randomly from 0..LIMHI,
+        where LIMHI is the LIM's first limb, except choose from
+        0..(LIMHI-1) if there is just one limb.  RHI == LIMHI is an
+        unlucky edge case as later limbs might cause the result to be
+        exceed or equal LIM; if this happens, it causes another
+        iteration in the outer loop.  */
+
+      mp_limb_t rhi = get_random_limb_lim (limhi + (1 < nlimbs));
+      edgy = rhi == limhi;
+      r_limb[nlimbs - 1] = rhi;
+
+      for (mp_size_t i = nlimbs - 1; 0 < i--; )
+       {
+         /* get_random_limb_lim (edgy ? limb_lim[i] + 1 : 0)
+            would be wrong here, as the full mp_limb_t range is
+            needed in later limbs for the edge case to have the
+            proper weighting.  */
+         mp_limb_t ri = get_random_limb ();
+         if (edgy)
+           {
+             if (lim_limb[i] < ri)
+               break;
+             edgy = lim_limb[i] == ri;
+           }
+         r_limb[i] = ri;
+       }
+    }
+  while (edgy);
+
+  mpz_limbs_finish (mpz[0], nlimbs);
+  return make_integer_mpz ();
+}
diff --git a/src/bignum.h b/src/bignum.h
index 5f94ce850c..de9ee17c02 100644
--- a/src/bignum.h
+++ b/src/bignum.h
@@ -51,6 +51,7 @@ extern void emacs_mpz_mul_2exp (mpz_t, mpz_t const, EMACS_INT)
 extern void emacs_mpz_pow_ui (mpz_t, mpz_t const, unsigned long)
   ARG_NONNULL ((1, 2));
 extern double mpz_get_d_rounded (mpz_t const) ATTRIBUTE_CONST;
+extern Lisp_Object get_random_bignum (struct Lisp_Bignum const *);
 
 INLINE_HEADER_BEGIN
 
diff --git a/src/fns.c b/src/fns.c
index e8cf185755..6e89fe3ca5 100644
--- a/src/fns.c
+++ b/src/fns.c
@@ -55,41 +55,24 @@ DEFUN ("identity", Fidentity, Sidentity, 1, 1, 0,
   return argument;
 }
 
+/* Return a random Lisp fixnum I in the range 0 <= I < LIM,
+   where LIM is taken from a positive fixnum.  */
 static Lisp_Object
-get_random_bignum (Lisp_Object limit)
+get_random_fixnum (EMACS_INT lim)
 {
-  /* This is a naive transcription into bignums of the fixnum algorithm.
-     I'd be quite surprised if that's anywhere near the best algorithm
-     for it.  */
-  while (true)
+  /* Return the remainder of a random integer R (in range 0..INTMASK)
+     divided by LIM, except reject the rare case where R is so close
+     to INTMASK that the remainder isn't random.  */
+  EMACS_INT difflim = INTMASK - lim + 1, diff, remainder;
+  do
     {
-      Lisp_Object val = make_fixnum (0);
-      Lisp_Object lim = limit;
-      int bits = 0;
-      int bitsperiteration = FIXNUM_BITS - 1;
-      do
-        {
-          /* Shift by one so it is a valid positive fixnum.  */
-          EMACS_INT rand = get_random () >> 1;
-          Lisp_Object lrand = make_fixnum (rand);
-          bits += bitsperiteration;
-          val = CALLN (Flogior,
-                      Fash (val, make_fixnum (bitsperiteration)),
-                      lrand);
-          lim = Fash (lim, make_fixnum (- bitsperiteration));
-        }
-      while (!EQ (lim, make_fixnum (0)));
-      /* Return the remainder, except reject the rare case where
-        get_random returns a number so close to INTMASK that the
-        remainder isn't random.  */
-      Lisp_Object remainder = Frem (val, limit);
-      if (!NILP (CALLN (Fleq,
-                       CALLN (Fminus, val, remainder),
-                       CALLN (Fminus,
-                              Fash (make_fixnum (1), make_fixnum (bits)),
-                              limit))))
-       return remainder;
+      EMACS_INT r = get_random ();
+      remainder = r % lim;
+      diff = r - remainder;
     }
+  while (difflim < diff);
+
+  return make_fixnum (remainder);
 }
 
 DEFUN ("random", Frandom, Srandom, 0, 1, 0,
@@ -103,32 +86,26 @@ With a string argument, set the seed based on the string's 
contents.
 See Info node `(elisp)Random Numbers' for more details.  */)
   (Lisp_Object limit)
 {
-  EMACS_INT val;
-
   if (EQ (limit, Qt))
     init_random ();
   else if (STRINGP (limit))
     seed_random (SSDATA (limit), SBYTES (limit));
-  if (BIGNUMP (limit))
+  else if (FIXNUMP (limit))
+    {
+      EMACS_INT lim = XFIXNUM (limit);
+      if (lim <= 0)
+        xsignal1 (Qargs_out_of_range, limit);
+      return get_random_fixnum (lim);
+    }
+  else if (BIGNUMP (limit))
     {
-      if (0 > mpz_sgn (*xbignum_val (limit)))
-        xsignal2 (Qwrong_type_argument, Qnatnump, limit);
-      return get_random_bignum (limit);
+      struct Lisp_Bignum *lim = XBIGNUM (limit);
+      if (mpz_sgn (*bignum_val (lim)) <= 0)
+        xsignal1 (Qargs_out_of_range, limit);
+      return get_random_bignum (lim);
     }
 
-  val = get_random ();
-  if (FIXNUMP (limit) && 0 < XFIXNUM (limit))
-    while (true)
-      {
-       /* Return the remainder, except reject the rare case where
-          get_random returns a number so close to INTMASK that the
-          remainder isn't random.  */
-       EMACS_INT remainder = val % XFIXNUM (limit);
-       if (val - remainder <= INTMASK - XFIXNUM (limit) + 1)
-         return make_fixnum (remainder);
-       val = get_random ();
-      }
-  return make_ufixnum (val);
+  return make_ufixnum (get_random ());
 }
 
 /* Random data-structure functions.  */
diff --git a/src/lisp.h b/src/lisp.h
index 8053bbc977..c90f901ebc 100644
--- a/src/lisp.h
+++ b/src/lisp.h
@@ -4926,6 +4926,7 @@ extern void child_setup_tty (int);
 extern void setup_pty (int);
 extern int set_window_size (int, int, int);
 extern EMACS_INT get_random (void);
+extern unsigned long int get_random_ulong (void);
 extern void seed_random (void *, ptrdiff_t);
 extern void init_random (void);
 extern void emacs_backtrace (int);
diff --git a/src/sysdep.c b/src/sysdep.c
index b5b18ee6c0..1632f46d13 100644
--- a/src/sysdep.c
+++ b/src/sysdep.c
@@ -2200,6 +2200,16 @@ get_random (void)
   return val & INTMASK;
 }
 
+/* Return a random unsigned long.  */
+unsigned long int
+get_random_ulong (void)
+{
+  unsigned long int r = 0;
+  for (int i = 0; i < (ULONG_WIDTH + RAND_BITS - 1) / RAND_BITS; i++)
+    r = random () ^ (r << RAND_BITS) ^ (r >> (ULONG_WIDTH - RAND_BITS));
+  return r;
+}
+
 #ifndef HAVE_SNPRINTF
 /* Approximate snprintf as best we can on ancient hosts that lack it.  */
 int