gamma function

$pi = 3.14159265358979323846;
$e = 2.718281828459045;

while (<>) {
  chomp;
  print "gamma($_) (aka (", $_+1. ")!: ", &gamma($_), "\n";
}
    
sub gamma {
  my $z = shift @_;
  my $x = sqrt(2*$pi/$z);
  my $sqbit = $z * sinh(1/$z + 1/(810 *($z**6)));
  return 2 * $x * (($z/$e * sqrt($sqbit)) ** $z);
}

sub sinh {
  $z = shift @_;
  return .5 * (exp($z) - exp(-$z));
}