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