Return to random2.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/random2.c, Revision 1.1.1.3

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