dieharder-3.31.1/libdieharder/bits.cruft

/*
 *========================================================================
 * $Id: bits.c 245 2006-10-05 16:34:42Z rgb $
 *
 * See copyright in copyright.h and the accompanying file COPYING
 *========================================================================
 */

/*
 *========================================================================
 * This should be a nice, big case switch where we add EACH test
 * we might want to do and either just configure and do it or
 * prompt for input (if absolutely necessary) and then do it.
 *========================================================================
 */

#include <dieharder/libdieharder.h>

/*
 * This should be the only tool we use to access bit substrings
 * from now on.  Note that we write our own bitstring access tools
 * instead of using ldap's LBER (Basic Encoding Rules) library call
 * to both retain control and because it would introduce a slightly
 * exotic dependency in the resulting application.
 *
 * bstring is a pointer to the uint string to be parsed.  It is a uint
 * pointer to make it easy to pass arbitrary strings which will generally
 * be e.g. unsigned ints in dieharder but might be other data types
 * in other applications (might as well make this semi-portable while I'm
 * writing it at all).  bslen is the length of bitstring in uints.  blen is
 * the length of the bitstring to be returned (as an unsigned int) and has
 * to be less than the length, in bits, of bitstring.  Finally, boffset
 * is the bit index of the point in the bitstring from which the result
 * will be returned.
 *
 * The only other thing that has to be defined is the direction from which
 * the bit offset and bit length are counted.  We will make the
 * LEAST significant bit offset 0, and take the string from the direction
 * of increasing signicance.  Examples:
 *
 *   bitstring:  10010110, length 1 (byte, char).
 * for offset 0, length 4 this should return: 0110
 *     offset 1, length 4: 1011
 *     offset 5, length 4: 0100 (note periodic wraparound)
 *     offset 6, length 4: 1010 (ditto)
 *
 * where of course the strings are actually returned as e.g.
 *     00000000000000000000000000000110  (unsigned int).
 */
