package umontreal.iro.lecuyer.hups;
/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
/**
* This class implements basic methods for working with radical
* inverses of integers in an arbitrary basis b.
* These methods are used in classes that implement
* point sets and sequences based on the van der Corput sequence
* (the Hammersley nets and the Halton sequence, for example).
*
*
* We recall that for a k-digit integer i whose digital
* b-ary expansion is
*
*
*
* i = a0 + a1b + ... + ak-1bk-1,
*
* the radical inverse in base b is
*
*
*
* ψb(i) = a0b-1 + a1b-2 + ... + ak-1b-k.
*
* The van der Corput sequence in base b is the sequence
*
* ψb(0), ψb(1), ψb(2),...
*
*
* Note that ψb(i) cannot always be represented exactly
* as a floating-point number on the computer (e.g., if b is not
* a power of two). For an exact representation, one can use the integer
*
*
*
* bkψb(i) = ak-1 + ... + a1bk-2 + a0bk-1,
*
* which we called the integer radical inverse representation.
* This representation is simply a mirror image of the digits of the
* usual b-ary representation of i.
*
*
* It is common practice to permute locally the values of the
* van der Corput sequence. One way of doing this is to apply a
* permutation to the digits of i before computing ψb(i).
* That is, for a permutation π of the digits
* {0,..., b - 1},
*
*
*
* ψb(i) = ∑r=0k-1arb-r-1
*
* is replaced by
*
*
*
* ∑r=0k-1π(ar)b-r-1.
*
* Applying such a permutation only changes the order in which the
* values of ψb(i) are enumerated. For every integer k, the first
* bk values that are enumerated remain the same (they are the values
* of ψb(i) for
* i = 0,..., bk - 1), but they are enumerated in a
* different order. Often, different permutations π will be applied
* for different coordinates of a point set.
*
*
* The permutation π can be deterministic or random. One (deterministic)
* possibility implemented here is the Faure permutation σb of
*
* {0,..., b - 1} defined as follows.
*
* For b = 2, take
* σ = I, the identical permutation. For even
* b = 2c > 2, take
*
*
*
*
* σ[i]
* =
* 2τ[i] i = 0, 1,…, c - 1
*
* σ[i + c]
* =
* 2τ[i] + 1 i = 0, 1,…, c - 1
*
*
*
*
* where τ[i] is the Faure
* permutation for base c. For odd b = 2c + 1,
* take
*
*
*
*
* σ[c]
* =
* c
*
* σ[i]
* =
* τ[i], if 0 <= τ[i] < c
*
* σ[i]
* =
* τ[i] + 1, if c <= τ[i] < 2c
*
*
*
*
* for
* 0 <= i < c, and take
*
*
*
*
* σ[i]
* =
* τ[i - 1], if 0 <= τ[i - 1] < c
*
* σ[i]
* =
* τ[i - 1] + 1, if c <= τ[i - 1] < 2c
*
*
*
*
* for
* c < i <= 2c, and where τ[i] is the Faure
* permutation for base c. The Faure permutations give very small
* discrepancies (amongst the best known today) for small bases.
*
*/
public class RadicalInverse {
private static final int NP = 168; // First NP primes in table below.
private static final int PLIM = 1000; // NP primes < PLIM
// The first NP prime numbers
private static final int[] PRIMES = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383,
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647,
653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743,
751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997 };
private static final int NRILIM = 1000; // For nextRadicalInverse
private int b; // Base
private double invb; // 1/b
private double logb; // natural log(b)
private int JMAX; // b^JMAX = 2^32
private int co; // Counter for nextRadicalInverse
private double xx; // Current value of x
private long ix; // xx = RadicalInverse (ix)
/*
// For Struckmeier's algorithm
private static final int PARTITION_MAX = 54; // Size of partitionL, bi
private int partM; // Effective size of partitionL, bi
private double[] bi; // bi[i] = (b + 1)/b^i - 1
private double[] partitionL; // L[i] = 1 - 1/b^i
// Boundaries of Struckmeier partitions Lkp of [0, 1]
*/
/**
* Initializes the base of this object to b
* and its first value of x to x0.
*
* @param b Base
*
* @param x0 Initial value of x
*
*/
public RadicalInverse (int b, double x0) {
co = 0;
this.b = b;
invb = 1.0 / b;
logb = Math.log (b);
JMAX = (int) (32.0 * 0.69314718055994530941 / logb);
xx = x0;
ix = computeI (x0);
// initStruckmeier (b);
}
private long computeI (double x) {
// Compute i such that x = RadicalInverse (i).
int[] digits = new int[JMAX]; // Digits of x
int j;
for (j = 0; (j < JMAX) && (x > 0.5e-15); j++) {
x *= b;
digits[j] = (int) x;
x -= digits[j];
}
long i = 0;
for (j = JMAX - 1; j >= 0; j--) {
i = i * b + digits[j];
}
return i;
}
/**
* Provides an elementary method for obtaining the first n prime
* numbers larger than 1.
* Creates and returns an array that contains
* these numbers. This is useful for determining the prime bases for
* the different coordinates of the Halton sequence and Hammersley nets.
*
* @param n number of prime numbers to return
*
* @return an array with the first n prime numbers
*
*/
public static int[] getPrimes (int n) {
// Allocates an array of size n filled with the first n prime
// numbers. n must be positive (n > 0). Routine may fail if not enough
// memory for the array is available. The first prime number is 2.
int i;
boolean moreTests;
int[] prime = new int[n];
int n1 = Math.min (NP, n);
for (i = 0; i < n1; i++)
prime[i] = PRIMES[i];
if (NP < n) {
i = NP;
for (int candidate = PLIM + 1; i < n; candidate += 2) {
prime[i] = candidate;
for (int j = 1; (moreTests = prime[j] <= candidate / prime[j])
&& ((candidate % prime[j]) > 0); j++);
if (! moreTests)
i++;
}
}
return prime;
}
/**
* Computes the radical inverse of i in base b.
* If
* i = ∑r=0k-1arbr, the method computes and returns
*
*
*
* x = ∑r=0k-1arb-r-1.
*
*
* @param b base used for the operation
*
* @param i the value for which the radical inverse will be computed
*
* @return the radical inverse of i in base b
*
*/
public static double radicalInverse (int b, long i) {
double digit, radical, inverse;
digit = radical = 1.0 / (double) b;
for (inverse = 0.0; i > 0; i /= b) {
inverse += digit * (double) (i % b);
digit *= radical;
}
return inverse;
}
/**
* A fast method that incrementally computes the radical inverse xi+1
* in base b from xi = x = ψb(i),
*
* using addition with rigthward carry.
* The parameter invb is equal to 1/b.
* Using long incremental streams (i.e., calling this method several times
* in a row) cause increasing inaccuracy in x. Thus the user should
* recompute the radical inverse directly by calling
* {@link #radicalInverse radicalInverse} every once in a while (i.e. in every few
* thousand calls).
*
* @param invb 1/b where b is the base
*
* @param x the inverse xi
*
* @return the radical inverse xi+1
*
*/
public static double nextRadicalInverse (double invb, double x) {
// Calculates the next radical inverse from x in base b.
// Repeated application causes a loss of accuracy.
// Note that x can be any number from [0,1).
final double ALMOST_ONE = 1.0 - 1e-10;
double nextInverse = x + invb;
if (nextInverse < ALMOST_ONE)
return nextInverse;
else {
double digit1 = invb;
double digit2 = invb * invb;
while (x + digit2 >= ALMOST_ONE) {
digit1 = digit2;
digit2 *= invb;
}
return x + (digit1 - 1.0) + digit2;
}
}
/**
* A fast method that incrementally computes the radical inverse xi+1
* in base b from xi = ψb(i),
* using addition with rigthward carry as described in
* Wang and Hickernell.
* Since using long incremental streams (i.e., calling this method several
* times in a row) cause increasing inaccuracy in x, the method
* recomputes the radical inverse directly from i by calling
* {@link #radicalInverse radicalInverse} once in every 1000 calls.
*
* @return the radical inverse xi+1
*
*
*/
public double nextRadicalInverse () {
// Calculates the next radical inverse from xx in base b.
// Repeated application causes a loss of accuracy.
// For each NRILIM calls, a direct calculation via radicalInverse
// is inserted.
co++;
if (co >= NRILIM) {
co = 0;
ix += NRILIM;
xx = radicalInverse (b, ix);
return xx;
}
final double ALMOST_ONE = 1.0 - 1e-10;
double nextInverse = xx + invb;
if (nextInverse < ALMOST_ONE) {
xx = nextInverse;
return xx;
} else {
double digit1 = invb;
double digit2 = invb * invb;
while (xx + digit2 >= ALMOST_ONE) {
digit1 = digit2;
digit2 *= invb;
}
xx += (digit1 - 1.0) + digit2;
return xx;
}
}
/*
private void initStruckmeier (int b) {
bi = new double[1 + PARTITION_MAX];
partitionL = new double[1 + PARTITION_MAX];
logb = Math.log (b);
partitionL[0] = 0.0;
bi[0] = 1.0;
int i = 0;
while ((i < PARTITION_MAX) && (partitionL[i] < 1.0)) {
++i;
bi[i] = bi[i - 1] / b;
partitionL[i] = 1.0 - bi[i];
}
partM = i - 1;
for (i = 0; i <= partM + 1; ++i)
bi[i] = (b + 1) * bi[i] - 1.0;
}
public double nextRadicalInverse (double x) {
int k;
if (x < partitionL[partM]) {
k = 1;
// Find k such: partitionL[k-1] <= x < partitionL[k]
while (x >= partitionL[k])
++k;
} else { // x >= partitionL [partM]
k = 1 + (int)(-Math.log(1.0 - x) / logb);
}
return x + bi[k];
} */
/**
* Given the k b-ary digits of i in bdigits, returns the
* k digits of the integer radical inverse of i in idigits.
* This simply reverses the order of the digits.
*
* @param k number of digits in arrays
*
* @param bdigits digits in original order
*
* @param idigits digits in reverse order
*
*
*/
public static void reverseDigits (int k, int bdigits[], int idigits[]) {
for (int l = 0; l < k; l++)
idigits[l] = bdigits[k-l];
}
/**
* Computes the integer radical inverse of i in base b,
* equal to
* bkψb(i) if i has k b-ary digits.
*
* @param b base used for the operation
*
* @param i the value for which the integer radical inverse will be computed
*
* @return the integer radical inverse of i in base b
*
*/
public static int integerRadicalInverse (int b, int i) {
// Simply flips digits of i in base b.
int inverse;
for (inverse = 0; i > 0; i /= b)
inverse = inverse * b + (i % b);
return inverse;
}
/**
* Given the k digits of the integer radical inverse of i in bdigits,
* in base b, this method replaces them by the digits of the integer
* radical inverse of i + 1 and returns their number.
* The array must be large enough to hold this new number of digits.
*
* @param b base
*
* @param k initial number of digits in arrays
*
* @param idigits digits of integer radical inverse
*
* @return new number of digits in arrays
*
*/
public static int nextRadicalInverseDigits (int b, int k, int idigits[]) {
int l;
for (l = k-1; l >= 0; l--)
if (idigits[l] == b-1)
idigits[l] = 0;
else {
idigits[l]++;
return k;
}
if (l == 0) {
idigits[k] = 1;
return ++k;
}
return 0;
}
/**
* Computes the Faure permutation σb of the set
*
* {0,…, b - 1} and puts it in array pi.
*
* See the description in the introduction above.
*
* @param b the base
*
* @param pi an array of size at least b,
* to be filled with the permutation
*
*
*/
public static void getFaurePermutation (int b, int[] pi) {
// This is a recursive implementation.
// Perhaps not the most efficient...
int i;
if (b == 2) {
pi[0] = 0;
pi[1] = 1;
}
else if ((b & 1) != 0) {
// b is odd.
b--;
getFaurePermutation (b, pi);
for (i = 0; i < b; i++)
if (pi[i] >= b / 2)
pi[i]++;
for (i = b; i > b / 2; i--)
pi[i] = pi[i - 1];
pi[b / 2] = b / 2;
}
else {
b /= 2;
getFaurePermutation (b, pi);
for (i = 0; i < b; i++) {
pi[i] *= 2;
pi[i + b] = pi[i] + 1;
}
}
}
/**
* Computes the radical inverse of i in base b, where the digits
* are permuted using the permutation π.
* If
* i = ∑r=0k-1arbr, the method will compute and return
*
*
*
* x = ∑r=0k-1π[ar]b-r-1.
*
*
* @param b base b used for the operation
*
* @param pi an array of length at least b containing the permutation
* used during the computation
*
* @param i the value for which the radical inverse will be computed
*
* @return the radical inverse of i in base b
*
*/
public static double permutedRadicalInverse (int b, int[] pi, int i) {
double digit, radical, inverse;
digit = radical = 1.0 / (double) b;
for (inverse = 0.0; i > 0; i /= b) {
inverse += digit * (double) pi[i % b];
digit *= radical;
}
return inverse;
}
}
Ads help maintain this website.