mercredi 13 juillet 2016

Gaussian in pascal

I tried porting the code directly from the java source code to pascal, however it is throwing a run time error.

How can I get a proper Gaussian curve? What about pascals built in functions?

Original source code

    synchronized public double nextGaussian() {
    // See Knuth, ACP, Section 3.4.1 Algorithm C.
    if (haveNextNextGaussian) {
        haveNextNextGaussian = false;
        return nextNextGaussian;
    } else {
        double v1, v2, s;
        do {
            v1 = 2 * nextDouble() - 1; // between -1 and 1
            v2 = 2 * nextDouble() - 1; // between -1 and 1
            s = v1 * v1 + v2 * v2;
        } while (s >= 1 || s == 0);
        double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
        nextNextGaussian = v2 * multiplier;
        haveNextNextGaussian = true;
        return v1 * multiplier;
    }
}

First attempt at pascal port: (throws runtime error)

  function log (n : double) : double; 
  begin 
    result := ln(n) / ln(10); 
  end; 

  var hgauss : boolean;
  var ngauss : double;

  function guass() : double;
  var x1, x2, w : double;
  begin
    if hgauss then
    begin
      result := ngauss;
      hgauss := false;
    end else
    begin
      repeat
        x1 := 2.0 * rand() - 1.0;
        x2 := 2.0 * rand() - 1.0;
        w := x1 * x1 + x2 * x2;
      until w >= 1.0;         

      w := sqrt( (-2.0 * log( w ) ) / w );
      result := x1 * w;
      ngauss := x2 * w;
      hgauss := true;
    end;
  end; 

Invalid floating point operation here:

w := sqrt((-2.0 * log( w ) ) / w);

Second attempt at conversion (runs but I am not sure the math is correct)

  function log (n : double) : double; 
  begin 
    result := ln(n) / ln(10); 
  end; 

  var hgauss : boolean;
  var ngauss : double;

  function guass() : double;
  var x1, x2, w, num : double;
  begin
    if hgauss then
    begin
      result := ngauss;
      hgauss := false;
    end else
    begin
      repeat
        x1 := 2.0 * rand() - 1.0;
        x2 := 2.0 * rand() - 1.0;
        w := x1 * x1 + x2 * x2;
      until w >= 1.0;         

      num := -2.0 * log( w )  / w;
      w := sqrt(abs(num));
      if num < 0 then w := -w;
      result := x1 * w;
      ngauss := x2 * w;
      hgauss := true;
    end;
  end;




Aucun commentaire:

Enregistrer un commentaire