1/*
2 * randbitrun - check rand bit run lengths of the a55 generator
3 *
4 * Copyright (C) 1999  Landon Curt Noll
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/13 03:43:11
21 * File existed as early as:	1995
22 *
23 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
24 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
25 */
26
27/*
28 * We will use randbit(1) to generate a stream if single bits.
29 * The odds that we will have n bits the same in a row is 1/2^n.
30 */
31
32
33define randbitrun(run_cnt)
34{
35    local i;			/* index */
36    local max_run;		/* longest run */
37    local long_run_cnt;		/* number of runs longer than MAX_RUN */
38    local run;			/* current run length */
39    local tally_sum;		/* sum of all tally values */
40    local last;			/* last random number */
41    local current;		/* current random number */
42    local MAX_RUN = 18;		/* max run we will keep track of */
43    local mat tally[1:MAX_RUN]; /* tally of length of a rise run of 'x' */
44    local mat prob[1:MAX_RUN];	/* prob[x] = probability of 'x' length run */
45
46    /*
47     * parse args
48     */
49    if (param(0) == 0) {
50	run_cnt = 65536;
51    }
52
53    /*
54     * run setup
55     */
56    max_run = 0;		/* no runs yet */
57    long_run_cnt = 0;		/* no long runs set */
58    current = randbit(1);	/* our first number */
59    run = 1;
60
61    /*
62     * compute the run length probabilities
63     *
64     * A bit run length of 'r' occurs with a probability of:
65     *
66     *			1/2^n;
67     */
68    for (i=1; i <= MAX_RUN; ++i) {
69	prob[i] = 1.0/(1<<i);
70    }
71
72    /*
73     * look at a number of random number trials
74     */
75    for (i=0; i < run_cnt; ++i) {
76
77	/* get our current number */
78	last = current;
79	current = randbit(1);
80
81	/* look for a run break */
82	if (current != last) {
83
84	    /*	record the stats */
85	    if (run > max_run) {
86		max_run = run;
87	    }
88	    if (run > MAX_RUN) {
89		++long_run_cnt;
90	    } else {
91		++tally[run];
92	    }
93
94	    /* start a new run */
95	    current = randbit(1);
96	    run = 1;
97
98	/* note the continuing run */
99	} else {
100	    ++run;
101	}
102    }
103    /* determine the number of runs found */
104    tally_sum = matsum(tally) + long_run_cnt;
105
106    /*
107     * print the stats
108     */
109    printf("rand runbit test used %d values to produce %d runs\n",
110      run_cnt, tally_sum);
111    for (i=1; i <= MAX_RUN; ++i) {
112	printf("length=%d\tprob=%9.7f\texpect=%d \tcount=%d \terr=%9.7f\n",
113	  i, prob[i], round(tally_sum*prob[i]), tally[i],
114	  (tally[i] - round(tally_sum*prob[i]))/tally_sum);
115    }
116    printf("length>%d\t\t\t\t\tcount=%d\n", MAX_RUN, long_run_cnt);
117    printf("max length=%d\n", max_run);
118}
119