Fortranが書きたくない！！

わたしは大学で天気とか地震とか海流といった「地球物理」と呼ばれるジャンルを学んでいるのですが、このジャンルではコンピュータをガンガン回して計算しまくります。気象庁が毎日やってる天気予報もそうですし、「地球シミュレータ」でエルニーニョや地球温暖化の予測とかをやっているのもそうです。

で、そういう学科なので、コンピュータの実習があるのですが、なぜか今でもFortranを教えています。

たしかにFortranには、LAPACKのような優れた数値計算ライブラリがあったり、スパコンで書く時はコンパイラのサポートが充実してたりします。なので、そういう特定の用途には悪くない言語だと思うのですが、デスクトップPCで十分扱えるようなデータサイズで、そもそも計算速度はさして重要ではなくてデータ解析手法を学ぶ実習で使うのには辛すぎます。陽には書きたくありません。1でも書かないと単位くれないんだってさ！

なので、Fortranを書かずにFortranを書く方法を検討して、実際に実装しました。今回のソースコードも、すべてgithub上に上げています。

Fortranを書かずにFortranを書くには？

最初、Fortran上で小さなVMのインタプリタを実装して、そのVMの仮想機械語にコンパイルされる言語のコンパイラを実装しようと思いました。

が、VMを書くのが普通にしんどいのと、言語とコンパイラをこのためだけに実装するのもかなりつらいので、もっとラクな方法は無いか考えてみました。

冷静に考えると、別にVMのインタプリタを実装しなくとも、コンパイルされたFortranのコードは既にCPUという名前のx86機械語のインタプリタ機械（？）で実行されています。このインタプリタ（？）を利用できないでしょうか。利用できれば自力でVMのインタプリタを実装しなくていいので手間を省略できますし、コンパイラもgccとかclang/llvmとかすごいやつがいっぱい使えます！

具体的にはそのためにどうすれば良いのかというと、Fortranのinteger配列にC言語のソースをコンパイルしたx86の機械語の数字を沢山並べて入れて、そのinteger配列のポインタをC言語の関数ポインタに無理やりキャストして、そして実行すればいいのですっ！！

module a;interface;subroutine f (z) bind(c);use iso_c_binding;type(c_ptr), value :: z;
end subroutine;end interface;type G;procedure(f), pointer, nopass :: x;end type;type J;integer, pointer :: x;
end type;end module;program b;use iso_c_binding;use a;implicit none;type (J), pointer :: x(:);
real(8),pointer :: dbl(:);type(G), pointer :: p;integer, pointer :: d(:);integer :: i;integer :: n;allocate( x(1) );
allocate( dbl(10000000) );allocate( d(10000000) );dbl = 0;x(1)%x => d(645);d(1)=-443987883;d(2)=1213580125;d(3)=267576713;
d(4)=1223181585;d(5)=1223181709;d(6)=48087689;d(7)=450755289;d(8)=-128612024;d(9)=-398095544;d(10)=-532313784;
d(11)=1158680562;d(12)=1438866912;d(13)=1223002440;d(14)=1223705997;d(15)=-338050423;d(16)=-1991763235;
d(17)=-1958152107;d(18)=-1991708603;d(19)=267577413;d(20)=1575503120;d(21)=-1991748157;d(22)=286257893;
d(23)=-1924601787;d(24)=-1991710651;d(25)=-654123582;d(26)=1209720318;d(27)=1224234377;d(28)=1223181707;
d(29)=-220183159;d(30)=-532344817;d(31)=1213580125;d(32)=267576713;d(33)=1223181585;d(34)=1223181709;d(35)=48087689;
d(36)=450756569;d(37)=-128612024;d(38)=-398095544;d(39)=-532313784;d(40)=1158680562;d(41)=1438866912;
d(42)=-219838136;d(43)=-398126833;d(44)=-398095032;d(45)=-574453432;d(46)=-572401406;d(47)=1435060250;
d(48)=1166756088;d(49)=1166625000;d(50)=269480672;d(51)=-1017257915;d(52)=-443987883;d(53)=-394426040;
d(54)=-1192987255;d(55)=0;d(56)=-129660600;d(57)=16008647;d(58)=-352321536;d(59)=-196768982;d(60)=-1924622264;
d(61)=50452;d(62)=-1958215680;d(63)=21555269;d(64)=269480656;d(65)=269480448;d(66)=267581517;d(67)=267567448;
d(68)=-2080881391;d(69)=-1962806203;d(70)=1161557061;d(71)=1221491940;d(72)=1224230283;d(73)=-220707447;
d(74)=-666562545;d(75)=1213580125;d(76)=-1991711351;d(77)=-1991710595;d(78)=1435099253;d(79)=47324;d(80)=-1991770112;
d(81)=1170733125;d(82)=244;d(83)=-1958286592;d(84)=-1740049339;d(85)=-988508856;d(86)=0;d(87)=-398095544;
d(88)=-221249208;d(89)=-1962405873;d(90)=-1740049339;d(91)=-988508856;d(92)=0;d(93)=-532313272;d(94)=-221249208;
d(95)=-234876913;d(96)=-222209777;d(97)=-129167345;d(98)=-1051193358;d(99)=1158746098;d(100)=-196770824;
d(101)=-196769023;d(102)=2094810427;d(103)=1166756018;d(104)=1166625016;d(105)=269480656;d(106)=-1017262011;

