RandomTest.cpp 5.7 KB

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