package Gaussian;
require Exporter;

@ISA = qw(Exporter);
@EXPORT = qw(gauss_rand gauss);

my $phase=0;
my $next = 0;

sub gauss_rand {
    # Knuth, Section 3.4.1 
    my $v1, $v2, $s;
    if($phase==1){
	$phase = 0;
	return $next;
    }
    do {
	$v1 = 2*rand()-1;
	$v2 = 2*rand()-1;
	$s = ($v1*$v1) + ($v2*$v2);
    } while ($s >= 1);
    my $mult = sqrt(-2*log($s)/$s);
    $next = $v2*$mult;
    $phase = 1;
    return ($v1*$mult);
}

sub gauss {
	my $v = shift;
	my $r = shift;
	return exp(-1*($r*$r)/($v*$v));
}

1;