home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #16 / NN_1992_16.iso / spool / sci / math / stat / 1513 < prev    next >
Encoding:
Text File  |  1992-07-25  |  2.2 KB  |  76 lines
  1. Path: sparky!uunet!zaphod.mps.ohio-state.edu!usc!news
  2. From: ajayshah@almaak.usc.edu (Ajay Shah)
  3. Newsgroups: sci.math.stat
  4. Subject: C source for Park-Miller "minimal std. generator":
  5. Date: 25 Jul 1992 14:05:23 -0700
  6. Organization: University of Southern California, Los Angeles, CA
  7. Lines: 62
  8. Sender: ajayshah@almaak.usc.edu (Ajay Shah)
  9. Message-ID: <l73ggjINNh52@almaak.usc.edu>
  10. References: <1992Jul25.191638.16109@samba.oit.unc.edu>
  11. NNTP-Posting-Host: almaak.usc.edu
  12. Keywords: pseudo-random, C , uniform, normal
  13.  
  14. #include <stdio.h>
  15.  
  16. /* This file implements the "minimal standard" RNG
  17.    from the paper "RNGs: Good Ones are Hard to Find" by Park and
  18.    Miller, CACM, Volume 31, Number 10, October 1988. */
  19.  
  20. static long int seed = 1;
  21. /* static so programs will be more "orderly" in modifying the
  22.    seed.  The setseed function is the prescribed way in which the seed
  23.    is altered.  A good default to use is your SSN.
  24.    Any seed from 1 to 2147483646 is equally ok. */
  25.  
  26. int setseed(long int newseed)
  27. {
  28.     if ((newseed < 1) || (newseed > 2147483646)) return 1;
  29.     seed = newseed; return 0;
  30. }
  31.  
  32. double ranpm()
  33.      /* page 1195-right, "Integer Version 2", works iff ints are 32 bits.*/
  34. {
  35.     long int a=16807, m=2147483647, q=127773, r=2836;
  36.     int lo, hi, test;
  37.  
  38.     hi = seed / q;
  39.     lo = seed % q;
  40.     test = a * lo - r * hi;
  41.     seed = (test > 0) ? test : test + m;
  42.     return (seed/ (double) m);
  43. }
  44.  
  45. #ifdef TESTING
  46. /* They give us one test (on page 1195-left): if seed = 1 then
  47.    the 10001'th random number should be 1043618065.
  48.    Apart from this, I do some superficial testing aimed at checking
  49.    the code above is not buggy.  I'm not trying to test the
  50.    statistical properties of the RNG. */
  51.  
  52. int main()
  53. {
  54.     int n; double u, sum=0.0, sum2=0.0;
  55.     
  56.     for (n=1; n<=10000; n++) {
  57.     u = ranpm();
  58.     if ((u < 0.0) || (u > 1.0)) {
  59.         printf("Garbage at the %dth number!\n");
  60.         return 1;
  61.     }
  62.     sum += u;
  63.     sum2 += (u*u);
  64.     }
  65.     printf("The current value of seed is : %10d\n", seed);
  66.     printf("The Truth is %17s 1043618065\n\n", " ");
  67.  
  68.     printf("Mean of 1st 10k nums is %f, variance is %f\n",
  69.        sum/10000, sum2/9999);
  70.     printf("Expect .5018 and .3355\n");
  71.     return 0;
  72. }
  73. #endif /* TESTING */
  74. -- 
  75. Ajay Shah, (213)749-8133, ajayshah@usc.edu
  76.