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 "include/utils/SkRandom.h"
9 #include "src/core/SkTSort.h"
10 #include "tests/Test.h"
11 
anderson_darling_test(double p[32])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 + 32);
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 
chi_square_test(int bins[256],int e)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
normal_cdf(double z)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 
test_random_byte(skiatest::Reporter * reporter,int shift)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 
test_random_float(skiatest::Reporter * reporter)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.
test_single_gorilla(skiatest::Reporter * reporter,int shift)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 
test_gorilla(skiatest::Reporter * reporter)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 
test_range(skiatest::Reporter * reporter)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 
DEF_TEST(Random,reporter)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