source code

/**
* Robert Sedgewick Algorithms in C++<br>
* Chapter 35 Random Numbers<br>
* Problem6<br>
* 加法的合同法を使って1000未満の整数を1000個発生し、それらがランダムかどうかを判定する検定法を設計し適用せよ。
*/
public class Chapter35RandomNumbersProblem6 {
public static void main(String[] args) {
Chapter35RandomNumbersProblem6 o = new Chapter35RandomNumbersProblem6();
for (int i = 0; i < 100; i++) {
o.run();
}
}

private void run() {
AdditiveCongruentialMethod rng = new AdditiveCongruentialMethod(Math.abs((int) System.nanoTime()));
int N = 1000;
int r = 1000;
int[] a = new int[N];
int test = 0;
int valid = 0;
for (; test < 1000; test++) {
for (int i = 0; i < a.length; i++) {
a[i] = rng.nextInt(r);
}
if (simpletest(a, r)) {
valid++;
}
}
System.out.format("%10.5f\n", (double) valid / test);
}

private boolean simpletest(int[] a, int r) {
return Math.abs(chisquare(a, r) - r) <= 2 * Math.sqrt(r);
}

private double chisquare(int[] a, int r) {
int N = a.length;
int[] count = new int[r];
for (int i = 0; i < N; i++) {
count[a[i]]++;
}
double t = 0;
for (int i = 0; i < count.length; i++) {
t += count[i] * count[i];
}
return t / ((double) N / r) - N;
}

private class AdditiveCongruentialMethod {
private int[] a = new int[55];
private int j = 0;
private int m = 100000000;
private int m2 = 10000;
private int b = 31415821;

public AdditiveCongruentialMethod(int seed) {
seed = Math.abs(seed) % m;
a[0] = seed;
if (seed < 0 || seed >= m) {
debug("seed", seed);
throw new AssertionError();
}
for (int i = 1; i < a.length; i++) {
a[i] = (mult(a[i - 1], b) + 1) % m;
if (a[i] < 0 || a[i] >= m) {
debug("a[i]", a[i], "a[i-1]", a[i - 1], "b", b);
throw new AssertionError();
}
}
}

public int nextInt(int r) {
j = (j + 1) % 55;
a[j] = (a[(j + 23) % 55] + a[(j + 54) % 55]) % m;
int res = ((a[j] / m2) * r) / m2;
if (res < 0 || res >= r) {
throw new AssertionError();
}
return res;
}

private int mult(int a, int b) {
int a1 = a / m2;
int a2 = a % m2;
int b1 = b / m2;
int b2 = b % m2;
if (a * b != (a1 * m2 + a2) * (b1 * m2 + b2)) {
throw new AssertionError();
}
if ((a1 * m2 + a2) * (b1 * m2 + b2) != (a1 * b1) * m + (a2 * b1 + a1 * b2) * m2 + (a2 * b2)) {
throw new AssertionError();
}
int res = (((a2 * b1 + a1 * b2) % m2) * m2 + (a2 * b2)) % m;
if (res < 0 || res >= m) {
debug("a", a, "b", b, "res", res);
throw new AssertionError();
}
return res;
}
}

private static final boolean DEBUG = true;

private static final void debug(Object... o) {
if (DEBUG)
System.err.println(Arrays.deepToString(o));
}

private static final void message(Object... o) {
System.err.println(Arrays.deepToString(o));
}

}