uint get_bit_ntuple(uint *bitstring,uint bslen,uint blen,uint boffset)
{

 unsigned int b,rlen;
 int ioffset;
 unsigned int result,carry,nmask;

 /*
  * Some tests for failure or out of bounds.  8*blen has to be less than
  * sizeof(uint).
  */
 if(blen > 8*sizeof(uint)) blen = 8*sizeof(uint);
 /*
  * Now that blen is sane, we can make a mask for the eventual return
  * value of length blen.
  */
 nmask = 1;
 /* dumpbits(&nmask,32); */
 for(b=0;b<blen-1;b++) {
   nmask = nmask <<1;
   nmask++;
   /* dumpbits(&nmask,32); */
 }
 /* dumpbits(&nmask,32); */

 if(verbose == D_BITS || verbose == D_ALL){
   printf("# get_bit_ntuple(): bslen = %u, blen = %u, boffset = %u\n",bslen,blen,boffset);
   printf("# get_bit_ntuple(): bitstring (uint) = %u\n",*bitstring);
   printf("# get_bit_ntuple(): bitstring = ");
   dumpbits(&bitstring[0],32);
   printf("# get_bit_ntuple(): nmask     = ");
   dumpbits(&nmask,32);
 }

 /*
  * This is the index of bitstring[] containing the first bit to
  * be returned, for example if boffset is 30, rmax_bits is 24,
  * and bslen (the length of the uint bitstring) is 4 (indices
  * run from 0-3) then ioffset is 4 - 1 -1 = 2 and the 30th bit
  * from the RIGHT is to be found in bitstring[2]. Put this uint
  * into result to work on further.
  */
 ioffset = bslen - (uint) boffset/rmax_bits - 1;
 result = bitstring[ioffset];
 if(verbose == D_BITS || verbose == D_ALL){
   printf("bitstring[%d] = %u\n",ioffset,result);
   printf("Initial result =          ");
   dumpbits(&result,32);
 }
 /*
  * Now, WHICH bit in result is the boffset relative to ITS
  * rightmost bit?  We have to find the modulus boffset%rmax_bits.
  * For example 30%24 = 6, so in the continuing example it would
  * be bit 6 in result.
  */
 boffset = boffset%rmax_bits;
 if(verbose == D_BITS || verbose == D_ALL){
   printf("Shifting to bit offset %u\n",boffset);
 }

 /*
  * Now for the easy part.  We shift right boffset bits.
  */
 for(b=0;b<boffset;b++) result = result >> 1;
 if(verbose == D_BITS || verbose == D_ALL){
   printf("Shifted result =          ");
   dumpbits(&result,32);
 }

 /*
  * rlen is the cumulated length of result.  Initially, it is set to
  * rmax_bits - boffset, the number of bits result can now contribute from
  * boffset.  We now have to loop, adding bits from uints up the list
  * (cyclic) until we get enough to return blen.  Note that we have to
  * loop because (for example) blen = 32, rmax_bits = 16, boffset = 30
  * would start in bitstring[2] and get 2 bits (30 and 31), get all 16 bits
  * from bitstring[1], and still need 12 bits of bitstring[0] to return.
  */
 rlen = rmax_bits - boffset;
 if(verbose == D_BITS || verbose == D_ALL){
   printf("Cumulated %u signifcant bits\n",rlen);
 }

 while(blen > rlen){
   /*
    * If we get here, we have to use either bitstring[ioffset-1] or
    * bitstring[bslen-1], shifted and filled into result starting
    * at the correct slot.  Put this into carry to work on.
    */
   ioffset--;
   if(ioffset < 0) ioffset = bslen-1;
   carry = bitstring[ioffset];
   if(verbose == D_BITS || verbose == D_ALL){
     printf("bitstring[%d] = %u\n",ioffset,carry);
     printf("Next carry =              ");
     dumpbits(&carry,32);
   }

   /*
    * This is tricky!  Shift carry left until the first bit lines
    * up with the first empty bit in result, which we'll hope is
    * the current rlen bit.
    */
   for(b=0;b<rlen;b++){
     carry = carry << 1;
   }
   if(verbose == D_BITS || verbose == D_ALL){
     printf("Shifted carry =           ");
     dumpbits(&carry,32);
   }

   /*
    * Now we simply add result and carry AND increment rlen by
    * rmax_bit (since this is the exact number of bits it adds
    * to rlen).
    */
   result += carry;
   rlen += rmax_bits;
   if(verbose == D_BITS || verbose == D_ALL){
     printf("Cumulated %u signifcant bits\n",rlen);
     printf("result+carry =            ");
     dumpbits(&result,32);
   }
 }

 result = result & nmask;
   if(verbose == D_BITS || verbose == D_ALL){
     printf("Returning Result =        ");
     dumpbits(&result,32);
     printf("==========================================================\n");
   }
 return(result);

}

/*
 * dumpbits only can dump 8*sizeof(unsigned int) bits at a time.
 */
void dumpbits(uint *data, unsigned int nbits)
{

 uint i,j;
 uint mask;

 if(nbits > 8*sizeof(unsigned int)) {
   nbits = 8*sizeof(unsigned int);
 }

 mask = (uint)pow(2,nbits-1);
 for(i=0;i<nbits;i++){
   if(verbose == -1){
     printf("\nmask = %u = %04x :",mask,mask);
   }
   j = (mask & *data)?1:0;
   printf("%1u",j);
   mask = mask >> 1;
 }
 printf("\n");

}

void cycle(unsigned int *data, unsigned int nbits)
{
 unsigned int i;
 unsigned int result,rbit,lmask,rmask;
 /*
  * We need two masks.  One is a mask of all the significant bits
  * in the rng, which may be as few as 24 and can easily be only
  * 31.
  *
  * The other is lmask, with a single bit in the position of the
  * leftmost significant bit.  We can make them both at once.
  */
 rmask = 1;
 lmask = 1;
 for(i=0;i<nbits-1;i++) {
  rmask = rmask << 1;
  rmask++;
  lmask = lmask << 1;
 }
 if(verbose){
   printf("%u bit rmask = ",nbits);
   dumpbits(&rmask,32);
   printf("%u bit lmask = ",nbits);
   dumpbits(&lmask,32);
 }
 /*
  * Next, we create mask the int being bit cycled into an internal
  * register.
  */
 result = *data & rmask;
 if(verbose){
   printf("Original int: ");
   dumpbits(data,32);
   printf("  Masked int: ");
   dumpbits(&result,32);
 }
 rbit = 1 & result;
 result = result >> 1;
 result += rbit*lmask;
 *data = result;
 if(verbose){
   printf(" Rotated int: ");
   dumpbits(data,32);
 }

}

