List of usage examples for org.apache.commons.math.random MersenneTwister nextDouble
public double nextDouble()
From source file:cz.paulrz.montecarlo.random.Sobol.java
/** * @param dim dimension of the Sobol sequence. *//* ww w . ja va2 s. com*/ public Sobol(int dim) throws Exception { this.dim = dim; if (dim > 300) { throw new Exception("Sobol sequence only implemented for dimension at most 300." + "\nExiting."); } // array of polynomial degrees and coefficient arrays g = new int[dim]; p = new int[dim][]; // degree zero g[0] = 0; p[0] = new int[1]; p[0][0] = 1; // positive degree polynomials read from pp int k = 1; for (int i = 0; i < pp.length; i++) for (int j = 0; j < pp[i].length; j++) if (k < dim) { g[k] = i + 1; read_prim_pol(g[k], pp[i][j], k); k++; } else break; /* DEBUG: print the primitive polynomials: System.out.println("\nPrimitive polynomials:\n"); for(k=0;k<50;k++) { for(int j=0;j<=g[k];j++)System.out.print(p[k][j]); System.out.println(""); } */ // initialize the array of direction integers v = new long[dim][bits]; for (int j = 0; j < bits; j++) v[0][j] = (1L << (bits - j - 1)); if (dim > 1) v[1][0] = (1L << bits - 1); if (dim > 2) { v[2][0] = (1L << bits - 1); v[2][1] = (1L << bits - 2); } if (dim > 3) { v[3][0] = (1L << bits - 1); v[3][1] = (3L << bits - 2); v[3][2] = (7L << bits - 3); } if (dim > 4) { v[4][0] = (1L << bits - 1); v[4][1] = (1L << bits - 2); v[4][2] = (5L << bits - 3); } if (dim > 5) { v[5][0] = (1L << bits - 1); v[5][1] = (3L << bits - 2); v[5][2] = (1L << bits - 3); v[5][3] = (1L << bits - 4); } if (dim > 6) { v[6][0] = (1L << bits - 1); v[6][1] = (1L << bits - 2); v[6][2] = (3L << bits - 3); v[6][3] = (7L << bits - 4); } if (dim > 7) { v[7][0] = (1L << bits - 1); v[7][1] = (3L << bits - 2); v[7][2] = (3L << bits - 3); v[7][3] = (9L << bits - 4); v[7][4] = (9L << bits - 5); } if (dim > 8) { v[8][0] = (1L << bits - 1); v[8][1] = (3L << bits - 2); v[8][2] = (7L << bits - 3); v[8][3] = (7L << bits - 4); v[8][4] = (21L << bits - 5); } if (dim > 9) { v[9][0] = (1L << bits - 1); v[9][1] = (1L << bits - 2); v[9][2] = (5L << bits - 3); v[9][3] = (11L << bits - 4); v[9][4] = (27L << bits - 5); } if (dim > 10) { v[10][0] = (1L << bits - 1); v[10][1] = (1L << bits - 2); v[10][2] = (7L << bits - 3); v[10][3] = (3L << bits - 4); v[10][4] = (29L << bits - 5); } if (dim > 11) { v[11][0] = (1L << bits - 1); v[11][1] = (3L << bits - 2); v[11][2] = (7L << bits - 3); v[11][3] = (13L << bits - 4); v[11][4] = (3L << bits - 5); } if (dim > 12) { v[12][0] = (1L << bits - 1); v[12][1] = (3L << bits - 2); v[12][2] = (5L << bits - 3); v[12][3] = (1L << bits - 4); v[12][4] = (15L << bits - 5); } if (dim > 13) { v[13][0] = (1L << bits - 1); v[13][1] = (1L << bits - 2); v[13][2] = (1L << bits - 3); v[13][3] = (9L << bits - 4); v[13][4] = (23L << bits - 5); v[13][5] = (37L << bits - 6); } if (dim > 14) { v[14][0] = (1L << bits - 1); v[14][1] = (1L << bits - 2); v[14][2] = (3L << bits - 3); v[14][3] = (13L << bits - 4); v[14][4] = (11L << bits - 5); v[14][5] = (7L << bits - 6); } if (dim > 15) { v[15][0] = (1L << bits - 1); v[15][1] = (3L << bits - 2); v[15][2] = (3L << bits - 3); v[15][3] = (5L << bits - 4); v[15][4] = (19L << bits - 5); v[15][5] = (33L << bits - 6); } if (dim > 16) { v[16][0] = (1L << bits - 1); v[16][1] = (1L << bits - 2); v[16][2] = (7L << bits - 3); v[16][3] = (13L << bits - 4); v[16][4] = (25L << bits - 5); v[16][5] = (5L << bits - 6); } if (dim > 17) { v[17][0] = (1L << bits - 1); v[17][1] = (1L << bits - 2); v[17][2] = (1L << bits - 3); v[17][3] = (13L << bits - 4); v[17][4] = (15L << bits - 5); v[17][5] = (39L << bits - 6); } if (dim > 18) { v[18][0] = (1L << bits - 1); v[18][1] = (3L << bits - 2); v[18][2] = (5L << bits - 3); v[18][3] = (11L << bits - 4); v[18][4] = (7L << bits - 5); v[18][5] = (11L << bits - 6); } if (dim > 19) { v[19][0] = (1L << bits - 1); v[19][1] = (3L << bits - 2); v[19][2] = (1L << bits - 3); v[19][3] = (7L << bits - 4); v[19][4] = (3L << bits - 5); v[19][5] = (23L << bits - 6); v[19][6] = (79L << bits - 7); } if (dim > 20) { v[20][0] = (1L << bits - 1); v[20][1] = (3L << bits - 2); v[20][2] = (1L << bits - 3); v[20][3] = (15L << bits - 4); v[20][4] = (17L << bits - 5); v[20][5] = (63L << bits - 6); v[20][6] = (13L << bits - 7); } if (dim > 21) { v[21][0] = (1L << bits - 1); v[21][1] = (3L << bits - 2); v[21][2] = (3L << bits - 3); v[21][3] = (3L << bits - 4); v[21][4] = (25L << bits - 5); v[21][5] = (17L << bits - 6); v[21][6] = (115L << bits - 7); } if (dim > 22) { v[22][0] = (1L << bits - 1); v[22][1] = (3L << bits - 2); v[22][2] = (7L << bits - 3); v[22][3] = (9L << bits - 4); v[22][4] = (31L << bits - 5); v[22][5] = (29L << bits - 6); v[22][6] = (17L << bits - 7); } if (dim > 23) { v[23][0] = (1L << bits - 1); v[23][1] = (1L << bits - 2); v[23][2] = (3L << bits - 3); v[23][3] = (15L << bits - 4); v[23][4] = (29L << bits - 5); v[23][5] = (15L << bits - 6); v[23][6] = (41L << bits - 7); } if (dim > 24) { v[24][0] = (1L << bits - 1); v[24][1] = (3L << bits - 2); v[24][2] = (1L << bits - 3); v[24][3] = (9L << bits - 4); v[24][4] = (5L << bits - 5); v[24][5] = (21L << bits - 6); v[24][6] = (119L << bits - 7); } if (dim > 25) { v[25][0] = (1L << bits - 1); v[25][1] = (1L << bits - 2); v[25][2] = (5L << bits - 3); v[25][3] = (5L << bits - 4); v[25][4] = (1L << bits - 5); v[25][5] = (27L << bits - 6); v[25][6] = (33L << bits - 7); } if (dim > 26) { v[26][0] = (1L << bits - 1); v[26][1] = (1L << bits - 2); v[26][2] = (3L << bits - 3); v[26][3] = (1L << bits - 4); v[26][4] = (23L << bits - 5); v[26][5] = (13L << bits - 6); v[26][6] = (75L << bits - 7); } if (dim > 27) { v[27][0] = (1L << bits - 1); v[27][1] = (1L << bits - 2); v[27][2] = (7L << bits - 3); v[27][3] = (7L << bits - 4); v[27][4] = (19L << bits - 5); v[27][5] = (25L << bits - 6); v[27][6] = (105L << bits - 7); } if (dim > 28) { v[28][0] = (1L << bits - 1); v[28][1] = (3L << bits - 2); v[28][2] = (5L << bits - 3); v[28][3] = (5L << bits - 4); v[28][4] = (21L << bits - 5); v[28][5] = (9L << bits - 6); v[28][6] = (7L << bits - 7); } if (dim > 29) { v[29][0] = (1L << bits - 1); v[29][1] = (1L << bits - 2); v[29][2] = (1L << bits - 3); v[29][3] = (15L << bits - 4); v[29][4] = (5L << bits - 5); v[29][5] = (49L << bits - 6); v[29][6] = (59L << bits - 7); } if (dim > 30) { v[30][0] = (1L << bits - 1); v[30][1] = (3L << bits - 2); v[30][2] = (5L << bits - 3); v[30][3] = (15L << bits - 4); v[30][4] = (17L << bits - 5); v[30][5] = (19L << bits - 6); v[30][6] = (21L << bits - 7); } if (dim > 31) { v[31][0] = (1L << bits - 1); v[31][1] = (1L << bits - 2); v[31][2] = (7L << bits - 3); v[31][3] = (11L << bits - 4); v[31][4] = (13L << bits - 5); v[31][5] = (29L << bits - 6); v[31][6] = (3L << bits - 7); } // random initialization in dimension bigger than 32 MersenneTwister mt = new MersenneTwister(); for (k = 32; k < dim; k++) { for (int l = 0; l < g[k]; l++) { double u = mt.nextDouble(); long f = (1L << l + 1), n = (int) (f * u); while (n % 2 == 0) { u = mt.nextDouble(); n = (int) (f * u); } v[k][l] = (n << (bits - l - 1)); } } // end direction integer initialization // computation of direction integer v_kl for k>=degree[k] for (k = 1; k < dim; k++) for (int l = g[k]; l < bits; l++) { long n = (v[k][l - g[k]] >> g[k]); for (int j = 1; j <= g[k]; j++) if (p[k][j] != 0) n = n ^ v[k][l - j]; v[k][l] = n; } // initialize the vector of Sobol integers and Sobol points index = 1; x_int = new long[dim]; for (k = 0; k < dim; k++) x_int[k] = v[k][0]; x = new double[dim]; for (k = 0; k < dim; k++) x[k] = ((double) x_int[k]) / N; }