（略）

d(761)=-800748728;d(762)=-389576376;d(763)=-1168;d(764)=1223710091;d(765)=1220572555;d(766)=-1178057333;d(767)=20;
d(768)=-389576376;d(769)=-683;d(770)=-1980742261;d(771)=535478722;d(772)=-120467455;d(773)=-223591031;
d(774)=-1404753393;d(775)=-2008708280;d(776)=1118194;d(777)=16008647;d(778)=-352321536;d(779)=-196768970;
d(780)=-2092394424;d(781)=-1924660800;d(782)=50452;d(783)=-1958215680;d(784)=21530693;d(785)=-196768830;
d(786)=-1924622264;d(787)=50444;d(788)=-1958215680;d(789)=21545029;d(790)=9128136;d(791)=-2096985784;
d(792)=-1962806203;d(793)=1161557061;d(794)=-1195213652;d(795)=0;d(796)=1213580233;d(797)=-1017256567;
call c_f_pointer(C_LOC( x(1) ), p);i=0;read (*,*), n;dbl(1) = real(n);do i=2,n+1;read (*,*), dbl(i);end do;call p%x( c_loc(dbl(1)) );
n = int(dbl(1));do i=2, n+1;write (*,*), i, dbl(i);end do;deallocate ( d );deallocate ( x );deallocate ( dbl );
contains;end program;

Fortran Side

Fortranでのキャストテクニック

今回はintegerの配列のポインタをCの関数へのポインタだと思わせたいため、（Cで言うと）int*をvoid (*ptr)(double*)のポインタにキャストする必要があります。しかし、Fortranは割りと型に厳しく、Cのように簡単にキャストしたり出来ません。なので、ちょっと撚（ひね）る必要があります。

!! 関数の型定義。
!! typedef void(*f)(void*);
interface
	subroutine f (z) bind(c)
		use iso_c_binding
		type(c_ptr), value :: z
	end subroutine
end interface
!! struct FPointer { f fptr; };
type FPointer
	procedure(f), pointer, nopass :: fptr
end type
!! struct IPointer { int* iptr; };
type IPointer
	integer, pointer :: iptr
end type

このように２つ構造体（IPointer,FPointer)を作っておき、この2つの構造体のポインタをC言語のポインタを経由することで無理やりキャストします。

type (IPointer), pointer :: ip(:);       !! IPointer* ip;
real(8),pointer :: dbl(:);               !! double* dbl;
type(FPointer), pointer :: fp            !! FPointer* fp;
type(C_PTR) :: cptr                      !! void * cptr;
integer, pointer :: d(:) !命令列用の配列  !! int* d;
allocate( ip(1) )                        !! ip = malloc(sizeof(IPointer));
allocate( dbl(10000000) )                !! dbl = malloc(sizeof(double)*10000000);
allocate( d(10000000) )                  !! d = malloc(sizeof(int)*10000000);
!! 命令列をセット
d(1)=-443987883;
d(2)=1213580125;
!! 略
d(790)=-443987883;
d(791)=50013;
!! 645番目の要素がC言語側エントリポイント関数の頭なので、そこへのアドレスをセット。
ip(1)%iptr => d(645);
!! IPointerへのポインタを、C言語のポインタ、いわばvoid*に落とす
cptr = C_LOC( ip(1) );
!! そこからさらにFPointerへのポインタにキャストする（専用の関数を使います）
call c_f_pointer(cptr, fp)
!! 無理やり変換した関数ポインタ経由でバイナリコードにジャンプ
call fp%fptr( c_loc(dbl(1)) )

IPointer 構造体へのポインタを一度Cのポインタ（いわばvoid*）にキャストし、それを再度void*からFPointerのポインタにキャストしています。構造体の中身の型が「intへのポインタ」と「関数へのポインタ」で違うので、構造体を通して、その中身についてキャストすることができました。

