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
12static 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
34static 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
51static 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
59static 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
71static 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.
102static 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 <<= 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
152static 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
162static 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
180DEF_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