/*
 * This is still a good idea, but we have to modify it so that it ONLY
 * gets VALID bits by their absolute index.
 */
int get_bit(uint *rand_uint, unsigned int n)
{

 int i;
 unsigned int index,offset,mask;
 static unsigned last_rand;

 /*
  * This routine is designed to get the nth VALID bit of an input uint
  * *rand_int.  The indexing is a bit tricky.  index tells us which vector
  * element contains the bit being sought:
  */
 index = (int) (n/rmax_bits);

 /*
  * Then we have to compute the offset of the bit desired, starting from
  * the first significant/valid bit in the unsigned int.
  */
 offset = (8*sizeof(unsigned int) - rmax_bits) + n%rmax_bits;
 mask = (int)pow(2,8*sizeof(unsigned int) - 1);
 mask = mask>>offset;
 if(mask & rand_uint[index]){
   return(1);
 } else {
   return(0);
 }

}

int get_int_bit(uint i, uint n)
{

 unsigned int mask;

 /*
  * This routine gets the nth bit of the unsigned int i.  It does very
  * limited checking to ensure that n is in the range 0-sizeof(uint)
  * Note
  */
 if(n < 0 || n > 8*sizeof(uint)){
   fprintf(stderr,"Error: bit offset %d exceeds length %d of uint.\n",n,8*sizeof(uint));
   exit(0);
 }


 /*
  * Then we have make a mask and shift it over from the first (least
  * significant) bit in the unsigned int.  AND the result with i and
  * we're done.
  */
 mask = 1;
 mask = mask<<n;
 /* dumpbits(&mask,32); */
 /* dumpbits(&i,32); */
 if(mask & i){
   return(1);
 } else {
   return(0);
 }

}

/*
 * dumpbits only can dump 8*sizeof(unsigned int) bits at a time.
 */
void dumpbits_left(unsigned int *data, unsigned int nbits)
{

 int i,j;
 unsigned int mask;

 if(nbits > 8*sizeof(unsigned int)) {
   nbits = 8*sizeof(unsigned int);
 }

 mask = 1;
 for(i=0;i<nbits;i++){
   if(mask & *data){
     printf("1");
   } else {
     printf("0");
   }
   mask = mask << 1;
 }
 printf("\n");
}


unsigned int bit2uint(char *abit,unsigned int blen)
{

 int i,bit;
 unsigned int result;

 /* Debugging
 if(verbose == D_BITS || verbose == D_ALL){
   printf("# bit2uint(): converting %s\n",abit);
 }
 */

 result = 0;
 for(i = 0; i < blen; i++){
   result = result << 1;
   bit = abit[i] - '0';
   result += bit;
   /* Debugging
   if(verbose == D_BITS || verbose == D_ALL){
     printf("# bit2uint(): bit[%d] = %d, result = %u\n",i,bit,result);
   }
   */
 }

 /* Debugging
 if(verbose == D_BITS || verbose == D_ALL){
   printf("# bit2uint(): returning %0X\n",result);
 }
 */
 return(result);

}

