🦕

古代言語FORTRANの解読について

2022/08/24に公開

計算科学に関わる研究室では, 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の使い方はこちらを参考にしてください.

Windowsの場合
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
Ubuntuの場合
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文も禁止です(doifには必ずend doend 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

出力は同じですが, 自由形式に書き換えると逆アセンブルした結果が一行だけ変わりました. 何による違いなのかが全く分かりません. 有識者のコメントをお待ちしております.

fotmat.od
-  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,210 PRINT *, Iを2回実行する」という意味なのか?と勘違いしそうになりますが「DO 10 I=1,210 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文と併用するせいでしょうか. 下記は, もしi=1なら間を飛ばして 10 CONTINUEまで移動するといものです.

固定形式
      IF (I.EQ.0) GO TO 10
      PRINT *, 'I /= 0'
   10 CONTINUE
自由形式
  if (i/=0) then
    print *, 'I /= 0'
  end if

これをifend 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

あまりお行儀は良くないですが, thenend 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)については, 次のような行列

H = \begin{pmatrix} 1 & 3 & 5\\ 2 & 4 & 6\\ \end{pmatrix}

を作っていることになります. 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行で初期値を指定することもでき, それが先ほどの例です.

おわりに

デバッグ方法について

  1. 出力の比較
  2. こちらを参考に逆アセンブル
  3. 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から始めるJuliaJuliaへの変換を支援するスクリプトを使ったJuliaへの移植についても別の記事でまとめたいです.

JuliaでFortranのコンパイル・逆コンパイルなどを実行できるようになれば, Juliaによるリファクタリングツールも充実しそうな気がします.

他に参考になりそうな情報

https://www.nag-j.co.jp/fortran/FI_17.html
https://qiita.com/cure_honey/items/e06b89e238c3df3df693
https://qiita.com/implicit_none/items/55c47407aa376277a531
https://qiita.com/implicit_none/items/bc113fe438cfb07e0a44#gotoの飛び先多重ループからの脱出-1

この記事で紹介したプログラムは出力が書き換えの前後で変わらないことを確認したものです. プログラムと一緒にビルド・出力の比較・逆アセンブルを行うスクリプトは全てこちらのリポジトリに置いてあります.
https://github.com/ohno/F77f90

Discussion