Return to rrandomb.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/rrandomb.c, Revision 1.1.1.2

1.1       maekawa     1: /* mpz_rrandomb -- Generate a positive random mpz_t of specified bit size, with
                      2:    long runs of consecutive ones and zeros in the binary representation.
                      3:    Meant for testing of other MP routines.
                      4:
1.1.1.2 ! ohara       5: Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
1.1       maekawa     6:
                      7: This file is part of the GNU MP Library.
                      8:
                      9: The GNU MP Library is free software; you can redistribute it and/or modify
                     10: it under the terms of the GNU Lesser General Public License as published by
                     11: the Free Software Foundation; either version 2.1 of the License, or (at your
                     12: option) any later version.
                     13:
                     14: The GNU MP Library is distributed in the hope that it will be useful, but
                     15: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     16: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
                     17: License for more details.
                     18:
                     19: You should have received a copy of the GNU Lesser General Public License
                     20: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     21: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     22: MA 02111-1307, USA. */
                     23:
                     24: #include "gmp.h"
                     25: #include "gmp-impl.h"
                     26:
                     27: static void gmp_rrandomb _PROTO ((mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits));
                     28:
                     29: void
                     30: mpz_rrandomb (mpz_ptr x, gmp_randstate_t rstate, unsigned long int nbits)
                     31: {
                     32:   mp_size_t nl = 0;
                     33:
                     34:   if (nbits != 0)
                     35:     {
                     36:       mp_ptr xp;
1.1.1.2 ! ohara      37:       nl = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
        !            38:       MPZ_REALLOC (x, nl);
1.1       maekawa    39:
                     40:       xp = PTR(x);
                     41:       gmp_rrandomb (xp, rstate, nbits);
                     42:       MPN_NORMALIZE (xp, nl);
                     43:     }
                     44:
                     45:   SIZ(x) = nl;
                     46: }
                     47:
1.1.1.2 ! ohara      48: /* It's a bit tricky to get this right, so please test the code well if you
        !            49:    hack with it.  Some early versions of the function produced random numbers
        !            50:    with the leading limb == 0, and some versions never made the most
        !            51:    significant bit set.
        !            52:
        !            53:    This code and mpn_random2 are almost identical, though the latter makes bit
        !            54:    runs of 1 to 32, and forces the first block to contain 1-bits.
        !            55:
        !            56:    The random state produces some number of random bits per underlying lc
        !            57:    invocation (BITS_PER_RANDCALL).  We should perhaps ask for that, instead of
        !            58:    asking for 32, presuming that limbs are at least 32 bits.  FIXME: Handle
        !            59:    smaller limbs, such as 4-bit limbs useful for testing purposes, or limbs
        !            60:    truncated by nailing.
        !            61:
        !            62:    For efficiency, we make sure to use most bits returned from _gmp_rand, since
        !            63:    the underlying random number generator is slow.  We keep returned bits in
        !            64:    ranm/ran, and a count of how many bits remaining in ran_nbits.  */
        !            65:
        !            66: #define LOGBITS_PER_BLOCK 4
        !            67:
        !            68: /* Ask _gmp_rand for 32 bits per call unless that's more than a limb can hold.
        !            69:    Thus, we get the same random number sequence in the common cases.
        !            70:    FIXME: We should always generate the same random number sequence!  */
        !            71: #if GMP_NUMB_BITS < 32
        !            72: #define BITS_PER_RANDCALL GMP_NUMB_BITS
        !            73: #else
        !            74: #define BITS_PER_RANDCALL 32
        !            75: #endif
1.1       maekawa    76:
                     77: static void
                     78: gmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
                     79: {
                     80:   int nb;
1.1.1.2 ! ohara      81:   int bit_pos;                 /* bit number of least significant bit where
        !            82:                                   next bit field to be inserted */
        !            83:   mp_size_t ri;                        /* index in rp */
        !            84:   mp_limb_t ran, ranm;         /* buffer for random bits */
        !            85:   mp_limb_t acc;               /* accumulate output random data here */
        !            86:   int ran_nbits;               /* number of valid bits in ran */
        !            87:
        !            88:   ran_nbits = 0;
        !            89:   bit_pos = (nbits - 1) % GMP_NUMB_BITS;
        !            90:   ri = (nbits - 1) / GMP_NUMB_BITS;
1.1       maekawa    91:
                     92:   acc = 0;
1.1.1.2 ! ohara      93:   while (ri >= 0)
1.1       maekawa    94:     {
1.1.1.2 ! ohara      95:       if (ran_nbits < LOGBITS_PER_BLOCK + 1)
        !            96:        {
        !            97:          _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);
        !            98:          ran = ranm;
        !            99:          ran_nbits = BITS_PER_RANDCALL;
        !           100:        }
        !           101:
        !           102:       nb = (ran >> 1) % (1 << LOGBITS_PER_BLOCK) + 1;
1.1       maekawa   103:       if ((ran & 1) != 0)
                    104:        {
1.1.1.2 ! ohara     105:          /* Generate a string of nb ones.  */
1.1       maekawa   106:          if (nb > bit_pos)
                    107:            {
1.1.1.2 ! ohara     108:              rp[ri--] = acc | (((mp_limb_t) 2 << bit_pos) - 1);
        !           109:              bit_pos += GMP_NUMB_BITS;
1.1       maekawa   110:              bit_pos -= nb;
1.1.1.2 ! ohara     111:              acc = ((~(mp_limb_t) 1) << bit_pos) & GMP_NUMB_MASK;
1.1       maekawa   112:            }
                    113:          else
                    114:            {
                    115:              bit_pos -= nb;
1.1.1.2 ! ohara     116:              acc |= (((mp_limb_t) 2 << nb) - 2) << bit_pos;
1.1       maekawa   117:            }
                    118:        }
                    119:       else
                    120:        {
1.1.1.2 ! ohara     121:          /* Generate a string of nb zeroes.  */
1.1       maekawa   122:          if (nb > bit_pos)
                    123:            {
1.1.1.2 ! ohara     124:              rp[ri--] = acc;
1.1       maekawa   125:              acc = 0;
1.1.1.2 ! ohara     126:              bit_pos += GMP_NUMB_BITS;
1.1       maekawa   127:            }
                    128:          bit_pos -= nb;
                    129:        }
1.1.1.2 ! ohara     130:       ran_nbits -= LOGBITS_PER_BLOCK + 1;
        !           131:       ran >>= LOGBITS_PER_BLOCK + 1;
1.1       maekawa   132:     }
                    133: }