Colby - 2 months ago 8
Pascal Question

# 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 gauss() : 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 gauss() : 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;
``````

the `rand()` is in `[0,1)` range ( `0 <= rand() < 1` )
so `2.0 * rand() - 1.0` is in `[-1,1)` range
so `x1` and `x2` are in `[-1,1)` range
so `w := x1 * x1 + x2 * x2` is in `[0,2)` range

and in `sqrt( -2.0 * ln( w ) / w )` the `w` is positive
so the natural logarithm: ln(w) should be negative
so w should be in `(0,1)` range
so that loop should not exit `until (w > 0.0)and (w < 1.0);`

working sample code (using SCAR Divi 3.41.00):

``````program New;

var hgauss : boolean;
var ngauss : double;

function gauss() : 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 > 0.0)and (w < 1.0);
w := sqrt( -2.0 * ln( w ) / w );
result := x1 * w;
ngauss := x2 * w;
hgauss := true;
end;
end;

begin
writeln( gauss() );
writeln( gauss() );
end.
``````
Source (Stackoverflow)