submodule でインターフェース分離

# submodule
submodule とは Fortran 2008 で導入された Module をインターフェース部と実装部に分離する仕組みです。

Fortran 90 で導入された module は、インターフェース部分をコンパイラ側が mod ファイルとして自動で作ってくれます。これは C 言語などのようにヘッダーファイルを人間が用意する必要がないので、同じような宣言部分をわざわざもう一度書く必要もないし、また修正するうちに矛盾が生じることも無いのではなはだ便利なのですが、いくつか問題もあります。

一つは実装部のソースコードを隠したいとき煩雑になる事で、いま一つはインターフェースは変えないまま実装部を修正してもインターフェース部も再構成されてしまうため、そのモジュールを使っているモジュールが雪崩的に不必要に再コンパイルされてしまうということです。

そこで Fortran 2008 では submodule という仕組みが新たに用意され、必要とあらばインターフェース部と実装部が分離できるようになりました。また、この時引数型などの宣言部を二度繰り返して書かなくて済むような仕組みも用意されました。（繰り返して書いていいですけど。）

以下にその実例を示します。以前簡単な分数計算ルーチンを作ったので、それを改造してみました。
https://qiita.com/cure_honey/items/9130bc26a20da7068c10 

今回は求めたベルヌーイ数を使ってオイラー・マクローリンの定理によってオイラー定数を計算してみます。逆数和は∞までのところ 10 で止めています。
$$\gamma=\lim_{n\rightarrow\infty}\left(\sum_{i=1}^n{1\over i}-\ln n\right)\sim\sum_{i=1}^{10}{1\over i}-\ln10 + \left[ {B_1 \over10} +\sum_{k=1}^{\infty}{B_{2k}
 \over 2k\,10^{2k}}\right]$$


![61jhAr2NjGL.jpg](https://qiita-image-store.s3.amazonaws.com/0/44522/0885af0d-1008-b0ac-65a6-1ca98b27c73d.jpeg)
サブちゃん

## 計算結果
以下ではまずベルヌーイ数を有理数表現での分子・分母を３０まで表示し、次に上の式に従って求めたオイラー定数と、誤差の目安のために最後に足し引きした項を表示しています。

```
                     0                     1                     1
                     1                    -1                     2
                     2                     1                     6
                     3                     0                     1
                     4                    -1                    30
                     5                     0                     1
                     6                     1                    42
                     7                     0                     1
                     8                    -1                    30
                     9                     0                     1
                    10                     5                    66
                    11                     0                     1
                    12                  -691                  2730
                    13                     0                     1
                    14                     7                     6
                    15                     0                     1
                    16                 -3617                   510
                    17                     0                     1
                    18                 43867                   798
                    19                     0                     1
                    20               -174611                   330
                    21                     0                     1
                    22                854513                   138
                    23                     0                     1
                    24            -236364091                  2730
                    25                     0                     1
                    26               8553103                     6
                    27                     0                     1
                    28          -23749461029                   870
                    29                     0                     1
                    30         8615841276005                 14322
 0.577215664901532860606515837096008
 0.000000000000000000000020052695797
続行するには何かキーを押してください . . .
```

文献値では $\gamma=0.57721 56649 01532 86060 65120 90082 $となっています。誤差の見積もり項のあたりまで正しい値が得られています。

逆数和を 10 で止めずにもう少し先まで計算すれば、もっと良い精度が得られると思います。

## ソース・プログラム
オイラー定数は分数のまま計算すると桁あふれするので、各項を四倍精度実数に直して計算しています。

```fortran
module m_mod
  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

  interface 
    module pure elemental type(t_rat) function reduce(a)
      type(t_rat), intent(in) :: a
    end function reduce  

    module pure elemental type(t_rat) function rat_add(a, b)
      type(t_rat), intent(in) :: a, b
    end function rat_add

    module pure elemental type(t_rat) function rat_sub(a, b)
      type(t_rat), intent(in) :: a, b
    end function rat_sub

    module pure elemental type(t_rat) function rat_mul(a, b)
      type(t_rat), intent(in) :: a, b
    end function rat_mul
    
    module pure elemental type(t_rat) function rat_div(a, b)
      type(t_rat), intent(in) :: a, b
    end function rat_div
  end interface 
end module m_mod

submodule (m_mod) m_mod_implement
contains
  module procedure reduce 
    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  (m(2) == 0) exit
         m = [m(2), mod(m(1), m(2))]
       end do   
       gcd = m(1)
     end function gcd 
  end procedure reduce  

  module procedure rat_add
    rat_add = reduce( t_rat(a%n * b%d + a%d * b%n, a%d * b%d) )
  end procedure rat_add

  module procedure rat_sub
    rat_sub = reduce( t_rat(a%n * b%d - a%d * b%n, a%d * b%d) )
  end procedure rat_sub

  module procedure rat_mul
    rat_mul = reduce( t_rat(a%n * b%n, a%d * b%d) )
  end procedure rat_mul

  module procedure rat_div
    rat_div = reduce( t_rat(a%n * b%d, a%d * b%n) )
  end procedure rat_div
end submodule m_mod_implement

program hello
  use m_mod
  implicit none
  integer, parameter :: n = 32, kd = kind(1.0q0)
  type(t_rat) :: b(n, n), bn(0:n - 1), g, f
  integer(ki) :: i, j
  real(kd) :: d

  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))
  end do 
  !
  ! Bernoulli number Bn 
  bn(0) = t_rat( 1, 1)
  bn(1) = t_rat(-1, 1) * b(1, 2)
  do j = 2, n - 1
    bn(j) = b(1, j + 1)
  end do
  ! Sum(1/i, i=1:10)
  f = t_rat(0, 1)
  do j = 1, 10 
    f = f + t_rat(1, j)
  end do
  !
  ! -sum(bn(i)/i 10^i, i = 1 & 2:30:2)
  d = bn(1)%n / real(bn(1)%d, kd) / 10.0_kd  
  do j = 2, 30, 2  
    d = d + bn(j)%n / real(bn(j)%d, kd) / real(j, kd) / 10.0_kd**j 
  end do
  !
  do j = 0, 30
    print *, j, bn(j)%n, bn(j)%d 
  end do  
  print '(f36.33)', f%n / real(f%d, kd) - log(10.0_kd) + d
  print '(f36.33)', bn(30)%n / real(bn(30)%d, kd) / 30 / 10.0_kd**30
end program Hello
```