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]$$
サブちゃん
計算結果
以下ではまずベルヌーイ数を有理数表現での分子・分母を３０まで表示し、次に上の式に従って求めたオイラー定数と、誤差の目安のために最後に足し引きした項を表示しています。
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 で止めずにもう少し先まで計算すれば、もっと良い精度が得られると思います。
ソース・プログラム
オイラー定数は分数のまま計算すると桁あふれするので、各項を四倍精度実数に直して計算しています。
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