#!/usr/bin/perl 

########################################
#Author: Jesse Lovegren 29-May-2011
#This script generates Rom critical values
#for a given alpha and number of observations
#based on equations given in Biometrika 77(3):663-5
#
#
# Usage: rom.pl [n] [alpha]
#
########################################

sub fac{
	my $out = 1;
	my $i = shift @_;
	for $j (1 .. $i){
		$out *= $j;
		}
	return $out;
	}
sub combine{
	my $n = shift @_;
	my $k = shift @_;
	my $out;
	$out = &fac($n) / (&fac($n-$k) * &fac($k));
	return $out;
}	

my $n = shift @ARGV;
my $alpha = shift @ARGV;

unless ($n, $alpha){ die ("Usage: rom.pl [n] [alpha]\n");}
##############################
print "="x50, "\nRom critical values, n=$n, alpha=$alpha\n" , "="x50, "\n";
###########################
my @powers;

$powers[1] = $alpha;
for $i ( 2 .. $n ){
	$powers[$i] = $powers[$i-1] + $alpha ** $i;
	}
my $xterm = $powers[$n-1];
###########################
my @rom;
my $yterm = 0;;
###########################
###########################
MAINLOOP : for $p ( 1 .. $n ) {
		if ($p == 1)	{	$rom[$p] = $alpha;
					next;
				}
		elsif ($p == 2)	{	$rom[$p] = $alpha/2;
					next;
				}
		else	{	#the calculable cases
			YLOOP : for $q ( 3 .. $p ){
						$yterm += &combine($n,$q-2) * ($rom[$q-1]) ** ($n-$q+2);
							}
			$rom[$p] = ($powers[$p-1] - $yterm) / $p;
			}
		}
##########################
PRINTLOOP : for $r ( 1 .. $#rom ){
	print "a($r,$alpha) = $rom[$r]\n";
		}

print "=" x 50 , "\n";