sub digamma { my ($x) = @_; my $result = 0; my($xx, $xx2, $xx4); #assert(x > 0); for ( ; $x < 7; ++$x) { $result -= 1/$x; } $x -= 1.0/2.0; $xx = 1.0/$x; $xx2 = $xx*$xx; $xx4 = $xx2*$xx2; $result += log($x)+(1.0/24.0)*$xx2-(7.0/960.0)*$xx4+(31.0/8064.0)*$xx4*$xx2-(127.0/30720.0)*$xx4*$xx4; return $result; }