home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e065 / 3.ddi / N2RAND.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-28  |  10.1 KB  |  325 lines
  1. /* n2rand.c - NeuralWare Standard Random Number Generator */
  2. /* Version 1.1 of 6 March 1990 */
  3.  
  4. /************************************************************************
  5.  *                  *
  6.  *  Random Number Generator for NeuralWare Products     *
  7.  *                  *
  8.  ************************************************************************
  9.     05-Dec-90 cck
  10.   Modified n2rand.c and n2rand.h for stand-alone operation and
  11.     easy integration with flash code.
  12.   Moved RRand() and GRand() from n2fl_sum.c to n2rand.c
  13.  */
  14.  
  15. /* Its purpose is to provide a standard random number generator across all  */
  16. /* products that NeuralWare sells and supports.                             */
  17.  
  18. /* This code derived from Knuth, The Art of Computer Programming, Vol2,     */
  19. /* Seminumerical Algorithms; Sedgewick, Algorithms; and Press, Flannery,    */
  20. /* Teukolsky, and Vetterling, Numerical Recipes in C                        */
  21.  
  22. /* The concept is to use one type of generator, a linear congruential       */
  23. /* generator, to fill an array with random numbers. We then shuffle that    */
  24. /* array, and then take the difference between two of those random numbers  */
  25. /* in the array. That value is returned to the user, and also used to       */
  26. /* replace one of the numbers in the random number array.                   */
  27.  
  28. /* NRAND Returns a uniform random deviate between 0 and MAXRAND-1.          */
  29. /* Call NRAND as so: (SL) new_deviate = nrand();                            */
  30.  
  31. /* SRAND Sets seed of random sequence to another part of the sequence.      */
  32. /* Call SRAND as so: SRAND( (SL) mi_random_seed );                          */
  33.  
  34. #include <math.h>
  35.  
  36. #include "n2rand.h" /* defines MAXRAND & others       */
  37.  
  38.   static long    rand_array[56];           /* Lookup table of random numbers */
  39.   static long    first_gen, second_gen;    /* Result for the LCG             */
  40.   static int     table_index, loop;        /* Index into random number table */
  41.   static int     first_time_flag = 1;      /* First time thru flag -
  42.                    initially true */
  43.   static int     cursor_1, cursor_2, i;    /* Cursors to chose random #'s from
  44.                 lookup table of random #'s    */
  45. /*   */
  46. /************************************************************************
  47.  *                  *
  48.  *  set_rand() - Random number seed         *
  49.  *                  *
  50.  ************************************************************************
  51.  */
  52.  
  53. #ifdef PROTOTYPING
  54. void set_rand( long seed )
  55. #else
  56. void set_rand( seed )
  57. long seed;
  58. #endif
  59. {
  60. /* Reset first_time_flag on the first call */
  61.   if ( first_time_flag ) {
  62.   /* Now reset first time flag */
  63.       first_time_flag = 0;                /* now false */
  64.       }
  65.  
  66. /* Now test & correct seed value */
  67.       first_gen = MSEED - (seed < 0 ? -seed : seed);
  68.       first_gen %= MAXRAND;
  69.       rand_array[55] = first_gen;
  70.       second_gen = 1;
  71.  
  72. /* Now initialize the shuffle table */
  73.    /* We use a linear congruential generator to set up initial table */
  74.       for ( table_index = 1; table_index <= 54; table_index++)
  75.     {
  76.     i = (21*table_index) % 55;
  77.     rand_array[i] = second_gen;
  78.     second_gen = first_gen - second_gen;
  79.     if ( second_gen < MZ ) second_gen += MAXRAND;
  80.     first_gen = rand_array[i];
  81.     } /* End of for table_index = 1 to 54 */
  82.  
  83.    /* Now randomize the table */
  84.       for (loop=1; loop<=4; loop++)
  85.       {
  86.     for ( table_index=1; table_index <= 55; table_index++)
  87.         {
  88.         rand_array[table_index] -= rand_array[1 +
  89.              (table_index+30) % 55];
  90.         if ( rand_array[table_index] < MZ ) rand_array[table_index]
  91.              += MAXRAND;
  92.         }
  93.       }
  94.  
  95.    /* Finish up any remaining bits */
  96.  
  97.       cursor_1 = 0;                       /* init cursors into random # table */
  98.       cursor_2 = 31;
  99.  
  100.  
  101. }   /* End of initialization/setup code (if first_time true or reinit) */
  102. /*   */
  103. /************************************************************************
  104.  *                  *
  105.  *  nrand() - NeuralWare Random Number Generator      *
  106.  *                  *
  107.  ************************************************************************
  108.  */
  109.  
  110. long nrand()
  111. {
  112. /* Retrieve random number & generate a new replacement value */
  113. /* Initialize on the first call, if not already initialized by set_rand */
  114.   if ( first_time_flag ) {
  115.       set_rand(0);
  116.   /* Now reset first time flag */
  117.       first_time_flag = 0;                /* now false */
  118.       }
  119.  
  120.   if (++cursor_1 == 56) cursor_1 = 1;     /* rollover so 1 to 55 */
  121.   if (++cursor_2 == 56) cursor_2 = 1;     /* rollover so 1 to 55 */
  122.  
  123.             /* now get diff of 2 random #'s */
  124.   first_gen = rand_array[cursor_1]-rand_array[cursor_2];
  125.  
  126.             /* correct for rollover */
  127.   if ( first_gen < MZ ) first_gen += MAXRAND;
  128.  
  129.             /* replace old random # */
  130.   rand_array[cursor_1] = first_gen;
  131.  
  132.             /* and return random # to user */
  133.   return first_gen;
  134. }
  135. /*   */
  136. /************************************************************************
  137.  *                  *
  138.  *  RRand(), GRand - real random number       *
  139.  *                  *
  140.  *  RRand - Uniform Distribution          *
  141.  *  GRand - Gaussian Distribution         *
  142.  *                  *
  143.  ************************************************************************
  144.  */
  145.  
  146. double RRand( lv, hv )
  147. double lv, hv;    /* low / high limits */
  148. {
  149.     double rv;    /* return value */
  150.  
  151.     rv = nrand() * RANDINV;   /* rand() / whatever value is needed */
  152.  
  153.     rv = (hv - lv) * rv + lv;   /* range adjusted */
  154.  
  155.     return( rv );
  156. }
  157.  
  158. double GRand( lv, hv )
  159. double lv, hv;    /* low / high limits */
  160. {
  161.     static int    k=0;      /* k=0 need to caluclate g0 & g1,
  162.            * k=1 g1 already calculated */
  163.     static double g0,g1;    /* Gaussian(0,1) random variables */
  164.     double u0,u1;     /* uniforms on -1 to 1 */
  165.     double r2;        /* radius from (0,0) to (u0,u1) */
  166.     double x,y;       /* intermediate calculations */
  167.  
  168.     if( k ) {
  169.       k = 0;
  170.       return( 0.5 * (hv - lv) * g1 + 0.5 * (hv + lv) );
  171.     } else {
  172.       k = 1;
  173.       /* find two U(-1,1)'s that are independent and lie within unit circle
  174.        * using NeuralWare's stanard uniform generator
  175.        */
  176.       u1 = 2.0 * (nrand() * RANDINV) - 1.0;
  177.       y = u1 * u1;
  178.       for(r2=2.0; r2>1.0; ) {
  179.         u0 = u1;
  180.         r2 = y;
  181.         u1 = 2.0 * (nrand() * RANDINV) - 1.0;
  182.         y = u1 * u1;
  183.         r2 += y;
  184.       }
  185.       /* now get two Gaussian(0,1) */
  186.       x = sqrt( ( -2.0 * log(r2) )/ r2 );
  187.       g0 = u0 * x;
  188.       g1 = u1 * x;
  189.       return( 0.5 * (hv - lv) * g0 + 0.5 * (hv + lv) );
  190.     }
  191. }
  192. /*   */
  193. /************************************************************************
  194.  *                  *
  195.  *  Test Routine for Random Number Generator      *
  196.  *                  *
  197.  ************************************************************************
  198.  */
  199. #ifdef TEST_RAN
  200.  
  201. /* The following code is a simple test harness for the random number        */
  202. /* generator. This was derived from Sedgewick, Algorithms, Addison-Wesley,  */
  203. /* 1984.                                                                    */
  204.  
  205. #include "usermath.h"
  206.  
  207. #define R_MAX 2000           /* This implies up to 2000 values in histogram */
  208.  
  209. SREAL chi_square( N, r )
  210. UL    N, r;
  211.  
  212. /*
  213.   This routine calculates the sum of the squares of the frequencies
  214.   of occurence of each random index produced via nrand. The sum is
  215.   then scaled by the expected frequency of occurence, and normalized
  216.   by the size of the sequence.
  217. */
  218.  
  219. /*
  220.   N is the total number of random numbers we wish to generate
  221.   r is the range we wish to produce random number in; The numbers
  222.   generated range from 0 thru r-1.
  223.  
  224.   A good fit is indicated by how close this 'Chi-Square' is to the
  225.   expected value, i.e. for 10000 random numbers ranging between 1 &
  226.   100, a 'good' Chi-square would be around 100 in value
  227. */
  228.  
  229. {
  230. static long f[R_MAX];
  231. static long i, index;
  232. static SREAL temp, t1, t2;
  233.  
  234. /* First clear out the histogram array */
  235. for (index=0; index<R_MAX; index++) f[index] = 0;
  236.  
  237. /* Now accumulate a histogram of times a given index is produced */
  238.   for (i=0; i<N; i++) {
  239.  
  240. #if 1                      /* set this 1 to test Neuralware generator */
  241.  
  242.   t1 = nrand();
  243.   t2 = r;
  244.   index = (t1*t2)*RANDINV;
  245. #endif
  246.  
  247. #if 0                      /* set this 1 to test alternate generator */
  248.   t1 = random();
  249.   t2 = r;
  250.   index = (t1*t2)/2147483647.0;
  251. #endif
  252.  
  253.   if ( index >= R_MAX || index < 0 ) {
  254.   printf("ERROR - Chi-square, index = %d\n",i);
  255.   index = 0;
  256.   };
  257.   f[index] = f[index] + 1;
  258.   }
  259.  
  260. /* Now calculate the Chi Square of those frequencies */
  261.   /* We use floating point calculations to see any patterns in low order bits */
  262.   t1 = 0.0;
  263.   t2 = r;
  264.  
  265.   for (i=0; i<=r-1 ; i++)
  266.      t1 = t1 + (f[i]-(N/r))*(f[i]-(N/r));
  267.  
  268.   temp = ((t2*t1)/N);
  269.  
  270. return temp;
  271.  
  272. }
  273.  
  274. #ifdef __ZTC__                            /* if Zortech C */
  275. #include <time.h>
  276. #endif
  277.  
  278. main()
  279. #define SAMPLE_SIZE 1000
  280. #define RANGE       100
  281. #define NTRIALS     10
  282. {
  283. NINT     i;            /* loop index */
  284. SREAL temp[NTRIALS];   /* temporary store */
  285. SL    seed;            /* seed for random number generator */
  286.  
  287. #ifdef __ZTC__                            /* if Zortech C */
  288. time_t start, finish;  /* data structures to measure run time of rand funct */
  289. #endif
  290.  
  291. /* First reseed random number generator */
  292.   seed = 0;
  293.   set_rand(seed);
  294.  
  295. /* Second, start the loop for Ntrials of random numbers */
  296.  
  297. #ifdef __ZTC__                            /* if Zortech C */
  298.   time (&start);
  299.   for (i=0; i< NTRIALS; i++) {
  300.      temp[i] = chi_square(SAMPLE_SIZE, RANGE);
  301.      }
  302.   time(&finish);
  303. #else
  304.   for (i=0; i< NTRIALS; i++) {
  305.      temp[i] = chi_square(SAMPLE_SIZE, RANGE);
  306.      }
  307. #endif
  308.  
  309. /* Third, print out the duration per call & CHI square values */
  310.  
  311. #ifdef __ZTC__                            /* if Zortech C */
  312.   printf("time per Random call %2.10f seconds\n",
  313.            (difftime(finish, start)/SAMPLE_SIZE));
  314. #endif
  315.  
  316.   for (i=0; i< NTRIALS; i++)
  317.      printf("Trial # %d    CHI Square = %f\n",i,temp[i]);
  318.  
  319.  
  320.  
  321. }
  322. #endif /* ifdef TEST_RAN */
  323.  
  324.  
  325.