aboutsummaryrefslogtreecommitdiffhomepage
path: root/tests/RandomTest.cpp
blob: 51408e960536b6ee608a2ab481960d305a4aee01 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

/*
 * Copyright 2013 Google Inc.
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */

#include "Test.h"
#include "SkRandom.h"
#include "SkTSort.h"

static bool anderson_darling_test(double p[32]) {
    // Min and max Anderson-Darling values allowable for k=32
    const double kADMin32 = 0.202;        // p-value of ~0.1
    const double kADMax32 = 3.89;         // p-value of ~0.99

    // sort p values
    SkTQSort<double>(p, p + 31);

    // and compute Anderson-Darling statistic to ensure these are uniform
    double s = 0.0;
    for(int k = 0; k < 32; k++) {
        double v = p[k]*(1.0 - p[31-k]);
        if (v < 1.0e-30) {
           v = 1.0e-30;
        }
        s += (2.0*(k+1)-1.0)*log(v);
    }
    double a2 = -32.0 - 0.03125*s;

    return (kADMin32 < a2 && a2 < kADMax32);
}

static bool chi_square_test(int bins[256], int e) {
    // Min and max chisquare values allowable
    const double kChiSqMin256 = 206.3179;        // probability of chance = 0.99 with k=256
    const double kChiSqMax256 = 311.5603;        // probability of chance = 0.01 with k=256

    // compute chi-square
    double chi2 = 0.0;
    for (int j = 0; j < 256; ++j) {
        double delta = bins[j] - e;
        chi2 += delta*delta/e;
    }

    return (kChiSqMin256 < chi2 && chi2 < kChiSqMax256);
}

// Approximation to the normal distribution CDF
// From Waissi and Rossin, 1996
static double normal_cdf(double z) {
    double t = ((-0.0004406*z*z* + 0.0418198)*z*z + 0.9)*z;
    t *= -1.77245385091;  // -sqrt(PI)
    double p = 1.0/(1.0 + exp(t));

    return p;
}

static void test_random_byte(skiatest::Reporter* reporter, int shift) {
    int bins[256];
    memset(bins, 0, sizeof(int)*256);

    SkMWCRandom rand;
    for (int i = 0; i < 256*10000; ++i) {
        bins[(rand.nextU() >> shift) & 0xff]++;
    }

    REPORTER_ASSERT(reporter, chi_square_test(bins, 10000));
}

static void test_random_float(skiatest::Reporter* reporter) {
    int bins[256];
    memset(bins, 0, sizeof(int)*256);

    SkMWCRandom rand;
    for (int i = 0; i < 256*10000; ++i) {
        float f = rand.nextF();
        REPORTER_ASSERT(reporter, 0.0f <= f && f < 1.0f);
        bins[(int)(f*256.f)]++;
    }
    REPORTER_ASSERT(reporter, chi_square_test(bins, 10000));

    double p[32];
    for (int j = 0; j < 32; ++j) {
        float f = rand.nextF();
        REPORTER_ASSERT(reporter, 0.0f <= f && f < 1.0f);
        p[j] = f;
    }
    REPORTER_ASSERT(reporter, anderson_darling_test(p));
}

// This is a test taken from tuftests by Marsaglia and Tsang. The idea here is that
// we are using the random bit generated from a single shift position to generate
// "strings" of 16 bits in length, shifting the string and adding a new bit with each
// iteration. We track the numbers generated. The ones that we don't generate will
// have a normal distribution with mean ~24108 and standard deviation ~127. By
// creating a z-score (# of deviations from the mean) for one iteration of this step
// we can determine its probability.
//
// The original test used 26 bit strings, but is somewhat slow. This version uses 16
// bits which is less rigorous but much faster to generate.
static double test_single_gorilla(skiatest::Reporter* reporter, int shift) {
    const int kWordWidth = 16;
    const double kMean = 24108.0;
    const double kStandardDeviation = 127.0;
    const int kN = (1 << kWordWidth);
    const int kNumEntries = kN >> 5;  // dividing by 32
    unsigned int entries[kNumEntries];

    SkMWCRandom rand;
    memset(entries, 0, sizeof(unsigned int)*kNumEntries);
    // pre-seed our string value
    int value = 0;
    for (int i = 0; i < kWordWidth-1; ++i) {
        value <<= 1;
        unsigned int rnd = rand.nextU();
        value |= ((rnd >> shift) & 0x1);
    }

    // now make some strings and track them
    for (int i = 0; i < kN; ++i) {
        value <<= 1;
        unsigned int rnd = rand.nextU();
        value |= ((rnd >> shift) & 0x1);

        int index = value & (kNumEntries-1);
        SkASSERT(index < kNumEntries);
        int entry_shift = (value >> (kWordWidth-5)) & 0x1f;
        entries[index] |= (0x1 << entry_shift);
    }

    // count entries
    int total = 0;
    for (int i = 0; i < kNumEntries; ++i) {
        unsigned int entry = entries[i];
        while (entry) {
            total += (entry & 0x1);
            entry >>= 1;
        }
    }

    // convert counts to normal distribution z-score
    double z = ((kN-total)-kMean)/kStandardDeviation;

    // compute probability from normal distibution CDF
    double p = normal_cdf(z);

    REPORTER_ASSERT(reporter, 0.01 < p && p < 0.99);
    return p;
}

static void test_gorilla(skiatest::Reporter* reporter) {

    double p[32];
    for (int bit_position = 0; bit_position < 32; ++bit_position) {
        p[bit_position] = test_single_gorilla(reporter, bit_position);
    }

    REPORTER_ASSERT(reporter, anderson_darling_test(p));
}

static void test_range(skiatest::Reporter* reporter) {
    SkMWCRandom rand;

    // just to make sure we don't crash in this case
    (void) rand.nextRangeU(0, 0xffffffff);

    // check a case to see if it's uniform
    int bins[256];
    memset(bins, 0, sizeof(int)*256);
    for (int i = 0; i < 256*10000; ++i) {
        unsigned int u = rand.nextRangeU(17, 17+255);
        REPORTER_ASSERT(reporter, 17 <= u && u <= 17+255);
        bins[u - 17]++;
    }

    REPORTER_ASSERT(reporter, chi_square_test(bins, 10000));
}

static void TestRandom(skiatest::Reporter* reporter) {
    // check uniform distributions of each byte in 32-bit word
    test_random_byte(reporter, 0);
    test_random_byte(reporter, 8);
    test_random_byte(reporter, 16);
    test_random_byte(reporter, 24);

    test_random_float(reporter);

    test_gorilla(reporter);

    test_range(reporter);
}

#include "TestClassDef.h"
DEFINE_TESTCLASS("Random", RandomTestClass, TestRandom)