c_cpp/etc/calc/cal/randrun.cal

   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  */
  29
  30 /*
  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  */
  43
  44
  45 define randrun(run_cnt)
  46 {
  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 */
  57
  58     /*
  59      * parse args
  60      */
  61     if (param(0) == 0) {
  62         run_cnt = 65536;
  63     }
  64
  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;
  72
  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     }
  83
  84     /*
  85      * look at a number of random number trials
  86      */
  87     for (i=0; i < run_cnt; ++i) {
  88
  89         /* get our current number */
  90         last = current;
  91         current = rand();
  92
  93         /* look for a run break */
  94         if (current < last) {
  95
  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             }
 105
 106             /* start a new run */
 107             current = rand();
 108             run = 1;
 109
 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;
 117
 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);
 130 }
 131
 132 if (config("resource_debug") & 3) {
 133     print "randrun([run_length]) defined";
 134 }