1# THIS PURPOSELY DOES NOT HAVE A !# LINE !!!!
2#
3# Date: Mon, 9 Sep 2013 14:49:43 -0700
4# From: Bob Jewett <jewett@bill.scs.agilent.com>
5# Message-Id: <201309092149.r89Lnh94010909@bill.scs.agilent.com>
6# To: arnold@skeeve.com
7# Subject: Re: [bug-gawk] Bug in random() in builtin.c
8#
9# Hi Arnold,
10#
11# Attached below is a script that tests gawk for this particular
12# rand() problem.  The pair-wise combinations show a strong
13# autocorrelation for a delay of 31 pairs of rand() samples.
14#
15# The script prints out the measured autocorrelation for a record
16# of NSAMPLES pairs.  It also prints a fail message at the end if
17# it fails.
18#
19# If you want to see the autocorrelation values, there is a print
20# statement that if uncommented will save them to a file.
21#
22# Please let me know if the mailer screws up the transfer or
23# if you have any questions about the test.
24#
25# Best regards,
26# Bob
27#
28# -------------- test_pair_power_autocorrelation -----------------------
29#
30#!/bin/ksh
31
32#GAWK=/bin/gawk
33
34if [ -z "$AWK" ]; then
35    printf '$AWK must be set\n' >&2
36    exit 1
37fi
38
39# ADR: Get GAWK from the environment.
40# Additional note: This wants ksh/bash for the use of $RANDOM below to
41# seed the generator. However, shells that don't provide it won't be
42# a problem since gawk will then seed the generator with the time of day,
43# as srand() will be called without an argument.
44
45# large NSAMPLES and NRUNS will bring any correlation out of the noise better
46NSAMPLES=1024; MAX_ALLOWED_SIGMA=5; NRUNS=50;
47
48$AWK 'BEGIN{
49    srand('$RANDOM');
50    nsamples=('$NSAMPLES');
51    max_allowed_sigma=('$MAX_ALLOWED_SIGMA');
52    nruns=('$NRUNS');
53    for(tau=0;tau<nsamples/2;tau++) corr[tau]=0;
54
55    for(run=0;run<nruns;run++) {
56	sum=0;
57
58	# Fill an array with a sequence of samples that are a
59	# function of pairs of rand() values.
60
61	for(i=0;i<nsamples;i++) {
62	   samp[i]=((rand()-0.5)*(rand()-0.5))^2;
63	   sum=sum+samp[i];
64	   }
65
66	# Subtract off the mean of the sequence:
67
68	mean=sum/nsamples;
69	for(i=0;i<nsamples;i++) samp[i]=samp[i]-mean;
70
71	# Calculate an autocorrelation function on the sequence.
72	# Because the values of rand() should be independent, there
73	# should be no peaks in the autocorrelation.
74
75	for(tau=0;tau<nsamples/2;tau++) {
76	    sum=0;
77	    for(i=0;i<nsamples/2;i++) sum=sum+samp[i]*samp[i+tau];
78	    corr[tau]=corr[tau]+sum;
79	    }
80
81	}
82    # Normalize the autocorrelation to the tau=0 value.
83
84    max_corr=corr[0];
85    for(tau=0;tau<nsamples/2;tau++) corr[tau]=corr[tau]/max_corr;
86
87    # OPTIONALLY Print out the autocorrelation values:
88
89    # for(tau=0;tau<nsamples/2;tau++) print tau, corr[tau] > "pairpower_corr.data";
90
91    # Calculate the sigma for the non-zero tau values:
92
93    power_sum=0;
94
95    for(tau=1;tau<nsamples/2;tau++) power_sum=power_sum+(corr[tau])^2;
96
97    sigma=sqrt(power_sum/(nsamples/2-1));
98
99    # See if any of the correlations exceed a reasonable number of sigma:
100
101    passed=1;
102    for(tau=1;tau<nsamples/2;tau++) {
103	if ( abs(corr[tau])/sigma > max_allowed_sigma ) {
104	    print "Tau=", tau ", Autocorr=", corr[tau]/sigma, "sigma";
105	    passed=0;
106	    }
107        }
108    if(!passed) {
109	print "Test failed."
110	exit(1);
111        }
112    else exit (0);
113    }
114
115function abs(abs_input) { return(sqrt(abs_input^2)) ; }
116'
117