home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume3 / xstatic / part01 / genbits.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-20  |  6.9 KB  |  199 lines
  1. /*
  2. ** Portions stolen from 'xphoon' with grateful acknowledgement AND
  3. ** Portions written by Tim Bray are
  4. **  Copyright (C) 1988 by Jef Poskanzer and Craig Leres.
  5. **
  6. ** Permission to use, copy, modify, and distribute this software and its
  7. ** documentation for any purpose and without fee is hereby granted, provided
  8. ** that the above copyright notice appear in all copies and that both that
  9. ** copyright notice and this permission notice appear in supporting
  10. ** documentation.  This software is provided "as is" without express or
  11. ** implied warranty.
  12. */
  13.  
  14. #include <X11/Xos.h>
  15.  
  16. /*
  17.  * genbits(p, bytecount, bits) - fills the char array 'bits', which is of size
  18.  *  'bytecount', with bits each of of which is on with probability not too far
  19.  *  from p.
  20.  *
  21.  * This is quite a lot of bits; at least a million.  Calling random() this
  22.  *  many times is slow, sez the profiler.  For this reason, the random bits 
  23.  *  are generated eight at a time.  It is easy to calculate the probability 
  24.  *  of N bits in a byte being on as a function of p.  An array of 1024 chars 
  25.  *  is filled with numbers from 0 to 8 inclusive in numbers proportional to 
  26.  *  the probability of a byte having that many bits on.  Ten bits of random 
  27.  *  number index this array and select the number of bits to be on in the next
  28.  *  byte.  Once the number of on-bits is selected, another random number is 
  29.  *  used to select which byte to use among all the bytes with that many bits 
  30.  *  on.  For fear that something might not be relatively prime with something
  31.  *  else, the same random number is not used for the two indexing operations - 
  32.  *  this is almost certainly unnecessary caution.
  33.  *
  34.  * There is a slight loss due to quantizing the probability domain into
  35.  *  1024 slots.  For example if p is less than about .43, this will never 
  36.  *  generate a byte of all 1's.  An earlier version, which called random() 
  37.  *  once per bit, would in fact do this a few times in a million, as the 
  38.  *  combinatorics say it should.  Note that the possible bit counts can be 
  39.  *  encoded in 4 bits, so you could quantize into 2048 or 4096 chunks without 
  40.  *  wasting too much memory.  However, I went cross-eyed alternating my screen
  41.  *  back and forth between the two versions and couldn't convince myself that
  42.  *  there was a difference.  And this version (with the tenbits() trick)
  43.  *  has 24 times fewer calls to random and is dramatically faster.  Have to
  44.  *  show this to a real mathematician who understands what a random number is
  45.  *  and dodge the flying vomit.
  46.  */
  47. /* Static precomputed data for genbits */
  48. /* bytes with this many bits on */
  49. static int BytesWithCount0[] = { 0 };
  50. static int BytesWithCount1[] = { 1, 2, 4, 8, 16, 32, 64, 128 };
  51. static int BytesWithCount2[] = 
  52. { 
  53.   3, 5, 6, 9, 10, 12, 17, 18, 20, 24, 33, 34, 36, 40, 48, 65, 66, 68, 72, 80,
  54.   96, 129, 130, 132, 136, 144, 160, 192
  55. };
  56. static int BytesWithCount3[] = 
  57. {
  58.   7, 11, 13, 14, 19, 21, 22, 25, 26, 28, 35, 37, 38, 41, 42, 44, 49, 50, 52,
  59.   56, 67, 69, 70, 73, 74, 76, 81, 82, 84, 88, 97, 98, 100, 104, 112, 131, 133,
  60.   134, 137, 138, 140, 145, 146, 148, 152, 161, 162, 164, 168, 176, 193, 194,
  61.   196, 200, 208, 224
  62. };
  63. static int BytesWithCount4[] = 
  64. {
  65.   15, 23, 27, 29, 30, 39, 43, 45, 46, 51, 53, 54, 57, 58, 60, 71, 75, 77, 78,
  66.   83, 85, 86, 89, 90, 92, 99, 101, 102, 105, 106, 108, 113, 114, 116, 120, 135,
  67.   139, 141, 142, 147, 149, 150, 153, 154, 156, 163, 165, 166, 169, 170, 172,
  68.   177, 178, 180, 184, 195, 197, 198, 201, 202, 204, 209, 210, 212, 216, 225,
  69.   226, 228, 232, 240
  70. };
  71. static int BytesWithCount5[] = 
  72. { 
  73.   31, 47, 55, 59, 61, 62, 79, 87, 91, 93, 94, 103, 107, 109, 110, 115, 117,
  74.   118, 121, 122, 124, 143, 151, 155, 157, 158, 167, 171, 173, 174, 179, 181,
  75.   182, 185, 186, 188, 199, 203, 205, 206, 211, 213, 214, 217, 218, 220, 227,
  76.   229, 230, 233, 234, 236, 241, 242, 244, 248
  77. };
  78. static int BytesWithCount6[] = 
  79. {
  80.   63, 95, 111, 119, 123, 125, 126, 159, 175, 183, 187, 189, 190, 207, 215, 219,
  81.   221, 222, 231, 235, 237, 238, 243, 245, 246, 249, 250, 252
  82. };
  83. static int BytesWithCount7[] = { 127, 191, 223, 239, 247, 251, 253, 254 };
  84. static int BytesWithCount8[] = { 255 };
  85.  
  86. static int * BytesWithCount[] = 
  87. {
  88.   BytesWithCount0, BytesWithCount1, BytesWithCount2, BytesWithCount3,
  89.   BytesWithCount4, BytesWithCount5, BytesWithCount6, BytesWithCount7,
  90.   BytesWithCount8 
  91. };
  92.  
  93. /* Eight choose whatever - also the sizes of the BytesWith arrays */
  94. static int EightChoose[] = { 1, 8, 28, 56, 70, 56, 28, 8, 1 };
  95.  
  96. extern time_t time();
  97. extern long   random();
  98.  
  99. void
  100. genbits(p, bytecount, bits)
  101.   double          p;
  102.   int             bytecount;
  103.   register char * bits;
  104. {
  105.   char          bitcounts[1024];
  106.   register int      count;
  107.   int          i;
  108.   int          j;
  109.   int          k;
  110.   register int      last;
  111.   int          maxbitcount;
  112.   int          PartOfK[9];
  113.   double      PToThe[9];
  114.   double      probability;
  115.   double      QToThe[9];
  116.   register int      selector;
  117.   register char * stop;
  118.  
  119.   /* Compute the probabilities of bit counts 0 through 8 */
  120.   /* first, work out the powers of p and 1-p */
  121.   PToThe[0] = QToThe[0] = 1.0;
  122.   for (i = 1; i <= 8; i++)
  123.   { /* for each bit count */
  124.     PToThe[i] = PToThe[i - 1] * p;
  125.     QToThe[i] = QToThe[i - 1] * (1.0 - p);
  126.   } /* for each bit count */
  127.  
  128.   /* the probability of i bits being on is (8 choose i) * p**i * q**(8-i) */
  129.   count = maxbitcount = 0;
  130.   for (i = 0; i <= 8; i++)
  131.   { /* for each bit count */
  132.     probability = ((double) EightChoose[i]) * PToThe[i] * QToThe[8 - i];
  133.     PartOfK[i] = (int) (1024.0 * probability);
  134.     if (PartOfK[i] > PartOfK[maxbitcount])
  135.       maxbitcount = i;
  136.     count += PartOfK[i];
  137.   } /* for each bit count */
  138.  
  139.   /* might not come to 1024 because of rounding error; correct largest */
  140.   if (count != 1024)
  141.     PartOfK[maxbitcount] += 1024 - count;
  142.  
  143.   /* Load up bitcount array with appropriate numbers of different counts */
  144.   for(i = 0, j = 0; i <= 8; i++)
  145.     for (k = 0; k < PartOfK[i]; k++)
  146.       bitcounts[j++] = i;
  147.     
  148.   /* Set up and let 'er rip */
  149.   srandom((int) time((time_t *) 0));
  150.   last = tenbits();
  151.   stop = bits + bytecount;
  152.   while (bits < stop)
  153.   { /* for each byte */
  154.  
  155.     /* select bit count depending on probability */
  156.     selector = tenbits();
  157.     count = bitcounts[selector];
  158.  
  159.     /* use a randomly selected byte with that many bits on */
  160.     if (count == 0)
  161.       count = 0;
  162.     else if (count == 8)
  163.       count = 0xff;
  164.     else
  165.       count = BytesWithCount[count][last % EightChoose[count]];
  166.     *bits++ = count;
  167.  
  168.     /* remember the random number to re-use next time */
  169.     last = selector;
  170.  
  171.   } /* for each byte */
  172. }
  173.  
  174. /*
  175.  * Returns ten random bits.  Since it is expensive to call random(), and we
  176.  *  only need ten bits per, and random gives us 31 bits per call, why waste
  177.  *  the other 20?
  178.  */
  179. static int 
  180. tenbits()
  181. {
  182.   static int onetwothree = 3;
  183.   static int ranbits;
  184.  
  185.   /* use random() every third call; shift 10 bits over on the other two */
  186.   if (onetwothree == 3)
  187.   {
  188.     ranbits = random();
  189.     onetwothree = 1;
  190.   }
  191.   else
  192.   {
  193.     ranbits >>= 10;
  194.     onetwothree++;
  195.   }
  196.  
  197.   return ranbits & 1023;
  198. }
  199.