modified: makefile
[GalaxyCodeBases.git] / c_cpp / etc / calc / cal / randrun.cal
blobc9c54de9393f69316b03011703e461fbe1c603f1
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  * @(#) $Revision: 30.1 $
21  * @(#) $Id: randrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randrun.cal,v $
23  *
24  * Under source code control:   1995/02/12 20:00:06
25  * File existed as early as:    1995
26  *
27  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
28  */
31  * If X(j) < X(j+1) < ... X(j+k) >= X(j+k+1), then we have a run of 'k'.
32  * We ignore the run breaker, X(j+k+1), and start with X(j+k+2) when
33  * considering a new run in order to make our runs chi independent.
34  *
35  * See Knuth's "Art of Computer Programming - 2nd edition",
36  * Volume 2 ("Seminumerical Algorithms"), Section 3.3.2.
37  * "G. Run test", pp. 65-68,
38  * "problem #14", pp. 74, 536.
39  *
40  * We use the suggestion in problem #14 to allow an application of the
41  * chi-square test and to make estimating the run length probs easy.
42  */
45 define randrun(run_cnt)
47     local i;                    /* index */
48     local max_run;              /* longest run */
49     local long_run_cnt;         /* number of runs longer than MAX_RUN */
50     local run;                  /* current run length */
51     local tally_sum;            /* sum of all tally values */
52     local last;                 /* last random number */
53     local current;              /* current random number */
54     local MAX_RUN = 9;          /* max run we will keep track of */
55     local mat tally[1:MAX_RUN]; /* tally of length of a rise run of 'x' */
56     local mat prob[1:MAX_RUN];  /* prob[x] = probability of 'x' length run */
58     /*
59      * parse args
60      */
61     if (param(0) == 0) {
62         run_cnt = 65536;
63     }
65     /*
66      * run setup
67      */
68     max_run = 0;                /* no runs yet */
69     long_run_cnt = 0;           /* no long runs set */
70     current = rand();           /* our first number */
71     run = 1;
73     /*
74      * compute the run length probabilities
75      *
76      * A run length of 'r' occurs with a probability of:
77      *
78      *                  1/r! - 1/(r+1)!
79      */
80     for (i=1; i <= MAX_RUN; ++i) {
81         prob[i] = 1.0/fact(i) - 1.0/fact(i+1);
82     }
84     /*
85      * look at a number of random number trials
86      */
87     for (i=0; i < run_cnt; ++i) {
89         /* get our current number */
90         last = current;
91         current = rand();
93         /* look for a run break */
94         if (current < last) {
96             /*  record the stats */
97             if (run > max_run) {
98                 max_run = run;
99             }
100             if (run > MAX_RUN) {
101                 ++long_run_cnt;
102             } else {
103                 ++tally[run];
104             }
106             /* start a new run */
107             current = rand();
108             run = 1;
110         /* note the continuing run */
111         } else {
112             ++run;
113         }
114     }
115     /* determine the number of runs found */
116     tally_sum = matsum(tally) + long_run_cnt;
118     /*
119      * print the stats
120      */
121     printf("rand run test used %d values to produce %d runs\n",
122       run_cnt, tally_sum);
123     for (i=1; i <= MAX_RUN; ++i) {
124         printf("length=%d\tprob=%9.7f\texpect=%d \tcount=%d \terr=%9.7f\n",
125           i, prob[i], round(tally_sum*prob[i]), tally[i],
126           (tally[i] - round(tally_sum*prob[i]))/tally_sum);
127     }
128     printf("length>%d\t\t\t\t\tcount=%d\n", MAX_RUN, long_run_cnt);
129     printf("max length=%d\n", max_run);
132 if (config("resource_debug") & 3) {
133     print "randrun([run_length]) defined";