なんでこんな七面倒臭いことをしているかというと、intのポインタと関数のポインタを直接キャストすることができないからです。intのポインタは type(C_PTR)、関数のポインタはtype(C_FUNPTR)なので互換性がなく、さらにC言語みたいにintに無理やり落としてもう一度ポイ ンタにするような操作もできません。しかし、一回構造体でくるんでしまえばどちらも値のポインタとして扱い、キャストすることができます！

さらに、現在のLinuxでは、通常データ領域をプログラムとして実行することはできない（DEP）ので、コンパイルする時にその制限をはずす必要もあります。これには、 -z execstackというオプションスイッチが使えます。

これで、任意のx86機械語の列をFortranから実行することがきました。あとは、その機械語の列をどうやって作るかという話になりますね。原理的には自分で書いてもよいのですが、一応課題をやっているコードなので、結構複雑です。なので、C言語で書いたコードをコンパイルして、その結果を利用したいところです。

Binary Side

コードのコンパイルにはGCCを使うことにしましょう。

この時問題になるのは、Fortranの動的に配置した配列の上に機械語を置くので、機械語の置かれるアドレスが実行時になるまでわからないことです。なので、機械語がどこに置かれるかがわかっていることが前提の通常のコードはうまく動いてくれません。

そのようなケースは実は珍しくなく、共有ライブラリもアドレス上のどこにロードされるかは分からないので、ちゃんとそのためのコードをコンパイラが出力します2。そのようなコードのことを、PIC（Position Indepentent Code）と呼びます！gccでは、-fPICというコンパイルオプションをつけると、PICとしてコンパイルしてくれます。

じゃ、-fPICを付けてコンパイルしたバイナリから機械語をもってくれば、目的は達成できるのでしょうか…？

Position Independent Code

ところがどっこい、これではうまく行きません。なんでかというと、一個前に翻訳した記事で解説されている、PLTとGOTがあるからです。

こんな非常に簡単なC言語コードを考えてみましょう。

void calledFunction(){
}

void function(){
  calledFunction();
}

これをコンパイルするためのMakefileはこんな感じで、

.PHONY: all

all:
  gcc -c -o test.o test.c -fPIC
  ld -shared -fPIC -o test.so test.o

コンパイルした結果が、こんな感じでーす。

test.so:     file format elf64-x86-64


Disassembly of section .plt:

0000000000000280 <calledFunction@plt-0x10>:
280:    ff 35 82 0d 20 00        push   QWORD PTR [rip+0x200d82]        # 201008 <_GLOBAL_OFFSET_TABLE_+0x8>
286:    ff 25 84 0d 20 00        jmp    QWORD PTR [rip+0x200d84]        # 201010 <_GLOBAL_OFFSET_TABLE_+0x10>
28c:    0f 1f 40 00              nop    DWORD PTR [rax+0x0]

0000000000000290 <calledFunction@plt>:
290:    ff 25 82 0d 20 00        jmp    QWORD PTR [rip+0x200d82]        # 201018 <_GLOBAL_OFFSET_TABLE_+0x18>
296:    68 00 00 00 00           push   0x0
29b:    e9 e0 ff ff ff           jmp    280 <calledFunction@plt-0x10>

Disassembly of section .text:

00000000000002a0 <calledFunction>:
2a0:    55                       push   rbp
2a1:    48 89 e5                 mov    rbp,rsp
2a4:    5d                       pop    rbp
2a5:    c3                       ret

00000000000002a6 <function>:
2a6:    55                       push   rbp
2a7:    48 89 e5                 mov    rbp,rsp
2aa:    b8 00 00 00 00           mov    eax,0x0
2af:    e8 dc ff ff ff           call   290 <calledFunction@plt>
2b4:    5d                       pop    rbp
2b5:    c3                       ret

こんな感じで、ばっちりpltを経由して実行してしまいます。この機械語をFortran上に載せて強制的に呼び出したとしても、PLTがちゃんと解決されないとうまく動いてくれません。デフォルトのPLT解決処理は共有ライブラリとして実行されていることが前提になっているので、今回は利用することはできません。

integer配列のアドレスから関数のアドレスを解決することは原理的には可能なので、Fortran上でPLTのアドレスの所を埋めてちゃんと実行できるようにすること自体は可能（なはず）です。

しかし、かなり意味もなく複雑になってしまうので、できればPLTにジャンプするのではなくて、呼びたい関数に直接ジャンプするようなコードを出力して欲しいです。PLTを経由して関数を呼び出す時に使ってる0xe8命令は相対アドレスによるcall命令なので、可能なはずです。

