186d7f5d3SJohn Marino /* mpz_rrandomb -- Generate a positive random mpz_t of specified bit size, with
286d7f5d3SJohn Marino long runs of consecutive ones and zeros in the binary representation.
386d7f5d3SJohn Marino Meant for testing of other MP routines.
486d7f5d3SJohn Marino
586d7f5d3SJohn Marino Copyright 2000, 2001, 2002, 2004 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
2586d7f5d3SJohn Marino static void gmp_rrandomb __GMP_PROTO ((mp_ptr, gmp_randstate_t, mp_bitcnt_t));
2686d7f5d3SJohn Marino
2786d7f5d3SJohn Marino void
mpz_rrandomb(mpz_ptr x,gmp_randstate_t rstate,mp_bitcnt_t nbits)2886d7f5d3SJohn Marino mpz_rrandomb (mpz_ptr x, gmp_randstate_t rstate, mp_bitcnt_t nbits)
2986d7f5d3SJohn Marino {
3086d7f5d3SJohn Marino mp_size_t nl;
3186d7f5d3SJohn Marino
3286d7f5d3SJohn Marino nl = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
3386d7f5d3SJohn Marino if (nbits != 0)
3486d7f5d3SJohn Marino {
3586d7f5d3SJohn Marino MPZ_REALLOC (x, nl);
3686d7f5d3SJohn Marino gmp_rrandomb (PTR(x), rstate, nbits);
3786d7f5d3SJohn Marino }
3886d7f5d3SJohn Marino
3986d7f5d3SJohn Marino SIZ(x) = nl;
4086d7f5d3SJohn Marino }
4186d7f5d3SJohn Marino
4286d7f5d3SJohn Marino /* Ask _gmp_rand for 32 bits per call unless that's more than a limb can hold.
4386d7f5d3SJohn Marino Thus, we get the same random number sequence in the common cases.
4486d7f5d3SJohn Marino FIXME: We should always generate the same random number sequence! */
4586d7f5d3SJohn Marino #if GMP_NUMB_BITS < 32
4686d7f5d3SJohn Marino #define BITS_PER_RANDCALL GMP_NUMB_BITS
4786d7f5d3SJohn Marino #else
4886d7f5d3SJohn Marino #define BITS_PER_RANDCALL 32
4986d7f5d3SJohn Marino #endif
5086d7f5d3SJohn Marino
5186d7f5d3SJohn Marino static void
gmp_rrandomb(mp_ptr rp,gmp_randstate_t rstate,mp_bitcnt_t nbits)5286d7f5d3SJohn Marino gmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, mp_bitcnt_t nbits)
5386d7f5d3SJohn Marino {
5486d7f5d3SJohn Marino mp_bitcnt_t bi;
5586d7f5d3SJohn Marino mp_limb_t ranm; /* buffer for random bits */
5686d7f5d3SJohn Marino unsigned cap_chunksize, chunksize;
5786d7f5d3SJohn Marino mp_size_t i;
5886d7f5d3SJohn Marino
5986d7f5d3SJohn Marino /* Set entire result to 111..1 */
6086d7f5d3SJohn Marino i = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS - 1;
6186d7f5d3SJohn Marino rp[i] = GMP_NUMB_MAX >> (GMP_NUMB_BITS - (nbits % GMP_NUMB_BITS)) % GMP_NUMB_BITS;
6286d7f5d3SJohn Marino for (i = i - 1; i >= 0; i--)
6386d7f5d3SJohn Marino rp[i] = GMP_NUMB_MAX;
6486d7f5d3SJohn Marino
6586d7f5d3SJohn Marino _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);
6686d7f5d3SJohn Marino cap_chunksize = nbits / (ranm % 4 + 1);
6786d7f5d3SJohn Marino cap_chunksize += cap_chunksize == 0; /* make it at least 1 */
6886d7f5d3SJohn Marino
6986d7f5d3SJohn Marino bi = nbits;
7086d7f5d3SJohn Marino
7186d7f5d3SJohn Marino for (;;)
7286d7f5d3SJohn Marino {
7386d7f5d3SJohn Marino _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);
7486d7f5d3SJohn Marino chunksize = 1 + ranm % cap_chunksize;
7586d7f5d3SJohn Marino bi = (bi < chunksize) ? 0 : bi - chunksize;
7686d7f5d3SJohn Marino
7786d7f5d3SJohn Marino if (bi == 0)
7886d7f5d3SJohn Marino break; /* low chunk is ...1 */
7986d7f5d3SJohn Marino
8086d7f5d3SJohn Marino rp[bi / GMP_NUMB_BITS] ^= CNST_LIMB (1) << bi % GMP_NUMB_BITS;
8186d7f5d3SJohn Marino
8286d7f5d3SJohn Marino _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);
8386d7f5d3SJohn Marino chunksize = 1 + ranm % cap_chunksize;
8486d7f5d3SJohn Marino bi = (bi < chunksize) ? 0 : bi - chunksize;
8586d7f5d3SJohn Marino
8686d7f5d3SJohn Marino mpn_incr_u (rp + bi / GMP_NUMB_BITS, CNST_LIMB (1) << bi % GMP_NUMB_BITS);
8786d7f5d3SJohn Marino
8886d7f5d3SJohn Marino if (bi == 0)
8986d7f5d3SJohn Marino break; /* low chunk is ...0 */
9086d7f5d3SJohn Marino }
9186d7f5d3SJohn Marino }
92