2016年9月19日月曜日

160919

Ruby


Bernoulli number(2)

次の3 通りで求めてみた。
①定義とBn = 0 (n が3 以上の奇数)から導く方法
②「ベルヌーイ数とゼータ関数」の第2章 2.2 に載っているアルゴリズム
③「ベルヌーイ数とゼータ関数」の命題1.15 を利用する方法

def c(n, r)
  return 1 if r == 0
  return c(n, n - r) if r > n - r
  (n - r + 1..n).inject(:*) / (1..r).inject(:*)
end

def bernoulli_0(n)
  return 1r / 2 if n == 1
  return 0 if n % 2 == 1
  a = [1r]
  (1..n / 2).each{|i| a << 1r / 2 - (0..i - 1).inject(0){|s, j| s + c(2 * i + 1, 2 * j) * a[j]} / (2 * i + 1)}
  a[-1]
end

def bernoulli_1(n)
  return 1r / 2 if n == 1
  return 0 if n % 2 == 1
  a = []
  (0..n).each{|i|
    a << 1r / (i + 1)
    i.downto(1){|j| a[j - 1] = j * (a[j - 1] - a[j])}
  }
  a[0]
end

def bernoulli_2(n)
  return 1r / 2 if n == 1
  return 0 if n % 2 == 1
  a = [1, 1r / 6]
  (2..n / 2).each{|i| a << - (1..i - 1).inject(0){|s, j| s + c(2 * i, 2 * j) * a[i - j] * a[j]} / (2 * i + 1)}
  a[n / 2]
end

n = 100
p ary0 = (0..n).map{|i| bernoulli_0(i)}
p ary0 == (0..n).map{|i| bernoulli_1(i)}
p ary0 == (0..n).map{|i| bernoulli_2(i)}

出力結果
[(1/1), (1/2), (1/6), 0, (-1/30), 0, (1/42), 0, (-1/30), 0, (5/66), 0, (-691/2730), 0, (7/6), 0, (-3617/510), 0, (43867/798), 0, (-174611/330), 0, (854513/138), 0, (-236364091/2730), 0, (8553103/6), 0, (-23749461029/870), 0, (8615841276005/14322), 0, (-7709321041217/510), 0, (2577687858367/6), 0, (-26315271553053477373/1919190), 0, (2929993913841559/6), 0, (-261082718496449122051/13530), 0, (1520097643918070802691/1806), 0, (-27833269579301024235023/690), 0, (596451111593912163277961/282), 0, (-5609403368997817686249127547/46410), 0, (495057205241079648212477525/66), 0, (-801165718135489957347924991853/1590), 0, (29149963634884862421418123812691/798), 0, (-2479392929313226753685415739663229/870), 0, (84483613348880041862046775994036021/354), 0, (-1215233140483755572040304994079820246041491/56786730), 0, (12300585434086858541953039857403386151/6), 0, (-106783830147866529886385444979142647942017/510), 0, (1472600022126335654051619428551932342241899101/64722), 0, (-78773130858718728141909149208474606244347001/30), 0, (1505381347333367003803076567377857208511438160235/4686), 0, (-5827954961669944110438277244641067365282488301844260429/140100870), 0, (34152417289221168014330073731472635186688307783087/6), 0, (-24655088825935372707687196040585199904365267828865801/30), 0, (414846365575400828295179035549542073492199375372400483487/3318), 0, (-4603784299479457646935574969019046849794257872751288919656867/230010), 0, (1677014149185145836823154509786269900207736027570253414881613/498), 0, (-2024576195935290360231131160111731009989917391198090877281083932477/3404310), 0, (660714619417678653573847847426261496277830686653388931761996983/6), 0, (-1311426488674017507995511424019311843345750275572028644296919890574047/61410), 0, (1179057279021082799884123351249215083775254949669647116231545215727922535/272118), 0, (-1295585948207537527989427828538576749659341483719435143023316326829946247/1410), 0, (1220813806579744469607301679413201203958508415202696621436215105284649447/6), 0, (-211600449597266513097597728109824233673043954389060234150638733420050668349987259/4501770), 0, (67908260672905495624051117546403605607342195728504487509073961249992947058239/6), 0, (-94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330)]
true
true

0 件のコメント:

コメントを投稿

注: コメントを投稿できるのは、このブログのメンバーだけです。