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