FaureSequence.java

/*
 * Class:        FaureSequence
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.hups;
 
import umontreal.ssj.util.PrintfFormat;

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;

   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 < k; i++)
         numPoints *= b;
 
      // Compute the normalization factors.
      normFactor = 1.0 / Math.pow((double) b, (double) outDigits);
//      EpsilonHalf = 0.5*normFactor;
      double invb = 1.0 / b;
      factor = new double[outDigits];
      factor[0] = invb;
      for (j = 1; j < outDigits; j++)
         factor[j] = factor[j - 1] * invb;
 
      genMat = new int[dim * numCols][numRows];
      initGenMat();
   }

   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 };
 
}