Bernoulli数（1）

ベルヌイ数、またはベルヌーイ数。
http://ja.wikipedia.org/wiki/%E3%83%99%E3%83%AB%E3%83%8C%E3%83%BC%E3%82%A4%E6%95%B0

定義は2通りあるらしいが、Wikipediaと違う定義を使う。

右辺の分母を両辺にかけて整理すると、

これをきのう作った有理数クラスを使って、ベルヌイ数を1000まで計算したところ、1分ほどかかった。

B₀ = 1
B₁ = 1/2
B₂ = 1/6
B₃ = 0
B₄ = -1/30
B₅ = 0
B₆ = 1/42
B₈ = -1/30
B₁₀ = 5/66
B₁₂ = -691/2730
B₁₄ = 7/6
B₁₆ = -3617/510
B₁₈ = 43867/798
B₂₀ = -174611/330
B₂₂ = 854513/138
B₂₄ = -236364091/2730
B₂₆ = 8553103/6
B₂₈ = -23749461029/870
B₃₀ = 8615841276005/14322
B₃₂ = -7709321041217/510
B₃₄ = 2577687858367/6
B₃₆ = -26315271553053477373/1919190
B₃₈ = 2929993913841559/6
B₄₀ = -261082718496449122051/13530
B₄₂ = 1520097643918070802691/1806
B₄₄ = -27833269579301024235023/690
B₄₆ = 596451111593912163277961/282
B₄₈ = -5609403368997817686249127547/46410
B₅₀ = 495057205241079648212477525/66
B₁₀₀₀ = -18243104738661887254572640256857788
8793383368670429060521971581576411262572624911158657472
5773210697096154899246274955229080874882995394551887918
5675822415516684926972441849140122425798309556170986299
2465225174097919156372263614283427805489710022810454653
0844116137235069692022011624417917606802626020196202602
5579005841653927133285280600096662846763906834342263807
0295122610811666617281581715702361188930366816683991915
6379768387784569011484312275342742688059179988378025533
8278664578660218504589596267044201144363032146025948676
4674312436994856054301765557425137115021340105105840867
9874766352952749178734973676859834707623881634625147148
9942512878190574323531299070406930309477389251738705417
6806531183648189451892725726445949589759600705334767585
3897699248576309729639976364832442643512622073858780110
7315398330998175557751360081111707976250597322951308884
9006701133391676419537939945123776103061984293109331214
6321416835426077466412320898548150646291295965369973806
0825642880197849098973016582688092035550306928461519170
6946560725764114918719765109055159668403124118455436505
9302140284922169134185281979123358930199410122917734417
9402749357465188105943227449435409223195489428074206847
2714619294213343605461147540486788631325011439968153275
3236429290625909341100039136833631213891562170153595481
4084208794241665492294270773347605587841576592758201421
4726584822236443691314366097570085473354584000998591519
0584047337934331297339403392719579093995842312746836871
1696749786460913411872527166990047126222109345933847358
9242309517183798837432563465487604170316077418754242710
0652698181905912716906954466338361203745255515267088218
3839963301642034037323653333521203382720213197180035994
2204589948764600183502703856341178077687451616229338340
631455056219106004731529642292049578901/342999030

from rational import cRational
def calcBernoulli(m):

    B = [ cRational(0) for x in range(m) ]

    for n in range(m):

        C = 1

        for k in range(n):

            if k == 0:

                C = 1

            else:

                C *= n + 2 - k

                C /= k

            B[n] += C * B[k]

        B[n] = 1 - B[n] / (n + 1)

    return B
n = 1001

B = calcBernoulli(n)

for k in range(n):

    print "B<sub>" + str(k) + "</sub> = " + str(B[k])