home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume25 / mrandom / mrandom.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-12  |  10.3 KB  |  317 lines
  1. /*
  2.  *        mrandom.c
  3.  *
  4.  *    Wrapper for random(), allowing safe restart.
  5.  *
  6.  *    Original Implementation:
  7.  *        Clark Thomborson, September 1991
  8.  *
  9.  *    This material is based upon work supported by the National
  10.  *    Science Foundation under grant number MIP-9023238.  The
  11.  *    Government has certain rights in this material.
  12.  *
  13.  *    Any opinions, findings, and conclusions or recommendations
  14.  *    expressed in this material are those of the author and do
  15.  *    not necessarily reflect the view of the National Science
  16.  *    Foundation.
  17.  *
  18.  *    This code is neither copyrighted nor patented.
  19.  */
  20.  
  21. #include <values.h> /* we need MAXLONG */
  22. #include <stdio.h>  /* we need FILE */
  23. #include "mrandom.h"
  24.  
  25. #define RNGSTATELENGTH    32 /* legal values are 2, 8, 16, 32, 64 */
  26. #define RNGSTATE0    3  /* rngstate[0] for table length 32 */
  27. #define BILLION    1000000000 /* convenient modulus for a 32-bit counter */
  28.  
  29. /* Data in our RNG statefile */
  30. long    RNGseed;    /* the seed originally used to initialize rngstate */
  31. long    RNGcount1;  /* mod-BILLION counter of random() calls */
  32. long    RNGcount2;  /* div-BILLION counter */
  33. long    rngstate[RNGSTATELENGTH]; /* an array passed to random() */
  34. long    RNGnextval; /* random()'s next output */
  35.  
  36. /* Format of our RNG statefile */
  37. #define RNGfileLINE1    "Initial seed = %d\n"
  38. #define RNGfileLINE2    \
  39.     "Number of mrandom() calls after seeding = %ld billion + %ld\n"
  40. #define RNGfileLINE3    "RNG state table =\n"
  41. #define RNGfileLINEn    "   %08lx %08lx %08lx %08lx\n"
  42. #define RNGfileLINE12    "Next value in this pseudorandom sequence = %08lx\n"
  43.  
  44. /* Write a RNG state identifier into the user-supplied string rngid,
  45.  * which must be of length at least RNGIDSTRLEN.  If the user has not
  46.  * initialized the rng with init_rng(), restart_rng(), or reconstruct_rng(),
  47.  * abort with an error message to stderr.  Otherwise return the value of rngid.
  48.  */
  49. extern char *describe_rng (rngid)
  50. char rngid[RNGIDSTRLEN];
  51. {
  52.   if (rngstate[0] != RNGSTATE0) {
  53.     fprintf(stderr, "RNG has not been initialized!\n");
  54.     fflush(stderr);
  55.     exit(1);
  56.   }
  57.   sprintf(rngid, "RNG state identifier is (%ld, %ld, %ld)\n",
  58.          RNGseed, RNGcount1, RNGcount2);
  59.   return(rngid);
  60. }
  61.  
  62. /* Create a random number statefile initialized with the given seed.
  63.  * Return 1 if file is successfully created, 0 otherwise.
  64.  */
  65. extern int init_rng (seed, filename)
  66. int seed;
  67. char *filename;
  68. {
  69.   reconstruct_rng(seed,0,0);
  70.   return(save_rng(filename));
  71. } /* end init_rng */
  72.  
  73. /* Rebuild a random() state by reseeding the generator, then making
  74.  * (count2*1e9 + count1) calls to mrandom().  Useful for error-checking
  75.  * and error-recovery routines, although it is very slow for large counts.
  76.  */
  77. extern void reconstruct_rng (seed, count1, count2)
  78. long seed, count1, count2;
  79. {
  80.   double a;
  81.   initstate(seed, rngstate, RNGSTATELENGTH*4);
  82.   RNGseed = seed; /* keep a record of the seed */
  83.   RNGcount1 = 0;  /* mod-billion counter */
  84.   RNGcount2 = 0;  /* div-billion counter */
  85.   /* watch out for infinite loops */
  86.   if (count1 >= 0 && count1 < BILLION && count2 >= 0) {
  87.     if (count2 != 0) {
  88.       fprintf(stderr, "Warning: this reconstruction will take a LONG time!\n");
  89.       fflush(stderr);
  90.     }
  91.     for ( ; (RNGcount1 != count1) || (RNGcount2 != count2) ; ) {
  92.       a = frandom();
  93.     }
  94.   }
  95. } /* end reconstruct_rng */
  96.  
  97. /* Return the next random() output without disturbing the RNG state.
  98.    This procedure will only work properly if random's state is in rngstate[].
  99.  */
  100. int nextval()
  101. {
  102.   long state[RNGSTATELENGTH];
  103.   long i, r, retval;
  104.   for (i=0; i<RNGSTATELENGTH; i++) {
  105.     state[i] = rngstate[i];
  106.   }
  107.   retval = random(); /* next value in this pseudorandom sequence */
  108.   for (i=1; i<RNGSTATELENGTH-1; i++) {
  109.     r = random(); /* these calls are necessary to restore random()'s state */
  110.   }
  111.   for (i=0; i<RNGSTATELENGTH; i++) {
  112.     rngstate[i] = state[i];
  113.   }
  114.   return(retval);
  115. }
  116.  
  117. /* Restart a generator from a statefile.  Print a message on stderr
  118.  * if the restart failed due to a garbled or non-existent statefile, 
  119.  * and return 0.  Otherwise return 1.
  120.  */
  121. extern int restart_rng(filename)
  122. char *filename;
  123. {
  124.   FILE *fp;
  125.   long i, m, newstate[RNGSTATELENGTH], r; /* temps */
  126.   int errflag; /* initially 0, becomes 1 if a problem is discovered */
  127.  
  128.   /* restore random()'s internal variables to their initial state */
  129.   m = (RNGcount1%(RNGSTATELENGTH-1)) + RNGcount2*(BILLION%(RNGSTATELENGTH-1));
  130.   for ( ; m%(RNGSTATELENGTH-1) != 0; m++) {
  131.     r = random();
  132.   }
  133.  
  134.   /* restore counter values, retrieve original RNG seed and current state */
  135.   fp = fopen(filename, "r");
  136.   if (!fp) {
  137.     fprintf(stderr, "There is no RNG statefile in this directory!\n");
  138.     fflush(stderr);
  139.     exit(1);
  140.   }
  141.   fscanf(fp, RNGfileLINE1, &RNGseed);
  142.   fscanf(fp, RNGfileLINE2, &RNGcount2, &RNGcount1);
  143.   fscanf(fp, RNGfileLINE3);
  144.   for (i=0; i<RNGSTATELENGTH; i++) {
  145.     fscanf(fp, "%lx", &newstate[i]);
  146.   }
  147.   fscanf(fp, "\n");
  148.   fscanf(fp, RNGfileLINE12, &RNGnextval);
  149.   fclose(fp);
  150.  
  151.   /* error check: the first word of a 32-word rngstate is always 3 */
  152.   if (newstate[0] != RNGSTATE0) {
  153.     errflag = 1;
  154.   } else {
  155.     errflag = 0;
  156.   }
  157.  
  158.   /* tell random() we have a 32-word rngstate */
  159.   rngstate[0] = RNGSTATE0;
  160.   setstate(rngstate);
  161.  
  162.   /* If reconstruction will be rapid, do it as an additional error check. */ 
  163.   if (RNGcount1 < 50 && RNGcount2 == 0) {
  164.     reconstruct_rng(RNGseed, RNGcount1, RNGcount2);
  165.     /* see if we got to the same state */
  166.     for (i=0; i<RNGSTATELENGTH; i++) {
  167.       if (newstate[i] != rngstate[i]) {
  168.         errflag = 1;
  169.       }
  170.     }
  171.   } else {
  172.     /* quickly modify random()'s internal state to line up with RNGcount info */
  173.     m = (RNGcount1%(RNGSTATELENGTH-1)) + RNGcount2*(BILLION%(RNGSTATELENGTH-1));
  174.     for (i=0; i < m%(RNGSTATELENGTH-1); i++) {
  175.       r = random();
  176.     }
  177.   }
  178.  
  179.   /* copy in the new state */
  180.   for (i=0; i<RNGSTATELENGTH; i++) {
  181.     rngstate[i] = newstate[i];
  182.   }
  183.  
  184.   /* Check nextval() operation, and verify RNGnextval */
  185.   if (nextval() != nextval() || RNGnextval != nextval()) {
  186.     errflag = 1;
  187.   }
  188.  
  189.   if (errflag) {
  190.     fprintf(stderr,
  191.     "Warning: RNG statefile is inconsistent.  Did you edit it?\n");
  192.     fprintf(stderr,
  193.     "If not, check your program to make sure you:\n");
  194.     fprintf(stderr,
  195.     "   1. use frandom() or mrandom(m), not random();\n");
  196.     fprintf(stderr,
  197.     "   2. use init_rng(seed, filename), not srandom(seed);\n");
  198.     fprintf(stderr,
  199.     "   3. use restart_rng(filename), not setstate(state); and\n");
  200.     fprintf(stderr,
  201.     "   4. don't overwrite the private storage of mrandom.\n");
  202.     fflush(stderr);
  203.     return(0);
  204.   } else {
  205.     return(1); /* everything looks ok */
  206.   }
  207. } /* end restart_rng */
  208.  
  209. /* Generate a uniformly-distributed number a, 0.0 <= a < 1.0, using
  210.  * the 4.3bsd additive congruential generator random().
  211.  *
  212.  * This routine keeps track of the number of calls to random().
  213.  *
  214.  * From inspection of the object code for random() in a Sun-3 release,
  215.  * I deduce the following.  A call to initstate() or setstate()
  216.  *   1. defines the location of the state array (rngstate in the code below);
  217.  *   2. defines the length of the state array (32 words in the code below;
  218.  *    other possibilities are 2, 8, 16, or 64 words);
  219.  *   3. initializes several static variables to point at the state array
  220.  *    (the indices j and k in the code below);
  221.  *   4. defines the randomization algorithm (a linear congruential
  222.  *    generator is used if the state array has length 2).
  223.  * Additionally, initstate(seed) fills the state array with the output
  224.  * of a linear congruential generator that is seeded with the given seed.
  225.  *
  226.  * Here is a disassembled and decompiled listing of random(), as
  227.  * it would operate when initialized with my init_rng routine:
  228.  *
  229.  * long random()
  230.  * {
  231.  *   static int j=4, k=1;
  232.  *   long r;   
  233.  *   rngstate[j] += rngstate[k];
  234.  *   r = (rngstate[j]>>1) & 0x7fffffff;
  235.  *   j++;
  236.  *   if (j>31) j=1;
  237.  *   k++;
  238.  *   if (k>31) k=1;
  239.  *   return(r);
  240.  * }
  241.  *
  242.  * Note that j and k return to their original values after every 31
  243.  * calls to random().  This property allows me to write a restart_rng()
  244.  * routine even though I can't read the values of j and k. 
  245.  *
  246.  * Also note that this is a very poor RNG algorithm if it is only
  247.  * called once between successive setstate() calls.  In this pathological
  248.  * case, each output of random() will differ from the prior random() output
  249.  * by rngstate[1]/2 or rngstate[1]/2+1, mod 2^31-1.  Even in non-pathological
  250.  * calling sequences the long-period property of an additive congruential
  251.  * generator is only guaranteed if j and k always move through the array
  252.  * in the manner shown.  Thus the need for counting random() calls
  253.  * (or for writing your own time-optimized additive congruential generator).
  254.  *
  255.  * Finally, note that integer overflow is a frequent occurrence in
  256.  * rngstate[k] += rngstate[j].  Thus random() could conceivably
  257.  * behave differently on different machines.
  258.  */
  259. extern double frandom ()
  260. {
  261.   RNGcount1++; /* mod-billion counter */
  262.   if (RNGcount1 == BILLION) {
  263.     RNGcount1 = 0;
  264.     RNGcount2++; /* an overflow is unlikely in my lifetime */
  265.   }
  266.   return( (double) random()/MAXLONG );
  267. } /* end frandom */
  268.  
  269. /* Generate a random integer, uniformly distributed in the range 0..m-1.
  270.  * We use the most-significant bits of the result of a random() call,
  271.  * although this may only marginally improve the quality of our output.
  272.  * You may wish to rewrite this routine if you don't have hardware
  273.  * floating point.  If you do, be sure to include the counter-bumping
  274.  * code from frandom(). 
  275.  */
  276. extern long mrandom (m)
  277. long m;
  278. {
  279.   return( (int) (frandom()*m) );
  280. } /* end mrandom */
  281.  
  282. /* Save the RNG state to a statefile, after calling random() enough
  283.  * times to reset its internal state variables to their initial values.
  284.  * Check to be sure the RNG can be restarted by calling restart_rng().
  285.  * Return 0 if a problem is detected, printing an error message on stderr.
  286.  * Otherwise return 1.
  287.  */
  288. extern int save_rng(filename)
  289. char *filename;
  290. {
  291.   FILE *fp;
  292.   long i;
  293.  
  294.   /* write the statefile */
  295.   fp = fopen(filename, "w");
  296.   if (!fp) {
  297.     fprintf(stderr, "Trouble opening RNG statefile %s for writing.", filename);
  298.     fflush(stderr);
  299.     return(0);
  300.   }
  301.   fprintf(fp, RNGfileLINE1, RNGseed);
  302.   fprintf(fp, RNGfileLINE2, RNGcount2, RNGcount1);
  303.   fprintf(fp, RNGfileLINE3);
  304.   for (i=0; i<RNGSTATELENGTH; i+=4) {
  305.     fprintf(fp, RNGfileLINEn,
  306.         rngstate[i], rngstate[i+1], rngstate[i+2], rngstate[i+3]);
  307.   }
  308.   fprintf(fp, RNGfileLINE12, nextval());
  309.   fclose(fp);
  310.  
  311.   /* Return after checking that state was saved correctly. */
  312.   return( restart_rng(filename) );
  313.  
  314. } /* end save_rng */
  315.  
  316. /* end mrandom.c */
  317.