void fill_uint_buffer(uint *data,uint buflength)
{

 /*
  * This routine fills *data with random bits from the current
  * random number generator.  Note that MANY generators return
  * fewer than 32 bits, making this a bit of a pain in the ass.
  * We need buffers like this for several tests, though, so it
  * is worth it to create a routine to do this once and for all.
  */

 uint bufbits,bcnt,bdelta;
 uint i,tmp1,tmp2,mask;

 /*
  * Number of bits we must generate.
  */
 bufbits = buflength*sizeof(uint)*CHAR_BIT;
 bdelta = sizeof(uint)*CHAR_BIT - rmax_bits;
 mask = 0;
 for(i=0;i<bdelta;i++) {
  mask = mask<<1;
  mask++;
 }
 if(verbose == D_BITS || verbose == D_ALL){
   printf("rmax_bits = %d  bdelta = %d\n",rmax_bits,bdelta);
 }

 for(i=0;i<buflength;i++){

   /* put rmax_bits into tmp1 */
   tmp1 = gsl_rng_get(rng);
   /* Cruft
   printf("tmp1 = %10u = ",tmp1);
   dumpbits(&tmp1,32);
   */

   /* Shift it over to the left */
   tmp1 = tmp1 << bdelta;
   /*
   printf("tmp1 = %10u = ",tmp1);
   dumpbits(&tmp1,32);
   */

   /* put rmax_bits into tmp2 */
   tmp2 = gsl_rng_get(rng);
   /*
   printf("tmp2 = %10u = ",tmp2);
   printf("mask = %10u = ",mask);
   dumpbits(&mask,32);
   */

   /* mask the second rand */
   tmp2 = tmp2&mask;
   /*
   printf("tmp2 = %10u = ",tmp2);
   dumpbits(&tmp2,32);
   */

   /* Fill in the rest of the uint */
   tmp1 = tmp1 + tmp2;
   /* Cruft
   printf("tmp1 = %10u = ",tmp1);
   dumpbits(&tmp1,32);
   */

   data[i] = tmp1;
 }

}

/*
 * OK, the routines above work but they suck.  We need a set of SMALL
 * bitlevel routines for manipulating, aggregating, windowing and so
 * on, especially if we want to write extended versions of the bit
 * distribution test.
 *
 * On that note, let's generate SMALL routines for:
 *
 *  a) creating a uint mask() to select bits from uints
 *
 *  b) grabbing a masked window of bits and left/right shifting it
 *     to an arbitrary offset within the uint (no wraparound).
 *
 *  c) conjoining two such masked objects to produce a new object.
 *
 *  d) filling a buffer of arbitrary length with sequential bits from
 *     the selected rng.
 *
 *  e) Selecting a given window from that buffer with a given bitwise offset
 *     and with periodic wraparound.
 *
 *  f).... (we'll see what we need)
 */


/*
 * This generates a uint-sized mask of 1's starting at bit position
 * bstart FROM THE LEFT (with 0 being the most significant/sign bit)
 * and ending at bit position bstop.
 *
 * That is:
 *
 * umask(3,9) generates a uint containing (first line indicates bit
 * positions only):
 *  01234567890123456789012345678901
 *  00011111110000000000000000000000
 */
uint b_umask(uint bstart,uint bstop)
{

 uint b,mask,blen;

 if(bstart < 0 || bstop > 31 || bstop < bstart){
   printf("b_umask() error: bstart <= bstop must be in range 0-31.\n");
   exit(0);
 }
 blen = bstop-bstart+1;

 /*
  * Create blen 1's on right
  */
 mask = 1;
 for(b=1;b<blen;b++) {
   mask = mask <<1;
   mask++;
   /* dumpbits(&mask,sizeof(uint)*CHAR_BIT); */
 }

 /*
  * Now shift them over to the correct starting point.
  */
 mask = mask << (32-blen-bstart);
 /* dumpbits(&mask,sizeof(uint)*CHAR_BIT); */

 return mask;

}

/*
 * This uses b_umask to grab a particular window's worth of bits
 * from an arbitrary uint and THEN shift it to the new desired offset.
 * bstart FROM THE LEFT (with 0 being the most significant/sign bit)
 * and ending at bit position bstop.
 *
 * That is:
 *
 * b_window(input,2,5,0) generates a uint containing (first line indicates bit
 * positions only):
 *  input 01101011010000111010101001110110
 *   mask 00111100000000000000000000000000
 *      & 00101000000000000000000000000000
 *  shift 10100000000000000000000000000000
 */
