1
# THIS PURPOSELY DOES NOT HAVE A !# LINE !!!!
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
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.
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
19
# If you want to see the autocorrelation values, there is a print
20
# statement that if uncommented will save them to a file.
22
# Please let me know if the mailer screws up the transfer or
23
# if you have any questions about the test.
28
# -------------- test_pair_power_autocorrelation -----------------------
34
# ADR: Get GAWK from the environment.
35
# Additional note: This wants ksh/bash for the use of $RANDOM below to
36
# seed the generator. However, shells that don't provide it won't be
37
# a problem since gawk will then seed the generator with the time of day,
38
# as srand() will be called without an argument.
40
# large NSAMPLES and NRUNS will bring any correlation out of the noise better
41
NSAMPLES=1024; MAX_ALLOWED_SIGMA=5; NRUNS=50;
45
nsamples=('$NSAMPLES');
46
max_allowed_sigma=('$MAX_ALLOWED_SIGMA');
48
for(tau=0;tau<nsamples/2;tau++) corr[tau]=0;
50
for(run=0;run<nruns;run++) {
53
# Fill an array with a sequence of samples that are a
54
# function of pairs of rand() values.
56
for(i=0;i<nsamples;i++) {
57
samp[i]=((rand()-0.5)*(rand()-0.5))^2;
61
# Subtract off the mean of the sequence:
64
for(i=0;i<nsamples;i++) samp[i]=samp[i]-mean;
66
# Calculate an autocorrelation function on the sequence.
67
# Because the values of rand() should be independent, there
68
# should be no peaks in the autocorrelation.
70
for(tau=0;tau<nsamples/2;tau++) {
72
for(i=0;i<nsamples/2;i++) sum=sum+samp[i]*samp[i+tau];
73
corr[tau]=corr[tau]+sum;
77
# Normalize the autocorrelation to the tau=0 value.
80
for(tau=0;tau<nsamples/2;tau++) corr[tau]=corr[tau]/max_corr;
82
# OPTIONALLY Print out the autocorrelation values:
84
# for(tau=0;tau<nsamples/2;tau++) print tau, corr[tau] > "pairpower_corr.data";
86
# Calculate the sigma for the non-zero tau values:
90
for(tau=1;tau<nsamples/2;tau++) power_sum=power_sum+(corr[tau])^2;
92
sigma=sqrt(power_sum/(nsamples/2-1));
94
# See if any of the correlations exceed a reasonable number of sigma:
97
for(tau=1;tau<nsamples/2;tau++) {
98
if ( abs(corr[tau])/sigma > max_allowed_sigma ) {
99
print "Tau=", tau ", Autocorr=", corr[tau]/sigma, "sigma";
110
function abs(abs_input) { return(sqrt(abs_input^2)) ; }