===================================================================
RCS file: /home/cvs/OpenXM_contrib/gmp/mpz/Attic/rrandomb.c,v
retrieving revision 1.1.1.1
retrieving revision 1.1.1.2
diff -u -p -r1.1.1.1 -r1.1.1.2
--- OpenXM_contrib/gmp/mpz/Attic/rrandomb.c	2000/09/09 14:12:56	1.1.1.1
+++ OpenXM_contrib/gmp/mpz/Attic/rrandomb.c	2003/08/25 16:06:33	1.1.1.2
@@ -2,7 +2,7 @@
    long runs of consecutive ones and zeros in the binary representation.
    Meant for testing of other MP routines.
 
-Copyright (C) 2000 Free Software Foundation, Inc.
+Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -27,23 +27,15 @@ MA 02111-1307, USA. */
 static void gmp_rrandomb _PROTO ((mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits));
 
 void
-#if __STDC__
 mpz_rrandomb (mpz_ptr x, gmp_randstate_t rstate, unsigned long int nbits)
-#else
-mpz_rrandomb (x, rstate, nbits)
-     mpz_ptr x;
-     gmp_randstate_t rstate;
-     unsigned long int nbits;
-#endif
 {
   mp_size_t nl = 0;
 
   if (nbits != 0)
     {
       mp_ptr xp;
-      nl = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
-      if (x->_mp_alloc < nl)
-	_mpz_realloc (x, nl);
+      nl = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
+      MPZ_REALLOC (x, nl);
 
       xp = PTR(x);
       gmp_rrandomb (xp, rstate, nbits);
@@ -53,65 +45,89 @@ mpz_rrandomb (x, rstate, nbits)
   SIZ(x) = nl;
 }
 
-#define BITS_PER_CHUNK 4
+/* It's a bit tricky to get this right, so please test the code well if you
+   hack with it.  Some early versions of the function produced random numbers
+   with the leading limb == 0, and some versions never made the most
+   significant bit set.
 
-static void
-#if __STDC__
-gmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
+   This code and mpn_random2 are almost identical, though the latter makes bit
+   runs of 1 to 32, and forces the first block to contain 1-bits.
+
+   The random state produces some number of random bits per underlying lc
+   invocation (BITS_PER_RANDCALL).  We should perhaps ask for that, instead of
+   asking for 32, presuming that limbs are at least 32 bits.  FIXME: Handle
+   smaller limbs, such as 4-bit limbs useful for testing purposes, or limbs
+   truncated by nailing.
+
+   For efficiency, we make sure to use most bits returned from _gmp_rand, since
+   the underlying random number generator is slow.  We keep returned bits in
+   ranm/ran, and a count of how many bits remaining in ran_nbits.  */
+
+#define LOGBITS_PER_BLOCK 4
+
+/* Ask _gmp_rand for 32 bits per call unless that's more than a limb can hold.
+   Thus, we get the same random number sequence in the common cases.
+   FIXME: We should always generate the same random number sequence!  */
+#if GMP_NUMB_BITS < 32
+#define BITS_PER_RANDCALL GMP_NUMB_BITS
 #else
-gmp_rrandomb (rp, rstate, nbits)
-     mp_ptr rp;
-     gmp_randstate_t rstate;
-     unsigned long int nbits;
+#define BITS_PER_RANDCALL 32
 #endif
+
+static void
+gmp_rrandomb (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
 {
   int nb;
-  int bit_pos;
-  mp_size_t limb_pos;
-  mp_limb_t ran, ranm;
-  mp_limb_t acc;
-  mp_size_t n;
+  int bit_pos;			/* bit number of least significant bit where
+				   next bit field to be inserted */
+  mp_size_t ri;			/* index in rp */
+  mp_limb_t ran, ranm;		/* buffer for random bits */
+  mp_limb_t acc;		/* accumulate output random data here */
+  int ran_nbits;		/* number of valid bits in ran */
 
-  bit_pos = nbits % BITS_PER_MP_LIMB;
-  limb_pos = nbits / BITS_PER_MP_LIMB;
-  if (bit_pos == 0)
-    {
-      bit_pos = BITS_PER_MP_LIMB;
-      limb_pos--;
-    }
+  ran_nbits = 0;
+  bit_pos = (nbits - 1) % GMP_NUMB_BITS;
+  ri = (nbits - 1) / GMP_NUMB_BITS;
 
   acc = 0;
-  while (limb_pos >= 0)
+  while (ri >= 0)
     {
-      _gmp_rand (&ranm, rstate, BITS_PER_CHUNK + 1);
-      ran = ranm;
-      nb = (ran >> 1) + 1;
+      if (ran_nbits < LOGBITS_PER_BLOCK + 1)
+	{
+	  _gmp_rand (&ranm, rstate, BITS_PER_RANDCALL);
+	  ran = ranm;
+	  ran_nbits = BITS_PER_RANDCALL;
+	}
+
+      nb = (ran >> 1) % (1 << LOGBITS_PER_BLOCK) + 1;
       if ((ran & 1) != 0)
 	{
-	  /* Generate a string of ones.  */
+	  /* Generate a string of nb ones.  */
 	  if (nb > bit_pos)
 	    {
-	      rp[limb_pos--] = acc | ((((mp_limb_t) 1) << bit_pos) - 1);
-	      bit_pos += BITS_PER_MP_LIMB;
+	      rp[ri--] = acc | (((mp_limb_t) 2 << bit_pos) - 1);
+	      bit_pos += GMP_NUMB_BITS;
 	      bit_pos -= nb;
-	      acc = (~(mp_limb_t) 0) << bit_pos;
+	      acc = ((~(mp_limb_t) 1) << bit_pos) & GMP_NUMB_MASK;
 	    }
 	  else
 	    {
 	      bit_pos -= nb;
-	      acc |= ((((mp_limb_t) 1) << nb) - 1) << bit_pos;
+	      acc |= (((mp_limb_t) 2 << nb) - 2) << bit_pos;
 	    }
 	}
       else
 	{
-	  /* Generate a string of zeroes.  */
+	  /* Generate a string of nb zeroes.  */
 	  if (nb > bit_pos)
 	    {
-	      rp[limb_pos--] = acc;
+	      rp[ri--] = acc;
 	      acc = 0;
-	      bit_pos += BITS_PER_MP_LIMB;
+	      bit_pos += GMP_NUMB_BITS;
 	    }
 	  bit_pos -= nb;
 	}
+      ran_nbits -= LOGBITS_PER_BLOCK + 1;
+      ran >>= LOGBITS_PER_BLOCK + 1;
     }
 }