package umontreal.iro.lecuyer.hups;
/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */
import umontreal.iro.lecuyer.util.PrintfFormat;
import umontreal.iro.lecuyer.rng.RandomStream;


/**
 * This class implements point sets defined by integration
 * lattices of rank 1, defined as follows.
 * One selects an arbitrary positive integer n and a s-dimensional
 * integer vector 
 * (a0,..., as-1), where 
 * 0 <= aj < n for each j.
 * Usually, a0 = 1.  The points are defined by
 * 

 * 
 * 

 * ui = (i/n)(a0, a1,…, as-1) mod 1
 * 

 * for 
 * i = 0,..., n - 1.
 * These n points are distinct provided that n and the aj's have
 * no common factor.
 * 
 */
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.
   protected double[] shift;       // Random shift, initially null.




   /**
    * Instantiates a {@link Rank1Lattice} with n points and lattice
    *    vector a of dimension s.
    *  
    * @param n there are n points
    * 
    *    @param a the lattice vector
    * 
    *    @param s dimension of the lattice vector a
    * 
    */
   public Rank1Lattice (int n, int[] a, int s)  {
      dim = s;
      numPoints = n;
      normFactor = 1.0 / (double) n;
      v = new double[s];
      genAs = new int[s];
      for (int j = 0; j < s; j++) {
         if (a[j] < 0 || a[j] >= n)
            throw new IllegalArgumentException
               ("Rank1Lattice must have 0 <= a[j] < n");
         v[j] = normFactor * a[j];
         genAs[j] = a[j];
      }
   }


   /**
    * Returns the generator aj of the lattice. Its components
    *   are returned as a[j], for 
    * j = 0, 1,…,(s - 1).
    * 
    */
   public int[] getAs()  {
      return genAs;
   }


   /**
    * Adds a random shift to all the points of the point set,
    *   using stream stream to generate the random numbers.
    *   For each coordinate j from d1 to d2-1,
    *   the shift dj is generated uniformly over [0, 1) and added modulo 1 to
    *   all the coordinates of all the points.
    * 
    * @param d1 lower dimension of shift
    * 
    *    @param d2 upper dimension of shift is d2 - 1
    * 
    *    @param stream random number stream used to generate uniforms
    * 
    * 
    */
   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;
   }


   /**
    * Clears the random shift.
    * 
    */
   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;
      }
   }
}