home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / g__lib / acg.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-23  |  8.1 KB  |  263 lines
  1. #include <ACG.h>
  2. #include <assert.h>
  3.  
  4. //
  5. //    This is an extension of the older implementation of Algorithm M
  6. //    which I previously supplied. The main difference between this
  7. //    version and the old code are:
  8. //
  9. //        + Andres searched high & low for good constants for
  10. //          the LCG.
  11. //
  12. //        + theres more bit chopping going on.
  13. //
  14. //    The following contains his comments.
  15. //
  16. //    agn@UNH.CS.CMU.EDU sez..
  17. //    
  18. //    The generator below is based on 2 well known
  19. //    methods: Linear Congruential (LCGs) and Additive
  20. //    Congruential generators (ACGs).
  21. //    
  22. //    The LCG produces the longest possible sequence
  23. //    of 32 bit random numbers, each being unique in
  24. //    that sequence (it has only 32 bits of state).
  25. //    It suffers from 2 problems: a) Independence
  26. //    isnt great, that is the (n+1)th number is
  27. //    somewhat related to the preceding one, unlike
  28. //    flipping a coin where knowing the past outcomes
  29. //    dont help to predict the next result.  b)
  30. //    Taking parts of a LCG generated number can be
  31. //    quite non-random: for example, looking at only
  32. //    the least significant byte gives a permuted
  33. //    8-bit counter (that has a period length of only
  34. //    256).  The advantage of an LCA is that it is
  35. //    perfectly uniform when run for the entire period
  36. //    length (and very uniform for smaller sequences
  37. //    too, if the parameters are chosen carefully).
  38. //    
  39. //    ACGs have extremly long period lengths and
  40. //    provide good independence.  Unfortunately,
  41. //    uniformity isnt not too great. Furthermore, I
  42. //    didnt find any theoretically analysis of ACGs
  43. //    that addresses uniformity.
  44. //    
  45. //    The RNG given below will return numbers
  46. //    generated by an LCA that are permuted under
  47. //    control of a ACG. 2 permutations take place: the
  48. //    4 bytes of one LCG generated number are
  49. //    subjected to one of 16 permutations selected by
  50. //    4 bits of the ACG. The permutation a such that
  51. //    byte of the result may come from each byte of
  52. //    the LCG number. This effectively destroys the
  53. //    structure within a word. Finally, the sequence
  54. //    of such numbers is permuted within a range of
  55. //    256 numbers. This greatly improves independence.
  56. //    
  57. //
  58. //  Algorithm M as describes in Knuths "Art of Computer Programming",
  59. //    Vol 2. 1969
  60. //  is used with a linear congruential generator (to get a good uniform
  61. //  distribution) that is permuted with a Fibonacci additive congruential
  62. //  generator to get good independence.
  63. //
  64. //  Bit, byte, and word distributions were extensively tested and pass
  65. //  Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity
  66. //  assumption holds with probability > 0.999)
  67. //
  68. //  Run-up tests for on 7E8 numbers confirm independence with
  69. //  probability > 0.97.
  70. //
  71. //  Plotting random points in 2d reveals no apparent structure.
  72. //
  73. //  Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i),
  74. //    i=1..512)
  75. //  results in no obvious structure (A(i) ~ const).
  76. //
  77. //  Except for speed and memory requirements, this generator outperforms
  78. //  random() for all tests. (random() scored rather low on uniformity tests,
  79. //  while independence test differences were less dramatic).
  80. //
  81. //  AGN would like to..
  82. //  thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help.
  83. //
  84. //  And I would (DGC) would like to thank Donald Kunth for AGN for letting me
  85. //  use his extensions in this implementation.
  86. //
  87.  
  88. //
  89. //    Part of the table on page 28 of Knuth, vol II. This allows us
  90. //    to adjust the size of the table at the expense of shorter sequences.
  91. //
  92.  
  93. static randomStateTable[][3] = {
  94. {3,7,16}, {4,9, 32}, {3,10, 32}, {1,11, 32}, {1,15,64}, {3,17,128},
  95. {7,18,128}, {3,20,128}, {2,21, 128}, {1,22, 128}, {5,23, 128}, {3,25, 128},
  96. {2,29, 128}, {3,31, 128}, {13,33, 256}, {2,35, 256}, {11,36, 256},
  97. {14,39,256}, {3,41,256}, {9,49,256}, {3,52,256}, {24,55,256}, {7,57, 256},
  98. {19,58,256}, {38,89,512}, {17,95,512}, {6,97,512}, {11,98,512}, {-1,-1,-1} };
  99.  
  100. //
  101. // spatial permutation table
  102. //    RANDOM_PERM_SIZE must be a power of two
  103. //
  104.  
  105. #define RANDOM_PERM_SIZE 64
  106. unsigned long randomPermutations[RANDOM_PERM_SIZE] = {
  107. 0xffffffff, 0x00000000,  0x00000000,  0x00000000,  // 3210
  108. 0x0000ffff, 0x00ff0000,  0x00000000,  0xff000000,  // 2310
  109. 0xff0000ff, 0x0000ff00,  0x00000000,  0x00ff0000,  // 3120
  110. 0x00ff00ff, 0x00000000,  0xff00ff00,  0x00000000,  // 1230
  111.  
  112. 0xffff0000, 0x000000ff,  0x00000000,  0x0000ff00,  // 3201
  113. 0x00000000, 0x00ff00ff,  0x00000000,  0xff00ff00,  // 2301
  114. 0xff000000, 0x00000000,  0x000000ff,  0x00ffff00,  // 3102
  115. 0x00000000, 0x00000000,  0x00000000,  0xffffffff,  // 2103
  116.  
  117. 0xff00ff00, 0x00000000,  0x00ff00ff,  0x00000000,  // 3012
  118. 0x0000ff00, 0x00000000,  0x00ff0000,  0xff0000ff,  // 2013
  119. 0x00000000, 0x00000000,  0xffffffff,  0x00000000,  // 1032
  120. 0x00000000, 0x0000ff00,  0xffff0000,  0x000000ff,  // 1023
  121.  
  122. 0x00000000, 0xffffffff,  0x00000000,  0x00000000,  // 0321
  123. 0x00ffff00, 0xff000000,  0x00000000,  0x000000ff,  // 0213
  124. 0x00000000, 0xff000000,  0x0000ffff,  0x00ff0000,  // 0132
  125. 0x00000000, 0xff00ff00,  0x00000000,  0x00ff00ff   // 0123
  126. };
  127.  
  128. //
  129. //    SEED_TABLE_SIZE must be a power of 2
  130. //
  131. #define SEED_TABLE_SIZE 32
  132. static unsigned long seedTable[SEED_TABLE_SIZE] = {
  133. 0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b,
  134. 0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf,
  135. 0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706,
  136. 0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10,
  137. 0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a,
  138. 0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32,
  139. 0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f,
  140. 0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf
  141. };
  142.  
  143. //
  144. //    The LCG used to scramble the ACG
  145. //
  146. //
  147. // LC-parameter selection follows recommendations in 
  148. // "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi.
  149. //
  150. // LC_A = 251^2, ~= sqrt(2^32) = 66049
  151. // LC_C = result of a long trial & error series = 3907864577
  152. //
  153.  
  154. static const unsigned long LC_A = 66049;
  155. static const unsigned long LC_C = 3907864577;
  156. inline LCG(unsigned long x)
  157. {
  158.     return( x * LC_A + LC_C );
  159. }
  160.  
  161.  
  162. ACG::ACG(unsigned long seed, int size)
  163. {
  164.     register unsigned long u;
  165.     
  166.     if (seed > -1 && seed < SEED_TABLE_SIZE) {
  167.     u = seedTable[seed];
  168.     } else {
  169.     u = seed ^ seedTable[seed & (SEED_TABLE_SIZE-1)];
  170.     }
  171.     
  172.     //
  173.     //    Determine the size of the state table
  174.     //
  175.     
  176.     for (register int l = 0;
  177.      randomStateTable[l][0] != -1 && randomStateTable[l][1] < size;
  178.      l++);
  179.     
  180.     if (randomStateTable[l][1] == -1) {
  181.     l--;
  182.     }
  183.     
  184.     stateSize = randomStateTable[l][1];
  185.     auxSize = randomStateTable[l][2];
  186.     
  187.     j = randomStateTable[l][0] - 1;
  188.     k = randomStateTable[l][1] - 1;
  189.     
  190.     //
  191.     //    Allocate the state table & the auxillary table in a single malloc
  192.     //
  193.     
  194.     state = new unsigned long[stateSize + auxSize];
  195.     auxState = &state[stateSize];
  196.     
  197.     //
  198.     //    Initialize the state
  199.     //
  200.     register int i;
  201.     for(i = 0; i < stateSize; i++) {
  202.     state[i] = u = LCG(u);
  203.     }
  204.     
  205.     for (i = 0; i < auxSize; i++) {
  206.     auxState[i] = u = LCG(u);
  207.     }
  208.     
  209.     k = u % stateSize;
  210.     int tailBehind = (stateSize - randomStateTable[l][0]);
  211.     j = k - tailBehind;
  212.     if (j < 0) {
  213.     j += stateSize;
  214.     }
  215.     
  216.     //    k = randomStateTable[l][1];
  217.     //    j = randomStateTable[l][0];
  218.     
  219.     lcgRecurr = u;
  220.     
  221.     assert(sizeof(double) == 2 * sizeof(long));
  222. }
  223.  
  224. ACG::~ACG()
  225. {
  226.     if (state) delete state;
  227.     state = 0;
  228.     // don't delete auxState, it's really an alias for state.
  229. }
  230.  
  231. //
  232. //    Returns 32 bits of random information.
  233. //
  234.  
  235. unsigned long ACG::asLong()
  236. {
  237.     unsigned long result = state[k] + state[j];
  238.     state[k] = result;
  239.     j = (j <= 0) ? (stateSize-1) : (j-1);
  240.     k = (k <= 0) ? (stateSize-1) : (k-1);
  241.     
  242.     short int auxIndex = (result >> 24) & (auxSize - 1);
  243.     register unsigned long auxACG = auxState[auxIndex];
  244.     auxState[auxIndex] = lcgRecurr = LCG(lcgRecurr);
  245.     
  246.     //
  247.     // 3c is a magic number. We are doing four masks here, so we
  248.     // do not want to run off the end of the permutation table.
  249.     // This insures that we have always got four entries left.
  250.     //
  251.     register unsigned long *perm = & randomPermutations[result & 0x3c];
  252.     
  253.     result =  *(perm++) & auxACG;
  254.     result |= *(perm++) & ((auxACG << 24)
  255.                | ((auxACG >> 8)& 0xffffff));
  256.     result |= *(perm++) & ((auxACG << 16)
  257.                | ((auxACG >> 16) & 0xffff));
  258.     result |= *(perm++) & ((auxACG <<  8)
  259.                | ((auxACG >> 24) &   0xff));
  260.     
  261.     return(result);
  262. }
  263.