Halton Sequences霍尔顿序列生成器，计算Pi

这个程序的原理是这样的。假如有一个边长为1的正方形。以正方形的一个端点为圆心，以1为半径，画一个圆弧，于是在正方形内就有了一个直角扇形。在正方形里随机生成若干的点，则有些点是在扇形内，有些点是在扇形外。正方形的面积是1，扇形的面积是0.25*Pi。设点的数量一共是n，扇形内的点数量是nc，在点足够多足够密集的情况下，会近似有nc/n的比值约等于扇形面积与正方形面积的比值，也就是nc/n = 0.25*Pi/1，即Pi = 4*nc/n。

如何生成随机点？最简单的方式是在[0, 1]的区间内每次生成两个随机小数作为随机点的x和y坐标。可惜这种生成方式效果不够好，随机点之间有间隙过大和重叠的可能，会让计算精度不够高。Halton序列算法生成样本点的效果要好得多，更均匀，计算精度比随机生成的点更高，因此这个例子用Halton序列算法生成点集。关于Halton序列可以参考这里

package com.tiger.hadoop;

public class Halton {


/**
* @param args
*/

private static class HaltonSequence {
static final int[] P = {2, 3};
static final int[] K = {63, 40};


private long index;
private double[] x;
private double[][] q;
private int[][] d;


//构造函数
HaltonSequence(long startindex) {
index = startindex;
x = new double[K.length];
q = new double[K.length][];
d = new int[K.length][];
for(int i = 0; i < K.length; i++) {
q[i] = new double[K[i]];
d[i] = new int[K[i]];
}


for(int i = 0; i < K.length; i++) {
long k = index;
x[i] = 0;
for(int j = 0; j < K[i]; j++) {
q[i][j] = (j == 0? 1.0: q[i][j-1])/P[i];
d[i][j] = (int)(k % P[i]);
k = (k - d[i][j])/P[i];
x[i] += d[i][j] * q[i][j];
}
}
}


//生成点
double[] nextPoint() {
index++;
for(int i = 0; i < K.length; i++) {
for(int j = 0; j < K[i]; j++) {
d[i][j]++;
x[i] += q[i][j];
if (d[i][j] < P[i]) {
break;
}
d[i][j] = 0;
x[i] -= (j == 0? 1.0: q[i][j-1]);
}
}
return x;
}
}
public static void main(String[] args) {
HaltonSequence haltonseq=new HaltonSequence(1000);
int innum=0;
long num=100000000;
for (int i=0;i<num;i++){
double x=haltonseq.nextPoint()[0]-0.5;
double y=haltonseq.nextPoint()[1]-0.5;
if(x*x+y*y<0.25)innum++; 

}
System.err.println((double)4*innum/num);
;



}


}