Empty

2013年4月から社会人になりました。

SONY GOF FOR IT ~2問目: 実数の階乗~

Programming

概要

前回エントリーに引き続きGO FOR ITに参戦中。 今日は2問目に挑戦。

問題

2) 実数の階乗
ある検索サイトに5!と入力するとその計算結果である120が表示されます。
その検索サイトに2.5!と入力するとなんと3.32335097と表示されます。
さらにその検索サイトに(-1.9)!と入力すると-10.5705641と表示されます。
きっとそれらの仕組みはとても難しくて企業秘密に違いないので是非ともこれらを実行するプログラムを作ってほしい。
ただし、君のPCは古いのでネットワークや便利で高度な数学関数は入っていません。
入っている数学関数はsin,cos,tan,log,pow,floorなどの初歩的な関数のみです、残念ながら。

i)入力された整数a(0<=a<=10)の階乗を求めるプログラムを作ってください。
ii)入力された実数a(0<=a<=10)の階乗を求めるプログラムを作ってください。
iii)入力された実数a(-1.9<=a<=-1.1)の階乗を求めるプログラムを作ってください。

http://www.sony.co.jp/SonyInfo/Jobs/newgrads/sus/q02.html

i) の整数を対象にした問題について。

def fact(n):
    """普通にfactorial."""
    ANSWERS = {0:1, 1:1}
    def _fact(_n):
        if not _n in ANSWERS.keys():
            ANSWERS[_n] =  _n * _fact(_n -1)
        return ANSWERS[_n]
    return _fact(n)

特に説明とかいらないと思うので割愛。

ii), iii) について。

ここからが、本題。この問いはどうやら、基本的な数学関数のみでガンマ関数を実装する問題に落ち着くみたい。実装する前に、ガンマ関数についてあまり知識がなかったのでいろいろと調べてみた。

基礎知識
すごい参考になったページ

上記のページを見ていたら、

日本語Wikipediaにはスターリングの近似なる項目があり、ここにある式は比較的簡単ながらそこそこの近似値が得られるようになっています。ただ、近似値の精度は（方法にもよるとはいえ）Excelのgammalnに劣り、納得できる水準ではありません。
　そこで使うのがLanczos approximation（ランチョス近似）なのですが、...

http://void.yamicha.com/notebook/tdist_gamma.htm

という近似手法を紹介してくれた。スターリングの近似はシンプルだけど、精度が 悪いらしい 若干他の手法に比べると劣る部分があるらしい。

シンプルさの参考に

ということで、Lanczos approximation について調べてみたら、ご丁寧なことにPythonコードが載っていて、おかげさまで、すっかり実装する気がなくなってしまった。

Simple implementation

The following implementation in the Python programming language works for complex arguments and typically gives 15 correct decimal places:

from cmath import *
 
# Coefficients used by the GNU Scientific Library
g = 7
p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
     771.32342877765313, -176.61502916214059, 12.507343278686905,
     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
 
def gamma(z):
    z = complex(z)
    # Reflection formula
    if z.real < 0.5:
        return pi / (sin(pi*z)*gamma(1-z))
    else:
        z -= 1
        x = p[0]
        for i in range(1, g+2):
            x += p[i]/(z+i)
        t = z + g + 0.5
        return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x
http://en.wikipedia.org/wiki/Lanczos_approximation

基本的に、よくわからないが、配列pの値を使ってうまーく近似しているんだと。その配列pはどうやって得られたのだろうか・・・。むしろあの配列pの値を求めるのが、本筋なのではなかろうかと思ったが、今日はここまでとした。その他にも、Spouge's approximationっていう近似手法もあるらしい。

WEBで見つけたコード

Wikipediaに載ってたアルゴリズムの精度確認

でも、ちゃんとした結果になるか、気になったので確認だけしてみた。

こんな入力(Wikipediaに答えが載っていたもの)
if __name__ == "__main__":
    print "gamma(-3/2): %s" % (gamma(-3/2),)
    print "gamma(-1/2): %s" % (gamma(-1/2),)
    print "gamma(1):    %s" % (gamma(1),)
    print "gamma(3/2):  %s" % (gamma(3/2),)
    print "gamma(2):    %s" % (gamma(2),)
    print "gamma(5/2):  %s" % (gamma(5/2),)
    print "gamma(3):    %s" % (gamma(3),)
    print "gamma(7/2):  %s" % (gamma(7/2),)
    print "gamma(4):    %s" % (gamma(4),)
そして出力
gamma(-3/2):    (6.413262697e+15+0j)
gamma(-1/2):    (-2.5653050788e+16-0j)
gamma(1):       (1+0j)
gamma(3/2):     (1+0j)
gamma(2):       (1+0j)
gamma(5/2):     (1+0j)
gamma(3):       (2+0j)
gamma(7/2):     (2+0j)
gamma(4):       (6+0j)

あれれー？なんかおかしいぞー。この結果は、見なかったことに...はい、このエントリーのこと忘れました。


ということで、明日元気だったら、3問目に挑戦します。