186d7f5d3SJohn Marino /* mpz_urandomm (rop, state, n) -- Generate a uniform pseudorandom
286d7f5d3SJohn Marino integer in the range 0 to N-1, using STATE as the random state
386d7f5d3SJohn Marino previously initialized by a call to gmp_randinit().
486d7f5d3SJohn Marino
586d7f5d3SJohn Marino Copyright 2000, 2002 Free Software Foundation, Inc.
686d7f5d3SJohn Marino
786d7f5d3SJohn Marino This file is part of the GNU MP Library.
886d7f5d3SJohn Marino
986d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
1086d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
1186d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1286d7f5d3SJohn Marino option) any later version.
1386d7f5d3SJohn Marino
1486d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
1586d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1686d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1786d7f5d3SJohn Marino License for more details.
1886d7f5d3SJohn Marino
1986d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
2086d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
2186d7f5d3SJohn Marino
2286d7f5d3SJohn Marino #include "gmp.h"
2386d7f5d3SJohn Marino #include "gmp-impl.h"
2486d7f5d3SJohn Marino #include "longlong.h" /* for count_leading_zeros */
2586d7f5d3SJohn Marino
2686d7f5d3SJohn Marino
2786d7f5d3SJohn Marino #define MAX_URANDOMM_ITER 80
2886d7f5d3SJohn Marino
2986d7f5d3SJohn Marino void
mpz_urandomm(mpz_ptr rop,gmp_randstate_t rstate,mpz_srcptr n)3086d7f5d3SJohn Marino mpz_urandomm (mpz_ptr rop, gmp_randstate_t rstate, mpz_srcptr n)
3186d7f5d3SJohn Marino {
3286d7f5d3SJohn Marino mp_ptr rp, np, nlast;
3386d7f5d3SJohn Marino mp_size_t nbits, size;
3486d7f5d3SJohn Marino int count;
3586d7f5d3SJohn Marino int pow2;
3686d7f5d3SJohn Marino int cmp;
3786d7f5d3SJohn Marino TMP_DECL;
3886d7f5d3SJohn Marino
3986d7f5d3SJohn Marino size = ABSIZ (n);
4086d7f5d3SJohn Marino if (size == 0)
4186d7f5d3SJohn Marino DIVIDE_BY_ZERO;
4286d7f5d3SJohn Marino
4386d7f5d3SJohn Marino nlast = &PTR (n)[size - 1];
4486d7f5d3SJohn Marino
4586d7f5d3SJohn Marino /* Detect whether n is a power of 2. */
4686d7f5d3SJohn Marino pow2 = POW2_P (*nlast);
4786d7f5d3SJohn Marino if (pow2 != 0)
4886d7f5d3SJohn Marino for (np = PTR (n); np < nlast; np++)
4986d7f5d3SJohn Marino if (*np != 0)
5086d7f5d3SJohn Marino {
5186d7f5d3SJohn Marino pow2 = 0; /* Mark n as `not a power of two'. */
5286d7f5d3SJohn Marino break;
5386d7f5d3SJohn Marino }
5486d7f5d3SJohn Marino
5586d7f5d3SJohn Marino count_leading_zeros (count, *nlast);
5686d7f5d3SJohn Marino nbits = size * GMP_NUMB_BITS - (count - GMP_NAIL_BITS) - pow2;
5786d7f5d3SJohn Marino if (nbits == 0) /* nbits == 0 means that n was == 1. */
5886d7f5d3SJohn Marino {
5986d7f5d3SJohn Marino SIZ (rop) = 0;
6086d7f5d3SJohn Marino return;
6186d7f5d3SJohn Marino }
6286d7f5d3SJohn Marino
6386d7f5d3SJohn Marino TMP_MARK;
6486d7f5d3SJohn Marino np = PTR (n);
6586d7f5d3SJohn Marino if (rop == n)
6686d7f5d3SJohn Marino {
6786d7f5d3SJohn Marino mp_ptr tp;
6886d7f5d3SJohn Marino tp = TMP_ALLOC_LIMBS (size);
6986d7f5d3SJohn Marino MPN_COPY (tp, np, size);
7086d7f5d3SJohn Marino np = tp;
7186d7f5d3SJohn Marino }
7286d7f5d3SJohn Marino
7386d7f5d3SJohn Marino /* Here the allocated size can be one too much if n is a power of
7486d7f5d3SJohn Marino (2^GMP_NUMB_BITS) but it's convenient for using mpn_cmp below. */
7586d7f5d3SJohn Marino rp = MPZ_REALLOC (rop, size);
7686d7f5d3SJohn Marino /* Clear last limb to prevent the case in which size is one too much. */
7786d7f5d3SJohn Marino rp[size - 1] = 0;
7886d7f5d3SJohn Marino
7986d7f5d3SJohn Marino count = MAX_URANDOMM_ITER; /* Set iteration count limit. */
8086d7f5d3SJohn Marino do
8186d7f5d3SJohn Marino {
8286d7f5d3SJohn Marino _gmp_rand (rp, rstate, nbits);
8386d7f5d3SJohn Marino MPN_CMP (cmp, rp, np, size);
8486d7f5d3SJohn Marino }
8586d7f5d3SJohn Marino while (cmp >= 0 && --count != 0);
8686d7f5d3SJohn Marino
8786d7f5d3SJohn Marino if (count == 0)
8886d7f5d3SJohn Marino /* Too many iterations; return result mod n == result - n */
8986d7f5d3SJohn Marino mpn_sub_n (rp, rp, np, size);
9086d7f5d3SJohn Marino
9186d7f5d3SJohn Marino MPN_NORMALIZE (rp, size);
9286d7f5d3SJohn Marino SIZ (rop) = size;
9386d7f5d3SJohn Marino TMP_FREE;
9486d7f5d3SJohn Marino }
95