Sample Output with uncommented line:
c = 10:         0.63331008207472
c = 16:         0.574732628084873
c = 27:         0.507598958744487
c = 46:         0.439468804563385
c = 77:         0.376031147677405
c = 129:        0.316746525404527
c = 359:        0.05 #rounding
c = 599:        1.1  #rounding
c = 1000:       0.04 #rounding

With commented line:
c = 10 pi = 1.6
c = 16 pi = 0.574732628084872737814588036475650211979002811242049849956301387578855115765977011850452313041695
c = 27 pi = 0.5
c = 46 pi = 0.44
c = 77 pi = 0.1
c = 129 pi = 0.32
c = 359 pi = 0.05


#invariant_dist.pl
use strict;
use Math::BigFloat;
use Math::Trig;

my @c = (10, 16, 27, 46, 77, 129, 359, 599, 1000);
my $rho = -5;
my @lambda;
for(0..@c-1) {
	$lambda[$_] = abs($c[$_]-($rho*sqrt($c[$_])));
} 
my $mu = 1;
my @pi = MMCQueueInvDisbn(\@lambda,$mu,\@c);
my $i = 0;
my $sum = 0;
my $mean = 0;

# Compute the invariant distribution for an M/M/c queue
sub MMCQueueInvDisbn { 
	my $lambda = shift;
	my $mu = shift;
	my $c = shift;
	my $rho;
	my $G = Math::BigFloat->new(0);
	my @pi;
	
	# Lets be efficent here!
	my @fact;
	foreach(0..$c[@c-1]) {
		$fact[$_] = 0;
	}
	
	for(my $i=0; $i < @c; $i++) {
		printf($i."\n");
		$rho = $lambda[$i] / $mu;
		$G = MMCQueueNormConst($lambda[$i],$mu,$c[$i],\@fact);
		for (0..$c[$i]-1) {
			$pi[$i][$_] = 0;
		}
		#for n > 250, n! = Inf so $pi = NaN so using BigFloat
		$pi[$i][$c[$i]] = $rho**$c[$i] / (Math::BigFloat->new(Factorial($c[$i], \@fact))->bstr()*$G->bstr());
		#Some Fucking Rounding is occuring
		#$pi[$i][$c[$i]] = (Math::BigFloat->new($rho)->bpow(Math::BigFloat->new($c[$i])))->bdiv(Math::BigFloat->new(Factorial($c[$i], \@fact))->bmul($G));
		printf("c = $c[$i]:\t $pi[$i][$c[$i]]\n"); 
	}
	
	return @pi;
}

# Compute the normalization constant for an M/M/c/c queue
sub MMCQueueNormConst { 
	my $lambda = shift;
	my $mu = shift;
	my $c = shift;
	my $fact = shift;
	
	my $rho = $lambda / $mu;
	my $G = Math::BigFloat->new(0);
	for (0..$c) {
		my $num = Math::BigFloat->new(Math::BigFloat->new($rho)->bpow(Math::BigFloat->new($_)));
		$G->badd($num->bdiv(Factorial($_, \@$fact)));
	}
	return $G; 
}

# Factorial subroutine
sub Factorial {
	my $n = shift;
	my $fact = shift;
	my $N = Math::BigFloat->new($n);
	
	if ($n == 0) { 
		return 1;
	} elsif($n > 300) {
		# Or could use Sterling's Approximation - trying BigFloat library first
		$N->bfac();
		return $N;
	} elsif (@$fact[$n] > 1) {
		return @$fact[$n]; 
	} elsif ($n < 0) { 
		die "$!\n";
	} else {
		$N = $N * Factorial($N-1, \@$fact);
		@$fact[$n] = $N;
		return $N;
	}
}