/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "gentst.p" */ #include /* generate a bunch of numbers, each the sum of random numbers to test genhis. Tom Schneider */ /* end of program */ /* begin module version */ #define version 3.13 /* of gentst.p 1994 sep 5 origin 1985 dec 19 tds */ /* end module version */ /* begin module describe.gentst */ /* name gentst: test random generator synopsis gentst(gentstp:in, data: out, output: out) files gentstp: parameter file controlling the program. Three numbers, one per line: seed: random seed to start the process total: the number of numbers to generate components: the number of random numbers between 0 and 1 to add together to generate the total data: the input file for genhis. this is a set of numbers which should have gaussian distribution if the random number generator is a reasonable one. It will be N(0,1), a normal distribution with mean 0 and standard deviation 1. genhisp: control file for genhis output: messages to the user description test of a random number generator by creating a gaussian distribution of numbers for plotting by genhis example seed := 0.5; total := 10000; components:= 100; see also tstrnd.p, genhis.p author thomas d. schneider bugs none known technical notes the constant n in procedure randomtest determines how many times the random number generator will be in a series of tests. if n is small, the the test will be poor, if it is large then the test may take a long time. */ /* end module describe.gentst */ /* begin module gentst.var */ Static _TEXT gentstp, data, genhisp; /* output files */ #define pow14 16384 #define pow15 32768L #define pow22 4194304L #define pow23 8388608L /* end module gentst.var */ /* begin module random */ Static Void random(seed) double *seed; { /* random generator 2. version = 1.00; 1986 December 31 Test this routine with the program tstrnd. written by David Masternarde */ /* This random number generator is based on a shift register with a single bit of feedback, as described in Electronics for Neurobiologists, by Brown, Maxfield and Moraff, MIT press 1973, referencing Random Process Simulation and Measurement by Korn, McGraw-Hill 1966. The random seed rand, a number between 0 and 1 exclusive, is converted to an integer between 1 and 2**23-1, inclusive. This 23-bit number is shifted right one bit and the output of the last (23rd) bit and the 9th bit are added modulo 2 (exclusive orred) and fed back into the new first bit. This is done between 4 and 11 times, depending on the last 3 bits of the original number. The result is converted back to a real number between 0 and 1 from which the 23 bit integer can be recovered on the next call. The 23-bit shift register goes through all 2**23-1 values before repeating; the repetition frequency of this algorithm could be less or greater depending on the seed, because of the random number of multiple shifts per call. */ /* powers of 2 */ long iseed; /* integer shift register */ long i, nrep; /* index, number of times to do shift */ iseed = (long)floor(*seed * pow23 + 0.5); /* convert to 23 bit number */ if (iseed < 1 || iseed >= pow23) iseed = 1; nrep = (iseed & 7) + 4; /* do it 4 to 11 times based on last 3 bits */ for (i = 1; i <= nrep; i++) { /* if last bit and 9th bit are equal, feed back a 0, otherwise a 1 */ if ((iseed & 1) == ((iseed & (pow15 - 1)) >= pow14)) iseed /= 2; else iseed = iseed / 2 + pow22; } *seed = (double)iseed / pow23; } /* random */ #undef pow14 #undef pow15 #undef pow22 #undef pow23 #define xnumber 20 /* number of points on x axis to plot */ /* end module random */ /* begin module gentst.themain */ Static Void themain(data, genhisp) _TEXT *data, *genhisp; { /* the main procedure of gentst */ long total; /* total numbers to make */ long components; /* each number is made of this many random numbers */ long indexfortotal, indexforcomponents; /* indices */ double seed; /* the current random number seed */ double sum; /* the sum of the seeds */ printf("gentst %4.2f\n", version); if (*gentstp.name != '\0') { if (gentstp.f != NULL) gentstp.f = freopen(gentstp.name, "r", gentstp.f); else gentstp.f = fopen(gentstp.name, "r"); } else rewind(gentstp.f); if (gentstp.f == NULL) _EscIO2(FileNotFound, gentstp.name); RESETBUF(gentstp.f, Char); fscanf(gentstp.f, "%lg%*[^\n]", &seed); getc(gentstp.f); fscanf(gentstp.f, "%ld%*[^\n]", &total); getc(gentstp.f); fscanf(gentstp.f, "%ld%*[^\n]", &components); getc(gentstp.f); if (*genhisp->name != '\0') { if (genhisp->f != NULL) genhisp->f = freopen(genhisp->name, "w", genhisp->f); else genhisp->f = fopen(genhisp->name, "w"); } else { if (genhisp->f != NULL) rewind(genhisp->f); else genhisp->f = tmpfile(); } if (genhisp->f == NULL) _EscIO2(FileNotFound, genhisp->name); SETUPBUF(genhisp->f, Char); /* set x scale to fit a page: */ fprintf(genhisp->f, "x n %ld\n", (long)xnumber); fprintf(genhisp->f, "p g\n"); /* plot gaussian */ if (*data->name != '\0') { if (data->f != NULL) data->f = freopen(data->name, "w", data->f); else data->f = fopen(data->name, "w"); } else { if (data->f != NULL) rewind(data->f); else data->f = tmpfile(); } if (data->f == NULL) _EscIO2(FileNotFound, data->name); SETUPBUF(data->f, Char); fprintf(data->f, "* gentst %4.2f\n", version); fprintf(data->f, "* %10.9f random seed\n", seed); fprintf(data->f, "* %ld numbers, each the sum of\n", total); fprintf(data->f, "* %ld random numbers\n", components); for (indexfortotal = 1; indexfortotal <= total; indexfortotal++) { sum = 0.0; for (indexforcomponents = 1; indexforcomponents <= components; indexforcomponents++) { random(&seed); sum += seed; } fprintf(data->f, "% .5E\n", (sum - components / 2.0) * sqrt(12.0 / components)); } } #undef xnumber /* end module gentst.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); genhisp.f = NULL; strcpy(genhisp.name, "genhisp"); data.f = NULL; strcpy(data.name, "data"); gentstp.f = NULL; strcpy(gentstp.name, "gentstp"); themain(&data, &genhisp); _L1: if (gentstp.f != NULL) fclose(gentstp.f); if (data.f != NULL) fclose(data.f); if (genhisp.f != NULL) fclose(genhisp.f); exit(EXIT_SUCCESS); } /* End. */