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

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

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

計算は中身を読みもせず、図を見てそのままプログラムにします。
![Beroulli.png](https://qiita-image-store.s3.amazonaws.com/0/44522/98b72fd6-9d2a-2fac-7553-bdccb555a95c.png)

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

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


##ソースコード
```fortran
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