package umontreal.iro.lecuyer.hups;
/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
import umontreal.iro.lecuyer.util.PrintfFormat;
/**
* This class implements digital nets or digital sequences formed by the
* first n = bk points of the Faure sequence in base b.
* Values of n up to 231 are allowed.
* One has r = k.
* The generator matrices are
*
*
*
* Cj = Pj mod b
*
* for
* j = 0,..., s - 1, where
* P is a k×k upper
* triangular matrix whose entry (l, c) is the number of combinations
* of l objects among c,
* for l <= c and is 0 for l > c.
* The matrix
* C0 is the identity,
* C1 = P,
* and the other
* Cj's can be defined recursively via
*
* Cj = PCj-1mod b.
* Our implementation uses the recursion
*
*
*
* Combination(c, l ) = Combination(c - 1, l ) + Combination(c - 1, l - 1)
*
* to evaluate the binomial coefficients in the matrices
* Cj,
* as suggested by Fox.
* The entries xj, l, c of
* Cj are computed as follows:
*
*
*
*
* xj, c, c
* =
* 1
* for c = 0,..., k - 1,
* xj, 0, c
* =
* jxj, 0, c-1
* for c = 1,..., k - 1,
* xj, l, c
* =
* xj, l-1, c-1 + jxj, l, c-1
* for 2 <= c < l <= k - 1,
* xj, l, c
* =
* 0
* for c > l or l >= k.
*
*
*
*
* For any integer m > 0 and ν >= 0, if we look at the
* vector
* (ui, j, 1,..., ui, j, m) (the first m digits
* of coordinate j of the output) when i goes from
* νbm to
* (ν +1)bm - 1, this vector takes each of its bm
* possible values exactly once.
* In particular, for ν = 0, ui, j visits each value in the
* set
* {0, 1/bm, 2/bm,...,(bm -1)/bm} exactly once, so all
* one-dimensional projections of the point set are identical.
* However, the values are visited in a different order for
* the different values of j (otherwise all coordinates would be identical).
* For j = 0, they are visited in the same
* order as in the van der Corput sequence in base b.
*
*
* An important property of Faure nets is that for any integers m > 0
* and ν >= 0, the point set
*
* {ui for
* i = νbm,...,(ν +1)bm -1}
* is a (0, m, s)-net in base b.
* In particular, for n = bk, the first n points form a
* (0, k, s)-net in base b.
* The Faure nets are also projection-regular and dimension-stationary.
*
*
* To obtain digital nets from the generalized Faure sequence
* , where
* Pj is left-multiplied by some
* invertible matrix
* Aj, it suffices to apply an appropriate
* matrix scramble (e.g., via {@link #leftMatrixScramble leftMatrixScramble}).
* This changes the order in which ui, j visits its different
* values, for each coordinate j, but does not change the set of values
* that are visited. The (0, m, s)-net property stated above remains valid.
*
*/
public class FaureSequence extends DigitalSequence {
// Maximum dimension for the case where b is not specified.
// Can be extended by extending the precomputed array prime[].
private static final int MAXDIM = 500;
// For storing the generator matrices for given dim and numPoints.
private int[][][] v;
/**
* Constructs a digital net in base b,
* with n = bk points and w output digits,
* in dim dimensions.
* The points are the first n points of the Faure sequence.
* The generator matrices
* Cj are r×k.
* Unless, one plans to apply a randomization on more than k digits
* (e.g., a random digital shift for w > k digits, or a linear
* scramble yielding r > k digits), one should
* take w = r = k for better computational efficiency.
* Restrictions: dim <= 500 and
* bk <= 231.
*
* @param b base
*
* @param k there will be b^k points
*
* @param r number of rows in the generator matrices
*
* @param w number of output digits
*
* @param dim dimension of the point set
*
*
*/
public FaureSequence (int b, int k, int r, int w, int dim) {
init (b, k, r, w, dim);
}
private void init (int b, int k, int r, int w, int dim) {
if (dim < 1)
throw new IllegalArgumentException
("Dimension for FaureSequence must be > 1");
if ((double)k * Math.log ((double)b) > (31.0 * Math.log (2.0)))
throw new IllegalArgumentException
("Trying to construct a FaureSequence with too many points");
if (r < k || w < r)
throw new IllegalArgumentException
("One must have k <= r <= w for FaureSequence");
this.b = b;
numCols = k;
numRows = r;
outDigits = w;
this.dim = dim;
int i, j;
numPoints = b;
for (i=1; i(b, k, w, w, dim)
* with base b equal to the smallest prime larger or equal to dim,
* and with at least n points.
*
* The values of k, r, and w are taken as
*
* k = ceil(logbn) and
*
* r = w = max(k, floor(30/log2b)).
*
* @param n minimal number of points
*
* @param dim dimension of the point set
*
*
*/
public FaureSequence (int n, int dim) {
if ((dim < 1) || (dim > MAXDIM))
throw new IllegalArgumentException
("Dimension for Faure net must be > 1 and < " + MAXDIM);
b = getSmallestPrime (dim);
numCols = (int) Math.ceil (Math.log ((double) n)
/ Math.log ((double) b));
outDigits = (int) Math.floor (Math.log ((double)(1 << (MAXBITS - 1)))
/ Math.log ((double)b));
outDigits = Math.max (outDigits, numCols);
numRows = outDigits;
init (b, numCols, numRows, outDigits, dim);
}
public String toString() {
StringBuffer sb = new StringBuffer ("Faure sequence:" +
PrintfFormat.NEWLINE);
sb.append (super.toString());
return sb.toString();
}
public void extendSequence (int k) {
init (b, k, numRows, outDigits, dim);
}
// Fills up the generator matrices in genMat for a Faure sequence.
// See Glasserman (2004), \cite{fGLA04a}, page 301.
private void initGenMat() {
int j, c, l;
// Initialize C_0 to the identity (for first coordinate).
for (c = 0; c < numCols; c++) {
for (l = 0; l < numRows; l++)
genMat[c][l] = 0;
genMat[c][c] = 1;
}
// Compute C_1, ... ,C_{dim-1}.
for (j = 1; j < dim; j++) {
genMat[j*numCols][0] = 1;
for (c = 1; c < numCols; c++) {
genMat[j*numCols+c][c] = 1;
genMat[j*numCols+c][0] = (j * genMat[j*numCols+c-1][0]) % b;
}
for (c = 2; c < numCols; c++) {
for (l = 1; l < c; l++)
genMat[j*numCols+c][l] = (genMat[j*numCols+c-1][l-1]
+ j * genMat[j*numCols+c-1][l]) % b;
}
for (c = 0; c < numCols; c++)
for (l = c+1; l < numRows; l++)
genMat[j*numCols+c][l] = 0;
}
/*
for (j = 0; j < dim; j++) {
for (l = 0; l < numRows; l++) {
for (c = 0; c < numCols; c++)
System.out.print (" " + genMat[j*numCols+c][l]);
System.out.println ("");
}
System.out.println ("");
}
*/
}
/*
// Fills up the generator matrices in genMat for a Faure net.
// See Glasserman (2004), \cite{fGLA04a}, page 301.
protected void initGenMatNet() {
int j, c, l, start;
// Initialize C_0 to the reflected identity (for first coordinate).
for (c = 0; c < numCols; c++) {
for (l = 0; l < numRows; l++)
genMat[c][l] = 0;
genMat[c][numCols-c-1] = 1;
}
// Initialize C_1 to the identity (for second coordinate).
for (c = 0; c < numCols; c++) {
for (l = 0; l < numRows; l++)
genMat[numCols+c][l] = 0;
genMat[numCols+c][c] = 1;
}
// Compute C_2, ... ,C_{dim-1}.
for (j = 2; j < dim; j++) {
start = j * numCols;
genMat[start][0] = 1;
for (c = 1; c < numCols; c++) {
genMat[start+c][c] = 1;
genMat[start+c][0] = ((j-1) * genMat[start+c-1][0]) % b;
}
for (c = 2; c < numCols; c++) {
for (l = 1; l < c; l++)
genMat[start+c][l] = (genMat[start+c-1][l-1]
+ (j-1) * genMat[start+c-1][l]) % b;
}
for (c = 0; c < numCols; c++)
for (l = c+1; l < numRows; l++)
genMat[start+c][l] = 0;
}
}
*/
// Returns the smallest prime larger or equal to d.
private int getSmallestPrime (int d) {
return primes[d-1];
}
// Gives the appropriate prime base for each dimension.
// Perhaps should be internal to getPrime, and non-static, to avoid
// wasting time and memory when this array is not needed ???
static final int primes[] =
{2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,
23,23,23,29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,
41,43,43,47,47,47,47,53,53,53,53,53,53,59,59,59,59,59,59,61,
61,67,67,67,67,67,67,71,71,71,71,73,73,79,79,79,79,79,79,83,
83,83,83,89,89,89,89,89,89,97,97,97,97,97,97,97,97,101,101,101,
101,103,103,107,107,107,107,109,109,113,113,113,113,127,127,127,127,127,127,127,
127,127,127,127,127,127,127,131,131,131,131,137,137,137,137,137,137,139,139,149,
149,149,149,149,149,149,149,149,149,151,151,157,157,157,157,157,157,163,163,163,
163,163,163,167,167,167,167,173,173,173,173,173,173,179,179,179,179,179,179,181,
181,191,191,191,191,191,191,191,191,191,191,193,193,197,197,197,197,199,199,211,
211,211,211,211,211,211,211,211,211,211,211,223,223,223,223,223,223,223,223,223,
223,223,223,227,227,227,227,229,229,233,233,233,233,239,239,239,239,239,239,241,
241,251,251,251,251,251,251,251,251,251,251,257,257,257,257,257,257,263,263,263,
263,263,263,269,269,269,269,269,269,271,271,277,277,277,277,277,277,281,281,281,
281,283,283,293,293,293,293,293,293,293,293,293,293,307,307,307,307,307,307,307,
307,307,307,307,307,307,307,311,311,311,311,313,313,317,317,317,317,331,331,331,
331,331,331,331,331,331,331,331,331,331,331,337,337,337,337,337,337,347,347,347,
347,347,347,347,347,347,347,349,349,353,353,353,353,359,359,359,359,359,359,367,
367,367,367,367,367,367,367,373,373,373,373,373,373,379,379,379,379,379,379,383,
383,383,383,389,389,389,389,389,389,397,397,397,397,397,397,397,397,401,401,401,
401,409,409,409,409,409,409,409,409,419,419,419,419,419,419,419,419,419,419,421,
421,431,431,431,431,431,431,431,431,431,431,433,433,439,439,439,439,439,439,443,
443,443,443,449,449,449,449,449,449,457,457,457,457,457,457,457,457,461,461,461,
461,463,463,467,467,467,467,479,479,479,479,479,479,479,479,479,479,479,479,487,
487,487,487,487,487,487,487,491,491,491,491,499,499,499,499,499,499,499,499,503};
}
Ads help maintain this website.