uint b_window(uint input,uint bstart,uint bstop,uint boffset)
{

 uint b,mask,output;
 int shift;

 if(bstart < 0 || bstop > 31 || bstop < bstart){
   printf("b_umask() error: bstart <= bstop must be in range 0-31.\n");
   exit(0);
 }
 if(boffset < 0 || boffset > 31){
   printf("b_window() error: boffset must be in range 0-31.\n");
   exit(0);
 }
 shift = bstart - boffset;

 /* dumpbits(&input,sizeof(uint)*CHAR_BIT); */
 mask = b_umask(bstart,bstop);
 /* dumpbits(&mask,sizeof(uint)*CHAR_BIT); */
 output = input & mask;
 /* dumpbits(&output,sizeof(uint)*CHAR_BIT); */
 if(shift>0){
   output = output << shift;
 } else {
   output = output >> (-shift);
 }
 /* dumpbits(&output,sizeof(uint)*CHAR_BIT); */

 return output;

}

/*
 * Rotate the uint left (with periodic BCs)
 */
uint b_rotate_left(uint input,uint shift)
{

 uint tmp;

 dumpbits(&input,sizeof(uint)*CHAR_BIT);
 tmp = b_window(input,0,shift-1,32-shift);
 dumpbits(&tmp,sizeof(uint)*CHAR_BIT);
 input = input << shift;
 dumpbits(&input,sizeof(uint)*CHAR_BIT);
 input += tmp;
 dumpbits(&input,sizeof(uint)*CHAR_BIT);

 return input;

}
/*
 * Rotate the uint right (with periodic BCs)
 */
uint b_rotate_right(uint input, uint shift)
{

 uint tmp;

 dumpbits(&input,sizeof(uint)*CHAR_BIT);
 tmp = b_window(input,32-shift,31,0);
 dumpbits(&tmp,sizeof(uint)*CHAR_BIT);
 input = input >> shift;
 dumpbits(&input,sizeof(uint)*CHAR_BIT);
 input += tmp;
 dumpbits(&input,sizeof(uint)*CHAR_BIT);

 return input;

}

/*
 * OK, with this all in hand, we can NOW write routines to return
 * pretty much any sort of string of bits from the prevailing rng
 * without too much effort. Let's get an ntuple from a uint vector
 * of arbitrary length and offset, with cyclic boundary conditions.
 *
 * We have to pack the answer into the LEAST significant bits in the
 * output vector, BTW, not the MOST.  That is, we have to fill the
 * output window all the way to the rightmost bit.  Tricky.
 */
