C言語とFortranの配列アクセスはどちらが速いか?

「Fortranの配列アクセスはC言語の配列アクセスよりも速い」
という噂を巷で聞いたので、本当にそうなのか調べてみました。

# 検証方法

以下の2次元配列を初期化するコードを用意し、`gcc`と`gfortran`でアセンブルします。

```c.c
int main(void) {
        int m[5][5], i, j;
        for (i = 0; i < 5; i++) {
                for (j = 0; j < 5; j++) {
                        m[i][j] = i * j;
                }
        }
}
```

```f.f90
program f
    integer m(5,5), i, j
    do i = 1, 5
        do j = 1, 5
            m(j, i) = i * j
        end do
    end do
end program f
```

なお、今回gcc 8.2.1を用いました。

# 逆アセンブル

どちらも以下のように最適化オプションを付けずにコンパイルしました。というのも、`gcc`と`gfortran`のバックエンドは共通であり、最適化後に生じうる差は最適化前の差に起因するものであると考えられるからです。(それに最適化前の方が両言語の設計思想の差も確認できて良いですね)

```terminal
$ gcc -S c.c
$ gfortran -S f.f90
```

逆アセンブル結果は以下のようになります。(一部省略)

```c.s
	.globl	main
main:
	pushq	%rbp
	movq	%rsp, %rbp
	addq	$-128, %rsp
	movq	%fs:40, %rax
	movq	%rax, -8(%rbp)
	xorl	%eax, %eax
	movl	$0, -120(%rbp)
	jmp	.L2
.L5:
	movl	$0, -116(%rbp)
	jmp	.L3
.L4:
	movl	-120(%rbp), %eax
	imull	-116(%rbp), %eax
	movl	%eax, %ecx
	movl	-116(%rbp), %eax
	movslq	%eax, %rsi
	movl	-120(%rbp), %eax
	movslq	%eax, %rdx
	movq	%rdx, %rax
	salq	$2, %rax
	addq	%rdx, %rax
	addq	%rsi, %rax
	movl	%ecx, -112(%rbp,%rax,4)
	addl	$1, -116(%rbp)
.L3:
	cmpl	$4, -116(%rbp)
	jle	.L4
	addl	$1, -120(%rbp)
.L2:
	cmpl	$4, -120(%rbp)
	jle	.L5
	movl	$0, %eax
	movq	-8(%rbp), %rdi
	xorq	%fs:40, %rdi
	je	.L7
	call	__stack_chk_fail@PLT
.L7:
	leave
	ret
```

```f.s
MAIN__:
	pushq	%rbp
	movq	%rsp, %rbp
	addq	$-128, %rsp
	movq	%fs:40, %rax
	movq	%rax, -8(%rbp)
	xorl	%eax, %eax
	movl	$1, -120(%rbp)
.L5:
	cmpl	$5, -120(%rbp)
	setg	%al
	movzbl	%al, %eax
	testl	%eax, %eax
	jne	.L7
	movl	$1, -116(%rbp)
.L4:
	cmpl	$5, -116(%rbp)
	setg	%al
	movzbl	%al, %eax
	testl	%eax, %eax
	jne	.L8
	movl	-120(%rbp), %eax
	movslq	%eax, %rdx
	movq	%rdx, %rax
	salq	$2, %rax
	addq	%rax, %rdx
	movl	-116(%rbp), %eax
	cltq
	addq	%rdx, %rax
	leaq	-6(%rax), %rdx
	movl	-120(%rbp), %eax
	imull	-116(%rbp), %eax
	movl	%eax, -112(%rbp,%rdx,4)
	addl	$1, -116(%rbp)
	jmp	.L4
.L8:
	nop
	addl	$1, -120(%rbp)
	jmp	.L5
.L7:
	nop
	nop
	movq	-8(%rbp), %rax
	xorq	%fs:40, %rax
	je	.L6
	call	__stack_chk_fail@PLT
.L6:
	leave
	ret
	.globl	main
main:
	pushq	%rbp
	movq	%rsp, %rbp
	subq	$16, %rsp
	movl	%edi, -4(%rbp)
	movq	%rsi, -16(%rbp)
	movq	-16(%rbp), %rdx
	movl	-4(%rbp), %eax
	movq	%rdx, %rsi
	movl	%eax, %edi
	call	_gfortran_set_args@PLT
	leaq	options.0.3784(%rip), %rsi
	movl	$7, %edi
	call	_gfortran_set_options@PLT
	call	MAIN__
	movl	$0, %eax
	leave
	ret
options.0.3784:
	.long	68
	.long	8191
	.long	0
	.long	1
	.long	1
	.long	0
	.long	31
```


