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