Akiyama-Tanigawa 法によるベルヌーイ数の計算

ベルヌーイ数

毎日暑くてやる気もしないので、ネットをプラプラ眺めていたら、ベルヌーイ数を求めるパスカルの三角形的な面白い式を見つけたので計算してみました。

オープンアクセスのようで、誰でも見れます。
https://cs.uwaterloo.ca/journals/JIS/vol3.html
Article 00.2.9:     Masanobu Kaneko, "The Akiyama-Tanigawa algorithm for Bernoulli numbers"

計算は中身を読みもせず、図を見てそのままプログラムにします。

漸化式 $b_{i,j}=i\cdot (b_{i, j-1}-b_{i+1,j-1})$ を計算し $b_{1, j}$ を出力します。

プログラム

簡単な分数計算用のルーチンを用いて有理数計算をしてみます。手抜きなのですぐにオーバーフローします。静的演算子オーバーロードを用いることで、メインルーチン側では、さほど意識せずに式を書けます。

ソースコード

module m_rat
  implicit none
  integer, parameter :: ki = selected_int_kind(15)

  type :: t_rat
    integer(ki) :: n, d
  end type t_rat

  interface operator(+)
    procedure :: rat_add
  end interface
  interface operator(-)
    procedure :: rat_sub
  end interface
  interface operator(*)
    procedure :: rat_mul
  end interface
  interface operator(/)
    procedure :: rat_div
  end interface

contains
  pure elemental type(t_rat) function reduce(a)
    type(t_rat), intent(in) :: a
    reduce = t_rat(0, 1)
    if (a%n == 0) return 
    reduce = t_rat(a%n / gcd(a%n, a%d), a%d / gcd(a%n, a%d))
    reduce%n = reduce%n * sign(1_ki, reduce%d)
    reduce%d = abs(reduce%d)
  contains
     pure integer(ki) function gcd(ia, ib)
       integer(ki), intent(in) :: ia, ib
       integer(ki) :: m(2)
       m = [ia, ib]
       do
         if  (mod(m(2), m(1)) == 0) exit
         m = [mod(m(2), m(1)), m(1)]
       end do   
       gcd = m(1)
     end function gcd 
  end function reduce  

  pure elemental type(t_rat) function rat_add(a, b)
    type(t_rat), intent(in) :: a, b
    rat_add = reduce( t_rat(a%n * b%d + a%d * b%n, a%d * b%d) )
  end function rat_add

  pure elemental type(t_rat) function rat_sub(a, b)
    type(t_rat), intent(in) :: a, b
    rat_sub = reduce( t_rat(a%n * b%d - a%d * b%n, a%d * b%d) )
  end function rat_sub

  pure elemental type(t_rat) function rat_mul(a, b)
    type(t_rat), intent(in) :: a, b
    rat_mul = reduce( t_rat(a%n * b%n, a%d * b%d) )
  end function rat_mul

  pure elemental type(t_rat) function rat_div(a, b)
    type(t_rat), intent(in) :: a, b
    rat_div = reduce( t_rat(a%n * b%d, a%d * b%n) )
  end function rat_div
end module m_rat

program Bernoulli
  use m_rat
  implicit none
  integer, parameter :: n = 32
  type(t_rat) :: b(n, n)
  integer(ki) :: i, j

  forall(i = 1:n) b(i, 1) = t_rat(1, i)
  do j = 2, n 
    forall(i = 1:n - j + 1) b(i, j) = t_rat(i, 1) * (b(i, j - 1) - b(i + 1, j - 1))
    print '(i2.2, a, 2i30)', j, ':', b(1, j) 
  end do 
end program Bernoulli

実行結果

分数の分子と分母を出力しています。
番号等が一般的なベルヌーイ数の定義と色々ずれていますが、気にしないでください。

02: 1 2
03: 1 6
04: 0 1
05: -1 30
06: 0 1
07: 1 42
08: 0 1
09: -1 30
10: 0 1
11: 5 66
12: 0 1
13: -691 2730
14: 0 1
15: 7 6
16: 0 1
17: -3617 510
18: 0 1
19: 43867 798
20: 0 1
21: -174611 330
22: 0 1
23: 854513 138
24: 0 1
25: -236364091 2730
26: 0 1
27: 8553103 6
28: 0 1
29: -23749461029 870
30: 0 1
31: 8615841276005 14322
32: 0 1
続行するには何かキーを押してください . . .

これ以上ではオーバーフローして出鱈目な値となります。

参考リンク：ベルヌーイ数
http://numbers.computation.free.fr/Constants/Miscellaneous/bernoulli.html