前回の続きです。
以下のコードを-O3でコンパイルすると、
pure function nth(v, n) result(out) implicit none integer, intent(in) :: v(:) integer, intent(in) :: n integer :: out out = v(n) return end function nth
こんなアセンブリになって驚いた。
nth_: mov rdx, QWORD PTR [rdi+40] mov eax, 1 mov rcx, QWORD PTR [rdi] test rdx, rdx cmove rdx, rax movsx rax, DWORD PTR [rsi] imul rax, rdx sub rax, rdx mov eax, DWORD PTR [rcx+rax*4] ret
参考： godbolt.org
これを解読していきたいと思う。
ここで、いくつか思い出さないといけないことがある。Fortranは多次元配列を組み込みで持っており、組み込み関数としてreshapeなんかも持っているということ。それと、Fortranでは特に指定されない限り配列の添字が1から始まることだ。これらを頭に入れて、このアセンブリを読んでいこう。
長いが、とりあえず何をしているかそのまま言葉にしてみる。
nth_: ; 第一引数から40バイト離れたところから64bitロードしてrdxに代入 mov rdx, QWORD PTR [rdi+40] ; raxの下位32bitに1を代入 mov eax, 1 ; rcxに第一引数の先から64bitロード mov rcx, QWORD PTR [rdi] ; rdxがゼロならZF=1に test rdx, rdx ; ZF=1ならrdxにrax(今1が入っている)を代入 cmove rdx, rax ; raxに第二引数の先から32bitロードし、符号拡張して代入 movsx rax, DWORD PTR [rsi] ; rax *= rdx imul rax, rdx ; rax -= rdx sub rax, rdx ; rcxからrax*4バイト先の32bitをロードしてeaxに代入 mov eax, DWORD PTR [rcx+rax*4] ; return eax; ret
最後、戻り値レジスタeaxにアドレスの先からDWORD(32bit)ロードしてreturnしていることから、このときrcxが配列の先頭を指していて、raxが添字になっていることがわかる（4はintegerのサイズ）。raxには色々な計算がされているが、rcxは単にrdiの先からポインタがロードされているだけだ。ここから、第一引数は先頭ポインタを指すポインタであることがわかる。
同時に、1行めからrdi+40の先を読み込んでいることに注目すると、第一引数は単なるポインタではなく、情報がつめ込まれた何らかの構造体の先頭を指しており、その1つめの要素が配列の先頭ポインタであることがわかる。
struct FortranIntArray { Int* array; // rdiの指す先 // ここに32バイト分何かが入っている // が、今回は使われていない QWORD unknown; // rdi+40の指す先 };
このrdi+40が指しているものは添字の計算で使われているので、この辺で配列の形状ではないかと予想がつく。これがどのように使われているか追いかけてみよう。一部抜粋する。
; rdx = FortranIntArray.unknown; mov rdx, QWORD PTR [rdi+40] ; rax = 1 mov eax, 1 ; if (rdx == 0) {ZF = 1;} else {ZF = 0} test rdx, rdx ; if (ZF == 1) {rdx = rax;} cmove rdx, rax
この時点で、もしunknown == 0ならrdx = 1になっている。この仮定（unknown == 0）を置いたまま、最後まで読んでみよう。
nth_: ; rdx = 1になったとする。 cmove rdx, rax ; rax に第二引数の値（添字）を符号拡張して代入 movsx rax, DWORD PTR [rsi] ; rax *= 1 imul rax, rdx ; rax -= 1 sub rax, rdx ; rcxからrax*4バイト先の32bitをロードしてeaxに代入 mov eax, DWORD PTR [rcx+rax*4] ; return eax; ret
この場合、渡した添字の一つ前を見ていることになる。最初に思い出したとおり、Fortranでは配列の添字が1から始まる。ポインタは先頭を指しているので0から始まることになる。この不一致を解決するために、1を引く必要がある。
では、unknown != 0だったケースを考えてみよう。1だったらどうなるだろうか。rdxは非ゼロなので1のままになる。っと、これは0の場合と同じだ。0のときはraxが代入されて1になるので。
じゃあ例えば、unknown == 2だったら？
; rdxが非ゼロなので、rdx == 2のまま cmove rdx, rax ; rax に第二引数の値（添字）を符号拡張して代入 movsx rax, DWORD PTR [rsi] ; rax *= 2 imul rax, rdx ; rax -= 2 sub rax, rdx ; rcxからrax*4バイト先の32bitをロードしてeaxに代入 mov eax, DWORD PTR [rcx+rax*4] ; return eax; ret
添字が2倍されてから2を引かれた。1番めの要素の場合、オフセットは0。2番めの要素の場合、オフセットは2。3番めの要素の場合、オフセットは4。これは、第一カラムに何要素あるかがunknownに入っていて、その分をスキップしているように見える。Fortranはcolumn-majorなので、配列のメモリ上での配置は以下のようになるからだ（わかりやすさのため、添字を0スタートにしている）。
0, 2, 4 1, 3, 5
と、ここまでで大分このコードが何をしているのかわかったわけだが、それはそれとしてFortranは配列の次元が違っていたらコンパイルエラーにしてくれる。なので実際にはa(2,3)みたいな配列をnthに渡すことはできない。多次元配列を渡すときは、多次元配列であるという型情報が関数のインターフェースに載るのだ。だからこの関数にはかならず1次元配列が来る。そう考えるとこの処理には何の意味もないのではないかという気もする。
（追記）後で思い出したが、Fortranにはスライスの時にストライドを設定できる。a(2:6:2)のようにすると、a(2), a(4), a(6)だけが使用可能になる。こういう場合だと、一次元配列なのにオフセットが必要になることがあり得る。（追記終わり）
本当にオフセットが入っているのかわからなくなってきたので、無理やりCのコードをリンクして配列の先を見てみることにした。まずは一次元配列を見てみる。比較しやすいように、デフォルトのa(6)、3から始まるa(3:6)、0から始まるa(0:6)を見てみる。ここで、コンパイラを手元のものに変えたため、詳細なレイアウトは変わっているかも知れないことに注意したい。
integer :: a(6) [rdi+0] 7ffc9ef56ee0 [rdi+8] ffffffffffffffff [rdi+16] 109 [rdi+24] 1 [rdi+32] 1 [rdi+40] 6 integer :: a2(3:6) [rdi+0] 7ffc9ef56ed0 [rdi+8] fffffffffffffffd [rdi+16] 109 [rdi+24] 1 [rdi+32] 3 [rdi+40] 6 integer :: j3(0:6) [rdi+0] 7fff7273e770 [rdi+8] 0 [rdi+16] 109 [rdi+24] 1 [rdi+32] 0 [rdi+40] 6
最初の部分は先頭ポインタだろう。2つめの64bitは添字から引き算するべき値に見える。配列が1から始まるときは-1（2の補数）で、3から始まるときは-3、0から始まるときは0だからだ。3つめは謎。4つめもよくわからない。5つめと6つめは最初と最後の添字だろう。これはわかりやすい。
二次元配列はどうだろう。
integer :: a(5,6) [rdi+0] 7fff7273e810 [rdi+8] fffffffffffffffa [rdi+16] 10a [rdi+24] 1 [rdi+32] 1 [rdi+40] 5 [rdi+48] 5 [rdi+56] 1 [rdi+64] 6
心の目に1,1,5と5,1,6という塊が見える。これは1次元配列の場合、a(6)で1,1,6、a(3:6)で1,3,6、a(0:6)で1,0,6だったものに相当するのだろう。2次元のときに2つめの次元で5,1,6となっていて、最初の数字が5になっている。多分これは前までの次元の積で、この次元の配列を見に行こうとした時にスキップするべき塊のサイズを先に計算しているのではないか。
三次元配列を覗いてみる。a(5,6,7)としておけば、2次元めまでは同じ1,1,5と5,1,6で、3次元めで30,1,7になるだろう。
integer :: a(5,6,7) [rdi+0] 7fff7273e420 [rdi+8] ffffffffffffffdc [rdi+16] 10b [rdi+24] 1 [rdi+32] 1 [rdi+40] 5 [rdi+48] 5 [rdi+56] 1 [rdi+64] 6 [rdi+72] 1e [rdi+80] 1 [rdi+88] 7
30は16 + 14なので16進で1e。そうなっているようだ。ところでポインタの次にある値は-36になっている。これは6*5+6だ。5x6x7の3次元配列の(1,1,1)番目の要素は、もし配列が0始まりなら5x6x1 + 6x1 + 1になる。この前半部分、5x6x1+6x1を前もって計算してあるのだろう。
あとは、配列の次元をどうやって記録しているかが知りたいところだ。これまでの情報から計算はできるが、毎回計算していたのでは遅かろう。どこかに即値があると思いたい。
よく見ると、[rdi+16]に入っている値が次元を上げると共に1増えていることに気づく。これっぽい。なぜそんな大きな値になっているのかはまだわからないが。
長くなってしまった。ここまでのことを踏まえると、最初の関数には配列の長さ情報を与えておくと命令数が2命令くらいになるに違いない。
pure function nth(v, n) result(out) implicit none integer, intent(in) :: v(3) ! 例えば integer, intent(in) :: n integer :: out out = v(n) return end function nth
なった。
nth_:
movsx rax, DWORD PTR [rsi]
mov eax, DWORD PTR [rdi-4+rax*4]
ret
Fortranで関数に配列を渡すときは、もしわかっているのなら、配列の長さを陽に渡していたほうが速度面からよさそうに思える。
C++版と比較して少し不思議に思う所は、配列の添字をポインタで渡しているところだ。添字はintegerなので、明らかにレジスタに収まる。実際、C++で似たようなものを書くとnはesiに乗った状態で（メモリアクセスせずに）持ってこられている。だがFortran版では[rsi]でメモリアクセスをしている。
int nth_int(const int* v, const std::int32_t n) noexcept { return v[n]; }
nth_int(int const*, int): movsx rsi, esi mov eax, DWORD PTR [rdi+rsi*4] ret
ちょっと調べたところ、Fortran2003でVALUE属性というものが入ったらしい。これがあれば明示的に値渡しにできる。それ以外の場合は-O3でも参照渡しになるのはちょっとビビった。
試してみよう。
値渡しになったようだ。レジスタの値がそのまま使われている。
pure function nth(v, n) result(out) implicit none integer, intent(in) :: v(3) integer, value, intent(in) :: n integer :: out out = v(n) return end function nth
nth_:
movsx rsi, esi
mov eax, DWORD PTR [rdi-4+rsi*4]
ret
細かなところだが意外と複雑で面白かった。わざわざFortranを書く状況なら速さが正義なはずなので、value属性や配列の長さに気をつけたほうがいいでしょうね。
