1/*
2 * randrun - perform a run test on rand()
3 *
4 * Copyright (C) 1999  David I. Bell
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1995/02/12 20:00:06
21 * File existed as early as:	1995
22 *
23 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26/*
27 * If X(j) < X(j+1) < ... X(j+k) >= X(j+k+1), then we have a run of 'k'.
28 * We ignore the run breaker, X(j+k+1), and start with X(j+k+2) when
29 * considering a new run in order to make our runs chi independent.
30 *
31 * See Knuth's "Art of Computer Programming - 2nd edition",
32 * Volume 2 ("Seminumerical Algorithms"), Section 3.3.2.
33 * "G. Run test", pp. 65-68,
34 * "problem #14", pp. 74, 536.
35 *
36 * We use the suggestion in problem #14 to allow an application of the
37 * chi-square test and to make estimating the run length probs easy.
38 */
39
40
41define randrun(run_cnt)
42{
43    local i;			/* index */
44    local max_run;		/* longest run */
45    local long_run_cnt;		/* number of runs longer than MAX_RUN */
46    local run;			/* current run length */
47    local tally_sum;		/* sum of all tally values */
48    local last;			/* last random number */
49    local current;		/* current random number */
50    local MAX_RUN = 9;		/* max run we will keep track of */
51    local mat tally[1:MAX_RUN]; /* tally of length of a rise run of 'x' */
52    local mat prob[1:MAX_RUN];	/* prob[x] = probability of 'x' length run */
53
54    /*
55     * parse args
56     */
57    if (param(0) == 0) {
58	run_cnt = 65536;
59    }
60
61    /*
62     * run setup
63     */
64    max_run = 0;		/* no runs yet */
65    long_run_cnt = 0;		/* no long runs set */
66    current = rand();		/* our first number */
67    run = 1;
68
69    /*
70     * compute the run length probabilities
71     *
72     * A run length of 'r' occurs with a probability of:
73     *
74     *			1/r! - 1/(r+1)!
75     */
76    for (i=1; i <= MAX_RUN; ++i) {
77	prob[i] = 1.0/fact(i) - 1.0/fact(i+1);
78    }
79
80    /*
81     * look at a number of random number trials
82     */
83    for (i=0; i < run_cnt; ++i) {
84
85	/* get our current number */
86	last = current;
87	current = rand();
88
89	/* look for a run break */
90	if (current < last) {
91
92	    /*	record the stats */
93	    if (run > max_run) {
94		max_run = run;
95	    }
96	    if (run > MAX_RUN) {
97		++long_run_cnt;
98	    } else {
99		++tally[run];
100	    }
101
102	    /* start a new run */
103	    current = rand();
104	    run = 1;
105
106	/* note the continuing run */
107	} else {
108	    ++run;
109	}
110    }
111    /* determine the number of runs found */
112    tally_sum = matsum(tally) + long_run_cnt;
113
114    /*
115     * print the stats
116     */
117    printf("rand run test used %d values to produce %d runs\n",
118      run_cnt, tally_sum);
119    for (i=1; i <= MAX_RUN; ++i) {
120	printf("length=%d\tprob=%9.7f\texpect=%d \tcount=%d \terr=%9.7f\n",
121	  i, prob[i], round(tally_sum*prob[i]), tally[i],
122	  (tally[i] - round(tally_sum*prob[i]))/tally_sum);
123    }
124    printf("length>%d\t\t\t\t\tcount=%d\n", MAX_RUN, long_run_cnt);
125    printf("max length=%d\n", max_run);
126}
127
128if (config("resource_debug") & 3) {
129    print "randrun([run_length]) defined";
130}
131