Akiyama-Tanigawa 法によるベルヌーイ数の計算 ref: http://qiita.com/cure_honey/items/9130bc26a20da7068c10

file0.f

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