From 2ef037c0dd3510a51ad73fdead1ded09848166f4 Mon Sep 17 00:00:00 2001 From: Paul Eggert Date: Wed, 16 Mar 2022 17:21:55 -0700 Subject: 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 cb5322f291a..e4e4d45d686 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 5f94ce850cf..de9ee17c027 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 e8cf1857550..6e89fe3ca5f 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 8053bbc9777..c90f901ebca 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 b5b18ee6c0f..1632f46d13e 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 -- cgit v1.2.3