#!/usr/bin/perl
#
# CS661, Homework 4, problem 3
# Author: Dmitri Tikhonov <dxt22@njit.edu>
# Date: February 22, 2009

use strict;
use warnings;

use constant DEFAULT_NUMBER_OF_ITEMS => 40;

use constant {
    A => 7141,
    C => 54773,
    M => 259200,
};

sub get_lcg {
    my $seed = shift;
    return sub {
        $seed = (A * $seed + C) % M;
        return $seed / M;
    }
}

sub sample_mean {
    my $sum;
    $sum += $_ for (@_);
    return $sum / @_;
}

sub sample_variance {
    my $mean = shift;
    my $sum;
    $sum += ($_ - $mean) ** 2 for (@_);
    return $sum / (@_ - 1);
}

sub usage {
    print <<USAGE;
Usage: $0 --part X arguments

Description: Performs parts (b), (c), and (d) of the homework, based on
the value of X.

Arguments:
    Part (b):   Seed (optional)
    Part (c):   Seed (optional)
    Part (d):   Seed (optional)

Default seed is zero
USAGE
    exit shift;
}

use Getopt::Long;
GetOptions(
    "part=s"    => \my $part,
    "help"      => sub { usage(0) },
) or usage(1);

my $seed = shift @ARGV || 0;
die "Seed must be between 0 and ${[M]}" unless $seed >= 0 && $seed < M;
my $lcg = get_lcg($seed);

if ('b' eq lc($part)) {
    # Part (b): just print the first forty numbers

    print $lcg->(), "\n" for (1 .. DEFAULT_NUMBER_OF_ITEMS);
} elsif ('c' eq lc($part)) {
    # Part (c): split the first forty numbers into bins

    my @bins;
    for (1 .. DEFAULT_NUMBER_OF_ITEMS) {
        my $x = $lcg->();
        if ($x < 0.2) {
            ++$bins[0];
        } elsif ($x < 0.5) {
            ++$bins[1];
        } else {
            ++$bins[2];
        }
    }

    for (my $i = 0; $i < @bins; ++$i) {
        printf("BIN %d: %2d items (%.2f%%)\n", $i, $bins[$i],
                100 * $bins[$i] / DEFAULT_NUMBER_OF_ITEMS);
    }
} elsif ('d' eq lc($part)) {
    # Part (d): calculate sample mean and sample variances for sequences
    # of size 10, 100, 1000, and 10000

    my @sequence_sizes = (10, 100, 1000, 10_000);
    my @sequence = map { $lcg->() } (1 .. $sequence_sizes[-1]);
    for my $size (@sequence_sizes) {
        my $mean = sample_mean(@sequence[0 .. $size - 1]);
        my $variance = sample_variance($mean, @sequence[0 .. $size - 1]);
        printf("SEQ of size %9d: sample mean %.4f, sample variance %.4f\n",
            $size, $mean, $variance);
    }
} else {
    die "Invalid part `$part'";
}