void get_ntuple_cyclic(uint *input,uint ilen,
    uint *output,uint olen,uint ntuple,uint offset)
{

 int istart,iend,istartbit,iendbit,ostart,oend,ostartbit,oendbit;
 int i,j,bcnt,bdelta,uintlen;

 uintlen = sizeof(uint)*CHAR_BIT;

 /*
  * We start somewhere in input[istart] and end in input[iend].
  * Ditto output[ostart], output[oend].  We also need the BIT offsets
  * in these uints.
  *
  * SO, suppose offset = 40, ntuple = 40.  Then istart = 40/32 = 1,
  * iend = 80/32 = 2, istartbit = 40 - 32 = 8, iendbit = 96 - 80 = 16
  * (where this is one bit PAST the end of the window!)
  * The slice thus runs from:
  *   index[1], bits 8-31 (total of 24 bits)
  * to:
  *   index[2], bits 0-15 (total of 16 bits)
  *
  * for a total of 40 bits starting at offset 1*32+8 = 40.  We can deal
  * with periodic wraparound as we work though index by using i = i%ilen
  * and trying to reach iend = iend%ilen;
  */
 istart = offset/uintlen;
 iend = (offset+ntuple)/uintlen;
 istartbit = offset - istart*uintlen;
 iendbit = (istartbit + ntuple)%uintlen;
 /* We want to point at the bit PAST the end */
 if(iendbit == 0){
   iend--;
   iendbit = uintlen;
 }
 /*
  * To accomplish periodic wraparound, we mod the
  * istart and iend markers (and later, i itself).
  */
 istart = istart%ilen;
 iend = iend%ilen;
 printf("istart = %u, iend = %u, istartbit = %u, iendbit = %u\n",
   istart,iend,istartbit,iendbit);
 /*
  * Now we have to do the same thing for the output vector.  Here life
  * is a BIT better, as we can just copy the windows above straight over
  * and shift things at the end.  So what we really have to do is figure
  * out ostart.  Assuming that olen = 4, ostart in our example should be:
  *  ostart = 4 - 40/32 - 1 = 2,
  * so the answer will go in output[2] and output[3].
  *
  *
  */
 if((ntuple%uintlen) == 0){
   ostart = olen - ntuple/uintlen;
 } else {
   ostart = olen - ntuple/uintlen - 1;
 }
 oend = olen - 1;
 printf("ostart = %u, oend = %u,\n",ostart,oend);

 /*
  * Now we do a double loop, iterating BACKWARDS.
  */
 i = iend;
 j = oend;

 /*
  * The only sane way to do this with periodic wraparound is to just
  * count the bits returned and the chunk sizes.  There are only four
  * (possible) chunks we might grab:
  */
 left = iendbit;
 right = uintlen - iendbit;  /* May be zero */
 /*
  * This is a possible chunk, depending on whether or not it is
  * less than right chunk.
  */
 if(istartbit > iendbit) {
   topright = uintlen - istartbit;
   topleft = 0;
 } else {
   topright = right;
   topleft = iendbit - istartbit;
 }

 bcnt = 0;
 while(bcnt < ntuple){
   /*
    * Second (least significant) part of output[j].
    */
   printf("Starting Fill Right: i = %d j = %d bcnt = %d\n",i,j,bcnt);

   /* Fill right to terminate */
   if((uintlen-istartbit) == (ntuple-bcnt)){
     bdelta = uintlen-istartbit;
     printf("3 Fill right from %d: output[%d] = ",i,j);
     output[j] = b_window(input[i],istartbit,uintlen-1,bdelta);
     bcnt += bdelta;
     dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
     break;
   }

   /* Fill right to terminate */
   if((ntuple-bcnt) == (iendbit-istartbit)){
     printf("2 Fill right from %d: output[%d] = ",i,j);
     output[j] = b_window(input[i],istartbit,iendbit-1,uintlen-iendbit+istartbit);
     bcnt += iendbit-istartbit;
     dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
     break;
   }

   /* Plain old fill right (all other cases) */
   printf("1 Fill right from %d: output[%d] = ",i,j);
   output[j] = b_window(input[i],0,iendbit-1,uintlen-iendbit);
   bcnt += iendbit;
   dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
   printf("Finished Fill Right: i = %d j = %d bcnt = %d\n",i,j,bcnt);
   if(bcnt == ntuple) break;

   /*
    * Now we check for the First (most significant) part of output[j].
    * It MUST come from the previous input line if any, so we decrement
    * i.
    */
   i--;
   if(i == -1) i = ilen-1;
   printf("Starting Fill Left: i = %d j = %d bcnt = %d\n",i,j,bcnt);

   if((ntuple-bcnt) > iendbit){
     printf("1 Fill left from %d: output[%d] = ",i,j);
     output[j] += b_window(input[i],iendbit,uintlen-1,0);
     bcnt += uintlen-iendbit;
     dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
   }

   /* Fill left to terminate */
   if((ntuple-bcnt) == (uintlen-iendbit)){
      printf("2 Fill left from %d: output[%d] = ",i,j);
      output[j] += b_window(input[i],iendbit,uintlen-1,0);
      bcnt += uintlen-iendbit;
      dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
      break;
   }

   /* Fill left to terminate */
   if((ntuple-bcnt) == (uintlen-istartbit)){
     printf("3 Fill left from %d: output[%d] = ",i,j);
     output[j] += b_window(input[i],istartbit,uintlen-1,istartbit-iendbit);
     bcnt += uintlen-istartbit;
     dumpbits(&output[j],sizeof(uint)*CHAR_BIT);
     break;
   }

   printf("Finished Fill Left: i = %d j = %d bcnt = %d\n",i,j,bcnt);
   j--;

 }

}