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