sub factorial_gosper {
  # sqrt(2*n*PI + 1/3)*n**n/exp(n)
  my $n = shift;
  return exp($n*log($n) - $n + 0.5*log(2*$n*PI + 1.0/3.0));
}