Rank1Lattice.java

/*
 * Class:        Rank1Lattice
 * Description:  Rank-1 lattice
 * 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;
import umontreal.ssj.rng.RandomStream;

public class Rank1Lattice extends PointSet {
 
   protected int[] genAs; // Lattice generator: a[i]
   protected double[] v; // Lattice vector: v[i] = a[i]/n
   protected double normFactor; // 1/n.
 
   private void initN(int n) {
      numPoints = n;
      normFactor = 1.0 / (double) n;
      for (int j = 0; j < dim; j++) {
         int amod = (genAs[j] % n) + (genAs[j] < 0 ? n : 0);
         v[j] = normFactor * amod;
      }
   }

   public Rank1Lattice(int n, int[] a, int s) {
      dim = s;
      v = new double[s];
      genAs = new int[s];
      for (int j = 0; j < s; j++) {
         genAs[j] = a[j];
      }
      initN(n);
   }

   public void setNumPoints(int n) {
      initN(n);
   }

   public int[] getAs() {
      return genAs;
   }

   public void addRandomShift(int d1, int d2, RandomStream stream) {
      if (null == stream)
         throw new IllegalArgumentException(PrintfFormat.NEWLINE + "   Calling addRandomShift with null stream");
      if (0 == d2)
         d2 = Math.max(1, dim);
      if (shift == null) {
         shift = new double[d2];
         capacityShift = d2;
      } else if (d2 > capacityShift) {
         int d3 = Math.max(4, capacityShift);
         while (d2 > d3)
            d3 *= 2;
         double[] temp = new double[d3];
         capacityShift = d3;
         for (int i = 0; i < d1; i++)
            temp[i] = shift[i];
         shift = temp;
      }
      dimShift = d2;
      for (int i = d1; i < d2; i++)
         shift[i] = stream.nextDouble();
      shiftStream = stream;
   }

   public void clearRandomShift() {
      super.clearRandomShift();
      shift = null;
   }
 
   public String toString() {
      StringBuffer sb = new StringBuffer("Rank1Lattice:" + PrintfFormat.NEWLINE);
      sb.append(super.toString());
      return sb.toString();
   }
 
   public double getCoordinate(int i, int j) {
      double x = (v[j] * i) % 1.0;
      if (shift != null) {
         if (j >= dimShift) // Extend the shift.
            addRandomShift(dimShift, j + 1, shiftStream);
         x += shift[j];
         if (x >= 1.0)
            x -= 1.0;
         if (x <= 0.0)
            x = EpsilonHalf; // avoid x = 0
      }
      return x;
   }
 
   // Recursive method that computes a^e mod m.
   protected long modPower(long a, int e, int m) {
      // If parameters a and m == numPoints could be omitted, then
      // the routine would run much faster due to reduced stack usage.
      // Note that a can be larger than m, e.g. in lattice sequences !
 
      if (e == 0)
         return 1;
      else if (e == 1)
         return a % m;
      else if ((e & 1) == 1)
         return (a * modPower(a, e - 1, m)) % m;
      else {
         long p = modPower(a, e / 2, m);
         return (p * p) % m;
      }
   }
 
   protected double radicalInverse(int base, int i) {
      double digit, radical, inverse;
      digit = radical = 1.0 / (double) base;
      for (inverse = 0.0; i > 0; i /= base) {
         inverse += digit * (double) (i % base);
         digit *= radical;
      }
      return inverse;
   }
 
   public PointSetIterator iterator() {
      return new Rank1LatticeIterator();
   }
 
// ************************************************************************
 
   protected class Rank1LatticeIterator extends PointSet.DefaultPointSetIterator {
      public double nextCoordinate() {
         // I tried with long's and with double's. The double version is
         // 4.5 times faster than the long version.
         if (curPointIndex >= numPoints || curCoordIndex >= dim)
            outOfBounds();
//      return (curPointIndex * v[curCoordIndex++]) % 1.0;
         double x = (curPointIndex * v[curCoordIndex]) % 1.0;
         if (shift != null) {
            if (curCoordIndex >= dimShift) // Extend the shift.
               addRandomShift(dimShift, curCoordIndex + 1, shiftStream);
            x += shift[curCoordIndex];
            if (x >= 1.0)
               x -= 1.0;
            if (x <= 0.0)
               x = EpsilonHalf; // avoid x = 0
         }
         curCoordIndex++;
         return x;
      }
   }
}