古代言語FORTRANの解読について
計算科学に関わる研究室では, FORTRANで書かれた秘伝のソースが醸成され, 継承されています.
教授「FORTRANわかる？」
学生「Fortranわかります！」
なんていう会話が聞こえてきますね. FORTRANは現代の最新エディタを以てしても扱うことが困難な古代技術ですから, この学生は苦しむことになります.
FORTRANとは何か
ここではFORTRAN 77以前のものをFORTRAN, Fortran 90以降のものをFortranと呼びます. FORTRANの時代遅れとされた機能はいくつかありますが, 中でも固定形式が特徴的です.
固定形式のプログラムの例としては, A.ザボ, N.S.オストランド『新しい量子化学（上） 電子構造の理論入門』の付録などをご覧ください. これは Computational Chemistry List (CCL) で配布されているプログラムです. かなり読みやすい方ですが, 固定形式に慣れていないと読めません.
図書館で30年前の数値計算の教科書を見ると, PascalやBasicが使われていた痕跡がありますが, もう誰も使っていません. Pythonの2系も動かないものばかりです. 一方, 後方互換性が守られているため, FORTRANは未だに現役で, 最新のコンパイラでも問題なくビルドできてしまいます. 実際, ザボの付録も警告は出ますがgfortranでビルド, 実行できます. なお, WSL上でのgfortranの使い方はこちらを参考にしてください.
Invoke-WebRequest -Uri http://www.ccl.net/cca/software/SOURCES/FORTRAN/szabo/szabo.f -OutFile "szabo.f"
gfortran szabo.f -o szabo.exe
./szabo.exe
wget http://www.ccl.net/cca/software/SOURCES/FORTRAN/szabo/szabo.f
gfortran szabo.f -o szabo
./szabo
書き換える必要性が生じない限りは無理にいじらないことが最善ですが, 不運にもFORTRANを書き換える必要性に直面してしまった哀れな子羊はFORTRANを学ぶ必要があります. しかし, 今どき講義でFORTRANを習うこともないでしょうから, 秘伝のソースを継承した際に苦しみながら調べることになります.
何がつらいのか？
プログラムが全部大文字だとか, インデントが無いという時点で読む気が失せますが, ツールで解決できますし, 実は大きな問題ではありません. グローバル変数もかなりヤバいですが, これはFORTRANに限った話ではありません. 最大の原因は文番号とGOTO文であり, この文法が理解できないことがつらいわけです. 早速, 文番号とGOTO文の扱いを中心に, FORTRANをFortranに書き換えていきましょう.
コーディングルール
これは飽くまでも個人的に提案するコーディングルールですが, このルールに適合するコードに変換していきます.
- 基本的に全て小文字で書きましょう. 大文字は読みにくいです.
- 文番号の使用は禁止です（do文, if文, write文のfmtなど）.
- continue文も禁止です（
doとifには必ずend doとend ifを対応させます）. - 全て
implicit noneにします. -
implicit integer (I-N)は禁止です. -
implicit double precision (a-h,o-z)は禁止です. - グローバル変数は禁止です.
- common文やmoduleによる変数の渡しも禁止です（グローバル変数と同様, 明示的に変数を引数を書くことを薦めます. 長くなる場合はオプショナル変数などを活用しましょう）.
- グローバル変数も原則禁止です.
- 空行のコメントアウトは不要です.
- 適切なインデントをしましょう.
Hello World!
これは特に解説は不要でしょう. 一応, こちらを参考に逆アセンブルした結果も完全に一致します.
SUBROUTINE HELLO()
PRINT *, "HELLO"
END SUBROUTINE
subroutine hello()
print *, "HELLO"
end subroutine
HELLO
FORMAT文
編集記述子についてはこちらを参考にしてください.
PRINT 10, 3.1415926535
10 FORMAT(F8.4)
PRINT 20, "HELLO", 1234, 3.1415926535
20 FORMAT(3X,A10,I8,F8.4)
print '(f8.4)', 3.1415926535
print '(3x,a10,i8,f8.4)', "HELLO", 1234, 3.1415926535
3.1416
HELLO 1234 3.1416
出力は同じですが, 自由形式に書き換えると逆アセンブルした結果が一行だけ変わりました. 何による違いなのかが全く分かりません. 有識者のコメントをお待ちしております.
- 96: c7 85 30 fe ff ff 04 movl $0x4,-0x1d0(%rbp)
+ 96: c7 85 30 fe ff ff 03 movl $0x3,-0x1d0(%rbp)
DO文
doには必ずend doを対応させましょう.
DO 10 I=1,2
10 PRINT *, I
do i=1,2
print *, i
end do
1
2
上記を見ると「DO 10 I=1,2は10 PRINT *, Iを2回実行する」という意味なのか？と勘違いしそうになりますが「DO 10 I=1,2は10 PRINT *, Iまでを2回実行する」という意味です. こちらの例を見てください.
DO 20 I=1,2
PRINT *, I
20 PRINT *, I
do i=1,2
print *, i
print *, i
end do
1
1
2
2
このように「DO 20 I=1,2は文番号20でまで行ったら戻る」という意味なんですね. これがなかなか理解できずに苦しみました. これを応用すると, CONTINUE文という何もしない命令をend doの代わりに使うことができます.
DO 30 I=1,2
PRINT *, I
30 CONTINUE
do i=1,2
print *, i
end do
1
2
さて, do文の意味は理解できたと思いますので, 2つのdo文が関わる例を見ていきましょう. 下記の例は単純にdo文が2つ並んでいるだけです.
DO 40 I=1,2
40 PRINT *, I
DO 50 J=3,5
50 PRINT *, J
do i=1,2
print *, i
end do
do j=3,5
print *, j
end do
1
2
3
4
5
上記の例と似ていますが, 下記の例は二重ループです. 60 CONTINUEが2つ分のend doに相当します.
DO 60 I=1,2
PRINT *, I
DO 60 J=3,5
PRINT *, I, J
60 CONTINUE
do i=1,2
print *, i
do j=3,5
print *, i, j
end do
end do
1
1 3
1 4
1 5
2
2 3
2 4
2 5
次は素直な二重ループですね. end doに相当する部分が無いのでわかりにくいです.
DO 70 I=1,2
DO 70 J=3,5
70 PRINT *, I, J
do i = 1,2
do j = 3,5
print *, i, j
end do
end do
1 3
1 4
1 5
2 3
2 4
2 5
最後は 80 CONTINUEが2つ分のend doに相当する素直な二重ループの例です.
DO 80 I = 1,2
DO 80 J = 3,5
PRINT *, I, J
80 CONTINUE
do i = 1,2
do j = 3,5
print *, i, j
end do
end do
1 3
1 4
1 5
2 3
2 4
2 5
IF文
恐らく, FORTRANからFortranへの変換の自動化を困難にしている要因の一つです. 昔のプログラムを見ていると「trueなら実行しない」という形でif文が使われているのをよく目にしますが, これはgoto文と併用するせいでしょうか. 下記は, もし 10 CONTINUEまで移動するといものです.
IF (I.EQ.0) GO TO 10
PRINT *, 'I /= 0'
10 CONTINUE
if (i/=0) then
print *, 'I /= 0'
end if
これをifとend ifで挟む書き方に直そうとするとi==0ではなくi/=0を判定しなければいけません. フローチャートを描いてから正しい判定文を逆算しないと間違える可能性があります. このあたりが自動化しにくいポイントでしょうか？ また, 結果は同じでもコードとしては変更が大きいせいか, 逆アセンブルの結果もかなり変わってしまいました.
IF (I.EQ.0) GO TO 20
PRINT *, 'I /= 0'
IF (I.EQ.1) GO TO 20
PRINT *, 'I /= 1'
20 CONTINUE
if (i/=0) then
print *, 'I /= 0'
if (i/=1) then
print *, 'I /= 1'
end if
end if
最後に, いわゆる算術IF文というラスボスです. これはかなり頭を使いますので, 最も簡単な例を示します. 分岐先がif文より後に出てくる素直に例を紹介します.
IF (I-1) 30, 40, 50
30 PRINT *, 'I < 1'
40 PRINT *, 'I < 1 or I = 1'
50 PRINT *, 'I < 1 or I = 1 or I > 1'
if (i<1) print *, 'I < 1'
if (i<1 .or. i==1) print *, 'I < 1 or I = 1'
if (i<1 .or. i==1 .or. i>1) print *, 'I < 1 or I = 1 or I > 1'
if (i<1) then
print *, 'I < 1'
end if
if (i<1 .or. i==1) then
print *, 'I < 1 or I = 1'
end if
if (i<1 .or. i==1 .or. i>1) then
print *, 'I < 1 or I = 1 or I > 1'
end if
あまりお行儀は良くないですが, thenとend ifを省略した書き方も示しました. 必ず同じように書き換えられるわけではなく, 分岐先の行番号によって違いますので気を付けてください. 分岐先がif文より前に出てくる例では, while文への書き換えが有効な場合があります.
do while文
FORTRANだとif文とgoto文でループが実装できてしまいますが, 時代遅れな書き方です. do whileを使いましょう.
I = 0
10 CONTINUE
PRINT *, I
I = I + 1
IF (I.LT.5) GO TO 10
i = 0
do while (i<5)
print *, i
i = i + 1
end do
0
1
2
3
4
ザボの付録の568行目に算術IF文による分岐を見つけてしまったので, 例示しておきます. 上記のコードに対して, if文が算術IF分に置き換わっただけです. IF (I-5) 20, 20, 30は「I<5のときとI==5のときは20番へ, I>5のときは30番へ」という意味なので, 単に<が<=に置き換わるだけで済みます.
I = 0
20 CONTINUE
PRINT *, I
I = I + 1
IF (I-5) 20, 20, 30
30 CONTINUE
i = 0
do while (i<=5)
print *, i
i = i + 1
end do
DIMENSION文とDATA文
2022/11/14に追記しました. 変数や配列に値を代入していきます. 配列構成子についてはこちらとこちらを参考にしてください.
IMPLICIT INTEGER (A-Z)
DIMENSION B(3)
DIMENSION C(3)
DIMENSION D(3)
DIMENSION E(5)
DIMENSION F(3)
DIMENSION G(3)
DIMENSION H(2,3)
DATA A/1/
DATA B/1,2,3/
DATA C/3,6,9/
DATA D/3*4/
DATA E/3*4,5,6/
DATA F,G/1,2,3,4,5,6/
DATA H/1,2,3,4,5,6/
implicit none
integer :: A = 1
integer :: B(3) = [1,2,3]
integer :: C(3) = [(3*i, i=1,3)]
integer :: D(3) = [(4, i=1,3)]
integer :: E(5) = [(4, i=1,3), 5, 6]
integer :: F(3) = [1,2,3], G(3) = [4,5,6]
integer :: H(2,3) = reshape([1,2,3,4,5,6], [2,3])
PRINT *, A -> 1
PRINT *, B -> 1 2 3
PRINT *, C -> 3 6 9
PRINT *, D -> 4 4 4
PRINT *, E -> 4 4 4 5 6
PRINT *, F -> 1 2 3
PRINT *, G -> 4 5 6
PRINT *, H(1,:) -> 1 3 5
PRINT *, H(2,:) -> 2 4 6
DIMENSION文は取り除きました. DATA文も時代遅れです. まず, Aについては型をIMPLICIT INTEGER (A-Z)で決めており, いきなりDATA A/1/と書いていますが, 明示的にinteger :: A = 1と書き換えましょう. B以降も同様に整数型ですが, 配列なのでDIMENSION文でサイズを指定しています. C(3)についてはDATA C/3,6,9/のように規則性があるのでinteger :: C(3) = [(3*i, i=1,3)]のように配列構成子の内部でdo文を使って置き換えることができます. D(3)についてはDATA D/3*4/とDATA D/4,4,4/が同じ意味になっています. E(5)についても同様の文法でDATA E/3*4,5,6/とDATA E/4,4,4,5,6/が同じ意味になっています. F(3),G(3)についてはDATA F,G/1,2,3,4,5,6/がDATA F/1,2,3/とDATA G/4,5,6/を合わせたものに対応しており, 2つの配列に値を代入することができます. さらにE(5)のような略記と組み合わせてDATA F,G/1,4*2,3/のようにされると, 初見では意味が分からないので苦戦しました（DATA F,G/1,4*2,3/はDATA F,G/1,2,2,2,2,3/と同じ意味で, DATA F/1,2,2/とDATA G/2,2,3/に分解できます）. 最後に2次元配列H(2,3)については, 次のような行列
を作っていることになります. DATA文はメモリの並び順の通りにデータを格納するようですので, DATA H/1,2,3,4,5,6/は
H(1,1) = 1
H(2,1) = 2
H(1,2) = 3
H(2,2) = 4
H(1,3) = 5
H(2,3) = 6
do j = 1,3
do i = 1,2
H(i,j) = i + 2*(j-1)
end do
end do
H(:,1) = [1,2]
H(:,2) = [3,4]
H(:,3) = [5,6]
などと同じことです. reshape(値, 配列のサイズ)を使って1行で初期値を指定することもでき, それが先ほどの例です.
おわりに
デバッグ方法について
- 出力の比較
- こちらを参考に逆アセンブル
- GDBによるステップ実行で1行ずつ確認する
などを考えましたが, 自動かつ確実な方法は無いかと探しています（それがあったら苦労してないんでしょうが...）. FORTRANからFortranの書き換えでは逆アセンブルは案外, 役に立たないかもしれないことがわかったので, GDBによるステップ実行を記録し, 書き換えの前後で比較するなどの方法が良いのでしょうか.
要は書き換えの前後でフローチャートが変わっていないことを保証する便利な方法があれば良いわけです. LLVMのFortranコンパイラも出たらしく, 制御フローグラフも出せるようなので, もしかすると使えるかもしれませんね. あるいは「全ての行にprint文を埋め込んでビルド・実行」するスクリプトなどを作れば, 変更前後でGoTo文を正しく排除できたか判定する一助となりそうです. (ホワイトボックステストの自動化もしたいところです)
自動リファクタリングは我々の悲願です. こちらやこちらなど, いくつかの論文を見つけたので別の記事で試してみる予定です.
また, 近年のHPC御三家（C, C++, Fortran）に加えてJuliaという選択肢も台頭してきました. ccall()による呼び出しでJuliaからCやFortranを呼び出すことはできますが, それはお上品なコードでの話であって, グローバル変数が入り乱れ, 機能がサブルーチンや関数として分離されていないようなプログラムでは使えません. 腹を括って, まずはFORTRANをFortranに, そしてFortranをJuliaに書き換えましょう. どうしてもコードを変更せずに呼び出したい場合はこちらで解説した方法でJuliaから呼び出すことができます. Juliaでネイティブ実装した方が, パッケージマネージャーやBigFloat型などの恩恵を受けることができますので, Fortranから始めるJuliaやJuliaへの変換を支援するスクリプトを使ったJuliaへの移植についても別の記事でまとめたいです.
JuliaでFortranのコンパイル・逆コンパイルなどを実行できるようになれば, Juliaによるリファクタリングツールも充実しそうな気がします.
他に参考になりそうな情報
この記事で紹介したプログラムは出力が書き換えの前後で変わらないことを確認したものです. プログラムと一緒にビルド・出力の比較・逆アセンブルを行うスクリプトは全てこちらのリポジトリに置いてあります.
Discussion