リンカースクリプトを自分で書けばPLTとGOTを使わない…！？

リンク時にPLTを経由するようにリンクしているので、じゃあすごく単純なリンカースクリプト自分で書けばPLTみたいな便利機能使わなくなるんじゃない？

ということで、実際に書いてみました。こうなりました。

OUTPUT_FORMAT("elf64-x86-64", "elf64-x86-64","elf64-x86-64")
OUTPUT_ARCH(i386:x86-64)

ENTRY(ep);
MEMORY
{
ROM(rxai) : ORIGIN = 0, LENGTH = 64k
}
SECTIONS
{
.text : {} > ROM
.rodata : {} > ROM
.data : {} > ROM
. = ALIGN(4);
__bss_start = .
; .bss : {} > ROM
__bss_end = . ;
}

Makefileは、

linkerscript:
    gcc -c -o test.o test.c -fPIC 
    ld test.so test.o --script ./linker.script

こんな感じにすると指定できます。結果は、

test-with-linkerscript:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <calledFunction>:
   0:    55                       push   rbp
   1:    48 89 e5                 mov    rbp,rsp
   4:    5d                       pop    rbp
   5:    c3                       ret    

0000000000000006 <function>:
   6:    55                       push   rbp
   7:    48 89 e5                 mov    rbp,rsp
   a:    b8 00 00 00 00           mov    eax,0x0
   f:    e8 ec ff ff ff           call   0 <calledFunction>
  14:    5d                       pop    rbp
  15:    c3                       ret

となって、見事、PLTを使わずに直接相対ジャンプするようになりました！

実は-pieだけで大丈夫

実をいうと、リンカースクリプトなんて書かなくても、-pieというオプションを付けて、通常の実行ファイルとしてコンパイルすればplt/gotを使わないコードを生成してくれます。

pie:
    gcc -c -o test.o test.c -pie
    ld -pie test.so test.o

-pieというオプションは、gcc –helpすると、

  -pie                     Create a position independent executable

ということで、実行時にどのアドレスに配置されても実行できる実行ファイルを生成してくれます。

実行ファイルは共有オブジェクトとは違って外部の他のプログラムからその中の関数が呼ばれるようなことはないので、内部のそれぞれの関数のアドレスを動的に解決する必要はありません。なので、PLTを経由せずに直接callするようなコードにしてくれる…のでしょう。

…でも冷静に考えると、共有ライブラリの中でならPLT経由せずに直接ジャンプでよくない？

全てが終わったあとに気づいたので時既に遅しという感じでしたが、一応書いておきます。

プログラムセクションを抽出する

残るところは、出来上がったファイルから機械語を抽出して、数字の列にするだけです。そのために、機械語が書いてある部分はファイルのなかでどこなのかを調べましょう。

セクションの一覧を出力するには、objdump -hが使えます。

% objdump -h test.so

test.so:     file format elf64-x86-64

Sections:
Idx Name          Size      VMA               LMA               File off  Algn
0 .interp       0000000f  00000000000001c8  00000000000001c8  000001c8  2**0
CONTENTS, ALLOC, LOAD, READONLY, DATA
1 .hash         00000028  00000000000001d8  00000000000001d8  000001d8  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
2 .dynsym       00000078  0000000000000200  0000000000000200  00000200  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
3 .dynstr       00000019  0000000000000278  0000000000000278  00000278  2**0
CONTENTS, ALLOC, LOAD, READONLY, DATA
4 .text         00000016  0000000000000294  0000000000000294  00000294  2**2
CONTENTS, ALLOC, LOAD, READONLY, CODE
5 .eh_frame     00000058  00000000000002b0  00000000000002b0  000002b0  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
6 .dynamic      000000c0  0000000000200f40  0000000000200f40  00000f40  2**3
CONTENTS, ALLOC, LOAD, DATA
7 .comment      0000002c  0000000000000000  0000000000000000  00001000  2**0
CONTENTS, READONLY

機械語本体が置かれているのは、.textというセクションです。この中のデータをFortranに書きましょう。File offの項目が、ファイルの先頭からの位置で、Sizeという項目がセクションのサイズで、今回の場合、要はプログラムコード全体のサイズです。

さらに、関数呼び出しするには、呼び出したい関数の、プログラムコード内での相対位置も必要です。これには、nmコマンドが使えます。

% nm test.so
0000000000201000 D __bss_start
0000000000000294 T calledFunction
0000000000200f40 d _DYNAMIC
0000000000201000 D _edata
0000000000201000 D _end
000000000000029a T function
0000000000000000 d _GLOBAL_OFFSET_TABLE_
                 U _start