# ループ部のアセンブリコードの解読

## C言語の場合

C言語の外側のループに対応するアセンブリコードの概略は以下のようになります。

```
	movl	$0, -120(%rbp)
	jmp	.L2
.L5:
# ...
	addl	$1, -120(%rbp)
.L2:
	cmpl	$4, -120(%rbp)
	jle	.L5
```

`i<5`の条件判定をループの中身の処理のあと(`.L2`)の部分で行っている点に注目。
これによって、ループの中で5回実行される分岐命令は `jle .L5` の1箇所のみとなっています。

## Fortranの場合

```
	movl	$1, -120(%rbp)
.L5:
	cmpl	$5, -120(%rbp)
	setg	%al
	movzbl	%al, %eax
	testl	%eax, %eax
	jne	.L7
# ...
	addl	$1, -120(%rbp)
	jmp	.L5
.L7:
	nop
```

`i`の比較処理(`.L5`)の部分がC言語と比べて複雑になっています。

```
    cmpl    $5, -120(%rbp)  # iの値を5と比較
    setg    %al             # ZFがクリアされていたら1を%alにセット
                            # (i != 5 なら%al == 1となる)
    movzbl  %al, %eax       # %alの値を%eaxに移す
    testl   %eax, %eax      # %eax の値を確認
    jne .L7                 # %eax == 0なら分岐
```

なぜこのような複雑なことをしているかは謎です。

# 演算`i * j`実行部

## C言語の場合

```
	movl	-120(%rbp), %eax
	imull	-116(%rbp), %eax
	movl	%eax, %ecx
```

# Fortranの場合
```
	movl	-120(%rbp), %eax
	imull	-116(%rbp), %eax
```
Fortranでは、先に代入先の配列要素のアドレスを計算しているため、`%eax`を`%ecx`に退避する必要がない分短くなっていますが、このような違いは最適化によって速攻解消されるので普通はあまり関係ないです。

# 配列要素への代入部

C言語とFortranでほぼ同一です。

```
	movl	%eax, -112(%rbp,%rax,4)
```

上のコードをC言語風に書くと以下のようになります。

```c
*(rbp - 112 + rax * 4) = eax
```

`rbp-112`は配列`m`のアドレスであると考えられます。
そこで、この`rax`をどのように計算するかが次の話題となります。

# アドレスの計算

## C言語の場合

```
	movl	-116(%rbp), %eax # j -> %eax
	movslq	%eax, %rsi
	movl	-120(%rbp), %eax # i -> %eax
	movslq	%eax, %rdx
	movq	%rdx, %rax
	salq	$2, %rax
	addq	%rdx, %rax
	addq	%rsi, %rax
```

このコードを実行すると、`%rax`には`j + 5 * i`を計算した結果が入ります。すなわち、多次元配列に対するアクセス`m[i][j]`が、1次元配列へのアクセス`m[j + 5 * i]`のように読み替えられています。

## Fortranの場合

```
	movl	-120(%rbp), %eax      # i -> %eax
	movslq	%eax, %rdx
	movq	%rdx, %rax
	salq	$2, %rax
	addq	%rax, %rdx
	movl	-116(%rbp), %eax     # j -> %eax
	cltq
	addq	%rdx, %rax       # rax = 5 * i + j
	leaq	-6(%rax), %rdx   # %rdx = %rax - 6
```

Fortranでは、C言語と同様にアドレス計算を行っているものの、最後に6要素分アドレスを減算しています。
これはFortranの添字が`1`から始まることに起因しています。

# 結論

- C言語のアドレッシングがFortranに比べて遅いという事はない
- Fortranでは添字が`1`から始まるため、その分アドレス計算の処理が余分に必要
- 実際のユースケースでは最適化を行うため、今回の比較がどれほど実用面での影響を及ぼすかは不明

次回、配列の初期化編(好評だったらやりたい)