SearcherKorobov.java

/*
 * Class:        SearcherKorobov
 * Description:  searches the best Korobov lattices with respect to a 
                 given discrepancy
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Richard Simard
 * @since        January 2009
 
 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.
 
 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 
 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.discrepancy;
 
import umontreal.ssj.util.Num;
import umontreal.ssj.hups.KorobovLattice;

public class SearcherKorobov extends Searcher {
   protected int bestA; // best generator a for Korobov lattices
 
   protected void print(int y) {
      // For debugging
      System.out.printf("  a = %d%n", y);
   }
 
   private void calcAs(int n, int s, int a) {
      long b = a;
      bestAs[0] = 1;
      for (int j = 1; j < s; ++j)
         bestAs[j] = (int) ((b * bestAs[j - 1]) % n);
   }
 
   private double exhaust(int s, boolean relPrime) {
      int n = disc.getNumPoints();
      double err;
      gamma = disc.getGamma();
      bestVal = Double.MAX_VALUE;
      bestA = -1;
 
      for (int i = 2; i < n; i++) {
         if (relPrime) {
            // Consider only values i relatively prime to n
            if (power2F) {
               if ((i & 1) == 0)
                  continue;
            } else {
               if (Num.gcd(n, i) != 1)
                  continue;
            }
         }
         // print (i);
         lat = new KorobovLattice(n, i, s);
         err = disc.compute(lat, gamma);
         if (err < bestVal) {
            bestVal = err;
            bestA = i;
         }
      }
 
      calcAs(n, s, bestA);
      return bestVal;
   }
 
   private double random(int s, int k, boolean relPrime) {
      int n = disc.getNumPoints();
      if (k >= n)
         return exhaust(s, relPrime);
 
      gamma = disc.getGamma();
      final int nm1 = n - 1;
      double err;
      bestVal = Double.MAX_VALUE;
      bestA = -1;
      int a;
      int i = 0;
 
      while (i < k) {
         if (power2F) {
            a = gen.nextInt(2, nm1);
            if (relPrime)
               a |= 1;
         } else {
            do {
               a = gen.nextInt(2, nm1);
            } while (relPrime && (Num.gcd(n, a) != 1));
         }
         // print (a);
         lat = new KorobovLattice(n, a, s);
         i++;
         err = disc.compute(lat, gamma);
         if (err < bestVal) {
            bestVal = err;
            bestA = a;
         }
      }
 
      calcAs(n, s, bestA);
      return bestVal;
   }

   public SearcherKorobov(Discrepancy disc, boolean primeN) {
      super(disc, primeN);
   }

   public double exhaust(int s) {
      return exhaust(s, false);
   }

   public double exhaustPrime(int s) {
      if (primeN)
         return exhaust(s, false); // all a are rel. prime to n
      else
         return exhaust(s, true);
   }

   public double random(int s, int k) {
      return random(s, k, false);
   }

   public double randomPrime(int s, int k) {
      if (primeN)
         return random(s, k, false);
      else
         return random(s, k, true);
   }

   public int getBestA() {
      return bestA;
   }
 
}