1 /* 2 * Copyright 2013 Google Inc. 3 * 4 * Use of this source code is governed by a BSD-style license that can be 5 * found in the LICENSE file. 6 */ 7 8 #include "SkRandom.h" 9 #include "SkTSort.h" 10 #include "Test.h" 11 12 static bool anderson_darling_test(double p[32]) { 13 // Min and max Anderson-Darling values allowable for k=32 14 const double kADMin32 = 0.202; // p-value of ~0.1 15 const double kADMax32 = 3.89; // p-value of ~0.99 16 17 // sort p values 18 SkTQSort<double>(p, p + 31); 19 20 // and compute Anderson-Darling statistic to ensure these are uniform 21 double s = 0.0; 22 for(int k = 0; k < 32; k++) { 23 double v = p[k]*(1.0 - p[31-k]); 24 if (v < 1.0e-30) { 25 v = 1.0e-30; 26 } 27 s += (2.0*(k+1)-1.0)*log(v); 28 } 29 double a2 = -32.0 - 0.03125*s; 30 31 return (kADMin32 < a2 && a2 < kADMax32); 32 } 33 34 static bool chi_square_test(int bins[256], int e) { 35 // Min and max chisquare values allowable 36 const double kChiSqMin256 = 206.3179; // probability of chance = 0.99 with k=256 37 const double kChiSqMax256 = 311.5603; // probability of chance = 0.01 with k=256 38 39 // compute chi-square 40 double chi2 = 0.0; 41 for (int j = 0; j < 256; ++j) { 42 double delta = bins[j] - e; 43 chi2 += delta*delta/e; 44 } 45 46 return (kChiSqMin256 < chi2 && chi2 < kChiSqMax256); 47 } 48 49 // Approximation to the normal distribution CDF 50 // From Waissi and Rossin, 1996 51 static double normal_cdf(double z) { 52 double t = ((-0.0004406*z*z* + 0.0418198)*z*z + 0.9)*z; 53 t *= -1.77245385091; // -sqrt(PI) 54 double p = 1.0/(1.0 + exp(t)); 55 56 return p; 57 } 58 59 static void test_random_byte(skiatest::Reporter* reporter, int shift) { 60 int bins[256]; 61 memset(bins, 0, sizeof(int)*256); 62 63 SkRandom rand; 64 for (int i = 0; i < 256*10000; ++i) { 65 bins[(rand.nextU() >> shift) & 0xff]++; 66 } 67 68 REPORTER_ASSERT(reporter, chi_square_test(bins, 10000)); 69 } 70 71 static void test_random_float(skiatest::Reporter* reporter) { 72 int bins[256]; 73 memset(bins, 0, sizeof(int)*256); 74 75 SkRandom rand; 76 for (int i = 0; i < 256*10000; ++i) { 77 float f = rand.nextF(); 78 REPORTER_ASSERT(reporter, 0.0f <= f && f < 1.0f); 79 bins[(int)(f*256.f)]++; 80 } 81 REPORTER_ASSERT(reporter, chi_square_test(bins, 10000)); 82 83 double p[32]; 84 for (int j = 0; j < 32; ++j) { 85 float f = rand.nextF(); 86 REPORTER_ASSERT(reporter, 0.0f <= f && f < 1.0f); 87 p[j] = f; 88 } 89 REPORTER_ASSERT(reporter, anderson_darling_test(p)); 90 } 91 92 // This is a test taken from tuftests by Marsaglia and Tsang. The idea here is that 93 // we are using the random bit generated from a single shift position to generate 94 // "strings" of 16 bits in length, shifting the string and adding a new bit with each 95 // iteration. We track the numbers generated. The ones that we don't generate will 96 // have a normal distribution with mean ~24108 and standard deviation ~127. By 97 // creating a z-score (# of deviations from the mean) for one iteration of this step 98 // we can determine its probability. 99 // 100 // The original test used 26 bit strings, but is somewhat slow. This version uses 16 101 // bits which is less rigorous but much faster to generate. 102 static double test_single_gorilla(skiatest::Reporter* reporter, int shift) { 103 const int kWordWidth = 16; 104 const double kMean = 24108.0; 105 const double kStandardDeviation = 127.0; 106 const int kN = (1 << kWordWidth); 107 const int kNumEntries = kN >> 5; // dividing by 32 108 unsigned int entries[kNumEntries]; 109 110 SkRandom rand; 111 memset(entries, 0, sizeof(unsigned int)*kNumEntries); 112 // pre-seed our string value 113 int value = 0; 114 for (int i = 0; i < kWordWidth-1; ++i) { 115 value <<= 1; 116 unsigned int rnd = rand.nextU(); 117 value |= ((rnd >> shift) & 0x1); 118 } 119 120 // now make some strings and track them 121 for (int i = 0; i < kN; ++i) { 122 value = SkLeftShift(value, 1); 123 unsigned int rnd = rand.nextU(); 124 value |= ((rnd >> shift) & 0x1); 125 126 int index = value & (kNumEntries-1); 127 SkASSERT(index < kNumEntries); 128 int entry_shift = (value >> (kWordWidth-5)) & 0x1f; 129 entries[index] |= (0x1 << entry_shift); 130 } 131 132 // count entries 133 int total = 0; 134 for (int i = 0; i < kNumEntries; ++i) { 135 unsigned int entry = entries[i]; 136 while (entry) { 137 total += (entry & 0x1); 138 entry >>= 1; 139 } 140 } 141 142 // convert counts to normal distribution z-score 143 double z = ((kN-total)-kMean)/kStandardDeviation; 144 145 // compute probability from normal distibution CDF 146 double p = normal_cdf(z); 147 148 REPORTER_ASSERT(reporter, 0.01 < p && p < 0.99); 149 return p; 150 } 151 152 static void test_gorilla(skiatest::Reporter* reporter) { 153 154 double p[32]; 155 for (int bit_position = 0; bit_position < 32; ++bit_position) { 156 p[bit_position] = test_single_gorilla(reporter, bit_position); 157 } 158 159 REPORTER_ASSERT(reporter, anderson_darling_test(p)); 160 } 161 162 static void test_range(skiatest::Reporter* reporter) { 163 SkRandom rand; 164 165 // just to make sure we don't crash in this case 166 (void) rand.nextRangeU(0, 0xffffffff); 167 168 // check a case to see if it's uniform 169 int bins[256]; 170 memset(bins, 0, sizeof(int)*256); 171 for (int i = 0; i < 256*10000; ++i) { 172 unsigned int u = rand.nextRangeU(17, 17+255); 173 REPORTER_ASSERT(reporter, 17 <= u && u <= 17+255); 174 bins[u - 17]++; 175 } 176 177 REPORTER_ASSERT(reporter, chi_square_test(bins, 10000)); 178 } 179 180 DEF_TEST(Random, reporter) { 181 // check uniform distributions of each byte in 32-bit word 182 test_random_byte(reporter, 0); 183 test_random_byte(reporter, 8); 184 test_random_byte(reporter, 16); 185 test_random_byte(reporter, 24); 186 187 test_random_float(reporter); 188 189 test_gorilla(reporter); 190 191 test_range(reporter); 192 } 193
