Gamma関数を使おう
ちょい書きメモです。
ようじょはある日perlでΓ函数を使いたくなりました。
するとCPANにMath::GammaFunction*1がありました。
でもRMath.hが無いとかで入りませんでした。
ようじょは考えるのさえ面倒なので実装することにしました。
精度は欲しいですがPPの実装なので計算も軽そうなLanczos近似*2による実装をしてみました。
sub gamma($){ my $x = shift; my $u = 2.102394798991390; my $a_inf = 3.062185443705942e-1; my $a_0 = 1.024166094985555; my $a_1 = 4.258010456317367e-1; my $b_1 = 1.000008131602802; return (( ( $a_1 / ( $x + $b_1 ) + $a_1 ) / ( $x + $b_1 ) + $a_0) / $x + $a_inf ) * exp(( $x - 0.5 ) * log( $x + $u ) - $x ); }
まあいいかなと思ったら、C99にtgamma関数なるものがすでに実装されていたのでした。
ということでXSでラップするだけ
double gamma(input) double input CODE: RETVAL = tgamma( input ); OUTPUT: RETVAL
アホらし。