セクションの所のVMAもしくはLMAというのがメモリ上での位置なので、機械語列の中での相対位置を調べるにはここでのアドレスから.textセクションのVMA/LMAの値を引けば求めることができます。

シェルスクリプトを使うともちろんこの作業は自動化できます。呼び出す関数名は”ep”としました。EntryPointの略ですね！

_start=$(objdump -h test.so | grep \\s.text\\s | awk '{ print $6 }')
_size=$(objdump -h test.so | grep \\s.text\\s | awk '{ print $3 }')
_ep=$(nm test.so | grep \\sep$ | awk '{ print $1 }')

Cのソースをコンパイルする所から、命令列を抽出してFortranのソースに変換する作業まで、これらの作業はすべて自動化しているので、よかったらgithubのソースを見てね！

今回書き込んだinteger型は4バイトなので、エントリポイントのオフセットが4の倍数でないとき（よくある）には、その整列条件を満たすために数バイトズラしたりしないといけません。そのコードとかも入ってます。

Programming Side

これで終わり…ではありません！もう一個、プログラムを書くという問題があります！！今回の環境では、

• 外部にあるライブラリ関数などは呼び出すことができない

という、非常に大きな制限があります。実行時に決定されるライブラリ関数の実際のアドレスを知らないし、知るためにはライブラリ関数を呼び出さないといけないので、「服を買う服がない」状態なんです。もちろん、Fotranでなら調べられるので、調べた結果をC言語の関数に投げればいいんですけど、Fortran書きたくないし（ぉ）。

科学計算すればいいだけなので、sin/cosとsqrtくらいが使えればなんとかなります。sin/cos/sqrt/absは、FPU命令を使って計算できるので、ライブラリは必要ありません。

そのためには、インラインアセンブリで陽にFPU命令を実行する必要があります。

double fsin(double v)
{
	intptr_t d0;
	__asm__ (
		".intel_syntax noprefix;"
		"fld qword ptr [%0];"
		"fsin;"
		"fstp qword ptr [%0];"
		".att_syntax;"
		: "=&r"(d0)
		: "0"((intptr_t) &v)
		: "eax", "ecx"
	);
	return v;
}

こんな感じのを、FPUの命令ごとに実装していきます。

関数だけでなく、piもFPUを使って取得しています。なんでかというと、今回は機械語だけを埋め込んで、それとは別にある定数の領域をFortranのソースに含めていないからです！

Cのソースに円周率の値を直接書いていてもすこし複雑な処理になると定数領域に飛ばされてしまって、.textセクションに含まれなくなってしまい、結果としてプログラムがおかしくなってしまいました。定数領域もFortranに埋め込むだけでその問題は回避できるのですが、その改修がちょっとめんどくさかったのでやめて、FPU命令を毎回呼ぶことで回避しました。

入力と出力ですが、実際のデータの読み込みと結果の出力のみ、Fortranが担当することにしました。C言語は関数のシグネチャをvoid (*fun)(double *)とし、doubleの配列を処理してその結果をもとの配列の領域に書き込むことだけに専念していただきます。

他の言語でやってみても面白いかもね！

と思った。Fortranはなんだかんだ任意の機械語を実行する難易度は低い言語だったと思います。

ここまでやるならFortranで書いたほうが早くね？

なんでFortranせなあかんねんと。その時間を使ってね、バイナリおじさんをする。そういう考え方もできると思うんですよ。だってスパコン向けでもないのにFortran書いてもこころぴょんぴょんしないでしょ。

— ψ（プサイ） (@tikal) June 14, 2014

ちなみに

書かなくて良くなった(╹◡╹)

高校宜しく職員室ならぬ教員室で面談したらFotran以外で書いてもOKって言われたよーやったー(╹◡╹) #プログラミング言語選択の自由は基本的人権

— ψ（プサイ） (@tikal) June 3, 2014

ひさびさに色々バイナリが触れてたのしかったです。まぁ結果としては良かったのではないでしょうか！

工学部「Fortran？書いたこと無いですね…」 文系「Fortranって何？」 理学部「（他にもっとうまい人いるしパソコンは計算機としか思ってないし）Fortranあんまり出来ないですね…」 わたし「（機械語をソースに直接書いて実行できるから）Fortranできます！！！！」

— ψ（プサイ） (@tikal) June 4, 2014

1. っていうか、いままでプログラム書いたことがない人にいきなりこれ教えるのってどうなんでしょう…。 [↩]
2. Windowsでは、64bit版からそうなりました。32bitでは固定だそうです。 [↩]