//
//  Ran0.java
//  Random deviate generator for the Huygens Landscapes
//
//  Created by Cara MacNish on 1/04/05.
//  Copyright 2005-2007 CSSE, UWA. All rights reserved.
//

public class Ran0 {
  //private final static long IA = 69621l;     // long ensures mults are done long 
  //IQ = 30845;
  //IR = 23902;
  
  final static int IM = 2147483647;           // modulus, m = 2^31 - 1
  private long IA;
  private double AM = 1.0/IM;
  private boolean gaussAvail = false;
  private double gaussVal = 0;
  private int idum;
  
  public Ran0 (long IA, long seed) {
    this.IA = IA;
    // clean seed - map into range [1..2^31-2]
    seed = Math.abs(seed);
    // if (seed == 0) System.out.println("Warning: zero seed moved");
    seed = ( (seed-1) % (IM-1) ) +1;   // map into [1..2^31-2]
        
    // scatter seed
    idum = (int)((IA*seed) % IM);    // one multiple by 'a'
    StringBuffer s = new StringBuffer(Integer.toHexString(idum));     // convert to hex
    s.reverse();                                    // flip so least sig becomes most sig
    while (s.length() < 8) s.append('0');           // pad out to 8 digits
    long temp = Long.parseLong(s.toString(),16);    // convert back to decimal
    idum = (int) ((temp % (IM-1)) + 1);             // may be larger than IM - map to [1, 2^31-1]
  }
  
  public double uniform () {
    idum = (int) ((IA * idum) % IM);
    return AM * idum;
  }
  
  public double gauss () {
    double gdev=0;
    if (!gaussAvail) {
      while (!gaussAvail) {
        double v1 = 2*uniform()-1;
        double v2 = 2*uniform()-1;
        double rsq = v1*v1 + v2*v2;
        if (rsq < 1.0 && rsq != 0) {
          double fac = Math.sqrt(-2*Math.log(rsq)/rsq);
          gaussVal = v1 * fac;
          gaussAvail = true;
          gdev = v2 * fac;
        }
      }
    }
    else {
      gaussAvail = false;
      gdev = gaussVal;
    }
    return gdev;
  }
  
  public double gauss (double mean, double stdev) {
     double gdev=0;
    if (!gaussAvail) {
      while (!gaussAvail) {
        double v1 = 2*uniform()-1;
        double v2 = 2*uniform()-1;
        double rsq = v1*v1 + v2*v2;
        if (rsq < 1.0 && rsq != 0) {
          double fac = Math.sqrt(-2*Math.log(rsq)/rsq);
          gaussVal = v1 * fac;
          gaussAvail = true;
          gdev = v2 * fac;
        }
      }
    }
    else {
      gaussAvail = false;
      gdev = gaussVal;
    }
    return mean + stdev * gdev;
  }
  
}