# Test gaussian()

use strict;
use warnings;

use Test::More 'no_plan';

my $CNT = 500000;

use_ok('Math::Random::MT::Auto', qw/gaussian/);
can_ok('Math::Random::MT::Auto', qw/gaussian/);
can_ok('main', qw/gaussian/);

# Cumulative distribution function
sub cdf
{
    my $z = $_[0];

    my $sum = 1.0 * $z;
    my $k = 1;
    my $factK = 1.0;

    my $term;
    do {
        $term = ($z**($k+$k+1)) / (($k+$k+1) * (2.0**$k) * $factK);
        $sum -= $term;
        $k++;
        $factK *= $k;

        $term = ($z**($k+$k+1)) / (($k+$k+1) * (2.0**$k) * $factK);
        $sum += $term;
        $k++;
        $factK *= $k;
    } while (abs($term) > 1e-50);

    $sum *= .398942280401432678;
    $sum += 0.5;

    return ($sum);
}

my %df;
for my $ii (-40 .. 39) {
    $df{$ii} = cdf(($ii+1)/10) * $CNT;
}
for (my $ii=39; $ii > -40; $ii--) {
    $df{$ii} -= $df{$ii-1};
}


### - Gaussian Function - ###

my (%bell, $x, $dev);

# Get random numbers and put them in bins
my $loops = 3;
LOOP1:
undef(%bell);
for (1 .. $CNT) {
    eval { $x = gaussian(10); };
    if ($@) {
        fail('gaussian(10) failed: ' . $@);
        exit(1);
    }

    # Handle 'rounding' using int()
    if ($x < 0) {
        $x = int($x) - 1;
    } else {
        $x = int($x);
    }

    # Make sure the tails don't overflow
    if ($x > 39) {
        $x = 40;
    } elsif ($x < -39) {
        $x = -40;
    }

    $bell{$x}++;
}

$dev = 0;
for my $ii (-40 .. 39) {
    my $bar1 = 3 * sqrt($df{$ii});
    my $bar2 = .025 * $df{$ii};
    my $bar = ($bar1 < $bar2) ? $bar2 : $bar1;
    if (($bell{$ii} < ($df{$ii}-$bar)) || (($df{$ii}+$bar) < $bell{$ii})) {
        $dev++;
    }
}

if (($dev > 3) && $loops--) {
    goto LOOP1;
}
ok($dev <= 3, 'Looks like a bell curve');


### - Gaussian OO - ###

my $prng;
eval { $prng = Math::Random::MT::Auto->new(); };
if (! ok(! $@, '->new() worked')) {
    diag('->new() died: ' . $@);
}
can_ok($prng, qw/rand irand gaussian exponential erlang poisson binomial
                 shuffle srand get_seed set_seed get_state set_state/);

# Get random numbers and put them in bins
$loops = 3;
LOOP2:
undef(%bell);
for (1 .. $CNT) {
    eval { $x = $prng->gaussian(10); };
    if ($@) {
        fail('$prng->gaussian(10) failed: ' . $@);
        exit(1);
    }

    # Handle 'rounding' using int()
    if ($x < 0) {
        $x = int($x) - 1;
    } else {
        $x = int($x);
    }

    # Make sure the tails don't overflow
    if ($x > 39) {
        $x = 40;
    } elsif ($x < -39) {
        $x = -40;
    }

    $bell{$x}++;
}

$dev=0;
for my $ii (-40 .. 39) {
    my $bar1 = 3 * sqrt($df{$ii});
    my $bar2 = .025 * $df{$ii};
    my $bar = ($bar1 < $bar2) ? $bar2 : $bar1;
    if (($bell{$ii} < ($df{$ii}-$bar)) || (($df{$ii}+$bar) < $bell{$ii})) {
        $dev++;
    }
}

if (($dev > 3) && $loops--) {
    goto LOOP2;
}
ok($dev <= 3, 'Looks like a bell curve');

exit(0);

# EOF