home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume28 / mrandom-3.0 / part01 / src / bentley.c next >
Encoding:
C/C++ Source or Header  |  1994-05-06  |  3.9 KB  |  159 lines
  1. /* bentley.c 3.1 5/28/93 */
  2. /*  This is mrandom's RNG Algorithm #2                    */
  3.  
  4. /*  Original Author: Jon Bentley                    */
  5. /*    (Source code obtained from DIMACS shared account, March 1992.)*/
  6. /*  Header file/mrandom interface by Clark Thomborson, May 1992     */
  7. /*  Vectorized interface by Clark Thomborson, June 1992            */
  8. /*  Updated mrandom interface by Robert Plotkin, Feb/March 1993     */
  9. /*                                                                  */
  10. /*  This file contains a set of c-language functions for generating */
  11. /*  uniform integers.   This is a COMPLETELY PORTABLE generator.    */
  12. /*  It will give IDENTICAL sequences of random numbers for any      */
  13. /*  architecture with at least 30-bit integers, regardless of the   */
  14. /*  integer representation, MAXINT value, or roundoff/truncation    */
  15. /*  method, etc.                                                    */
  16.  
  17. /*  This Truly Remarkable RNG is described more fully in            */
  18. /*  J. Bentley's column, ``The Software Exploratorium ''            */
  19. /*  to appear in Unix Review in 1991.                               */ 
  20. /*  It is based on one in Knuth, Vol 2, Section 3.2.2 (Algorithm A) */ 
  21.  
  22. #include "bentley.h"
  23.  
  24. /*----RNG Global Variables-------*/ 
  25.  
  26. #define arr rngstate
  27. #define a (rngstate[RNGstatesize_2 - 2])
  28. #define b (rngstate[RNGstatesize_2 - 1])
  29.  
  30. /*********************************/
  31. /* External interface to mrandom */
  32. /*********************************/
  33. /* Generating procedure */
  34. long _lprand(RNGdata *rng)
  35. {
  36. return(lprand(RNGstate));
  37. }
  38.  
  39. /* Seeding procedure */
  40. void _lprand_seed(rng, seed)
  41. RNGdata *rng;
  42. long *seed;
  43. {
  44. sprand(*seed, RNGstate);
  45. }
  46.  
  47. /* Checking procedure */
  48. int _lprand_check(rng)
  49. RNGdata *rng;
  50. {
  51. return(ckrand(RNGstate));
  52. }
  53.  
  54.  
  55. /*----RNG Initializer------------*/
  56. /* Call once before using lprand */ 
  57.  
  58. extern void sprand (seed, rngstate)
  59. long seed; 
  60. long rngstate[RNGstatesize_2];
  61. {
  62.   int i, ii;
  63.   long last, next, v[55];
  64.   arr[0] = last = seed; 
  65.   next = 1;
  66.   for (i=1; i < 55; i++) {
  67.     ii = ( 21 * i ) % 55;
  68.     arr[ii] = next;
  69.     next = last - next; 
  70.     if (next < 0)
  71.       next += PRANDMAX;
  72.     last = arr[ii];
  73.   }
  74.  
  75.   a = 0; /* invariant: (b-a-24)%55 = 0 */
  76.   b = 24; 
  77.   /* cycle 165 times (more would probably be better -cdt) */
  78.   for (i=0;i<165;i++) {
  79.     lprand(rngstate);
  80.   }
  81. }
  82.  
  83. /*---------RNG---------------------*/
  84. /* Returns a long integer in the range 0...PRANDMAX-1 */
  85.  
  86. extern long lprand(rngstate)
  87. long rngstate[RNGstatesize_2];
  88. {
  89.   long t,ta,tb,retval;
  90.   
  91.   /* get local copies of a,b from state array */ 
  92.   ta = a;
  93.   tb = b;
  94.   
  95.   if (ta-- == 0) ta = 54;
  96.   if (tb-- == 0) tb = 54;
  97.   
  98.   t = arr[ta] - arr[tb];
  99.   if (t < 0) t+= PRANDMAX;
  100.   
  101.   arr[ta] = t;
  102.   retval = t;
  103.   
  104.   /* update state array */
  105.   a = ta;
  106.   b = tb;
  107.   
  108.   return(retval);
  109. }
  110.  
  111. /*---------CKRAND---------------------*/
  112. /* Returns 0 if rngstate is corrupted */
  113.  
  114. extern int ckrand(rngstate)
  115. long rngstate[RNGstatesize_2];
  116. {
  117.   if ( (b-a-24)%55 != 0 ) return(0); else return(1);
  118. }
  119.  
  120. /*-----------------------------------------------*/
  121. /* This is a little driver program so you can    */
  122. /* test the code.              */
  123. /* Typing: a.out 0 3 1         */
  124. /* should produce              */
  125. /*     921674862               */
  126. /*     250065336               */
  127. /*     377506581               */
  128. /*  Typing: a.out 1000000 1 2  */ 
  129. /*  should produce             */
  130. /*     57265995                */
  131.  
  132. #ifdef DEBUG
  133. main(argc, argv)
  134. int argc;
  135. int *argv;
  136. {
  137.   int i;
  138.   long j;
  139.   int n;
  140.   int m; 
  141.   long seed; 
  142.   long rngstate[RNGstatesize_2];
  143.  
  144.   m = atoi(argv[1]);    /* Number to discard initially */ 
  145.   n = atoi(argv[2]);    /* Number to print */ 
  146.   seed = atoi(argv[3]); /* Seed */ 
  147.  
  148.   sprand(seed,rngstate);
  149.  
  150.   if ( !ckrand(rngstate) ) printf("Error detected by ckrand()\n");
  151.  
  152.   for (i=0; i < m; i++) j = lprand(rngstate);
  153.   for (i=0; i < n; i++) printf("%ld\n", lprand(rngstate));
  154.  
  155.   if ( !ckrand(rngstate) ) printf("Error detected by ckrand()\n");
  156. } 
  157. #endif               
  158.  
  159.