Gaussian in pascal

Question

Я попытался портировать код непосредственно из исходного кода java на pascal, однако он бросает ошибку времени выполнения.

Как я могу получить правильную гауссову кривую? Как насчет паскалей, встроенных в функции?

Исходный код:

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; 
    } 
} 

Первая попытка паскаль порт (броски ошибки во время выполнения):

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; 

Недопустимая операция с плавающей точкой здесь:

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

Вторая попытка преобразования (работает, но я не уверен, что математика верна):

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; 

Answer 1

rand() находится в [0,1) диапазоне (0 <= rand() < 1)
так 2.0 * rand() - 1.0 находится в [-1,1) диапазоне
так x1 и x2 находятся в [-1,1) диапазоне
так w := x1 * x1 + x2 * x2 находится в [0,2) диапазоне

и в sqrt(-2.0 * ln(w)/w)w является положительным
поэтому натуральный логарифм: ln (w) Должен быть отрицательным
так ж должна быть в диапазоне (0,1)
так, чтобы петля не должна выйти until (w > 0.0)and (w < 1.0);

работает образец кода (с использованием SCAR Дива 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. 

Answer 2

Ваш порт от Java Паскаля неисправна для важной части

do {...} while (s >= 1 || s == 0);

должны быть переведены на

repeat {...} until ((s<1) and (s<>0));

Так вы ошиблись условие завершения. Java завершает цикл, если 0 < s < 1, но ваш цикл завершен, если w >= 1.

Если w > 1 у вас есть -2*ln(w) < 0, а исключение с плавающей запятой происходит от квадратного корня отрицательного числа!

И для большинства Pascal версии вашего именования стандартных функций необычно, ИМО следует читать

repeat 
    x1 := 2.0 * random - 1.0; 
    x2 := 2.0 * random - 1.0; 
    w := x1 * x1 + x2 * x2; 
until (w<1.0) and (w>0.0);   

w := sqrt(-2.0*ln(w)/w); 
result := x1 * w; 
ngauss := x2 * w; 

И заметьте, что вы действительно должны использовать ln не ваш самодельный основанию 10 логарифма log. Используемый метод Marsaglia's polar method.