FFT: std::complex の代替は必要か？

浮動小数点演算を使うプログラムの多くは、無限大や非数を使用しない前提でしょう。

FFT を実装するとき、複素数型 std::complex を使いたくなります。 std::complex を使ったものを、コンパイルしてアセンブリ出力を確認すると、数命令での展開を期待できる部分が（コードが大きいので）関数呼び出しになったりします。

std::complex を使ったサンプル プログラム

```c++:sample.cpp
#include <complex>

typedef std::complex<float> complex_t;

complex_t sample(complex_t &a, complex_t &b)
{
  return a * b;
}
```

Apple clang version 11.0.0 (clang-1100.0.33.17)
Target: x86_64-apple-darwin19.3.0
$ g++ -O3 -S sample.cpp

<details><summary>のアセンブリ出力（長いので折りたたみ）</summary><div>


```text:sample.s
	.section	__TEXT,__text,regular,pure_instructions
	.build_version macos, 10, 15	sdk_version 10, 15
	.globl	__Z6sampleRNSt3__17complexIfEES2_ ## -- Begin function _Z6sampleRNSt3__17complexIfEES2_
	.p2align	4, 0x90
__Z6sampleRNSt3__17complexIfEES2_:      ## @_Z6sampleRNSt3__17complexIfEES2_
	.cfi_startproc
## %bb.0:
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset %rbp, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register %rbp
	popq	%rbp
	jmp	__ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_ ## TAILCALL
	.cfi_endproc
                                        ## -- End function
	.section	__TEXT,__literal16,16byte_literals
	.p2align	4               ## -- Begin function _ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_
LCPI1_0:
	.long	2147483647              ## float NaN
	.long	2147483647              ## float NaN
	.long	2147483647              ## float NaN
	.long	2147483647              ## float NaN
LCPI1_3:
	.long	2147483648              ## float -0
	.long	2147483648              ## float -0
	.long	2147483648              ## float -0
	.long	2147483648              ## float -0
LCPI1_4:
	.long	2139095040              ## float +Inf
	.long	2139095040              ## float +Inf
	.long	2139095040              ## float +Inf
	.long	2139095040              ## float +Inf
	.section	__TEXT,__literal4,4byte_literals
	.p2align	2
LCPI1_1:
	.long	2139095040              ## float +Inf
LCPI1_2:
	.long	1065353216              ## float 1
	.section	__TEXT,__text,regular,pure_instructions
	.globl	__ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_
	.weak_def_can_be_hidden	__ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_
	.p2align	4, 0x90
__ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_: ## @_ZNSt3__1mlIfEENS_7complexIT_EERKS3_S5_
	.cfi_startproc
## %bb.0:
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset %rbp, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register %rbp
	movsd	(%rdi), %xmm4           ## xmm4 = mem[0],zero
	movaps	%xmm4, %xmm0
	shufps	$65, %xmm4, %xmm0       ## xmm0 = xmm0[1,0],xmm4[0,1]
	movsd	(%rsi), %xmm10          ## xmm10 = mem[0],zero
	movddup	%xmm10, %xmm6           ## xmm6 = xmm10[0,0]
	mulps	%xmm0, %xmm6
	movaps	%xmm6, %xmm0
	unpckhpd	%xmm6, %xmm0    ## xmm0 = xmm0[1],xmm6[1]
	movaps	%xmm6, %xmm1
	shufps	$231, %xmm6, %xmm1      ## xmm1 = xmm1[3,1],xmm6[2,3]
	subss	%xmm1, %xmm0
	movshdup	%xmm6, %xmm5    ## xmm5 = xmm6[1,1,3,3]
	addss	%xmm6, %xmm5
	ucomiss	%xmm0, %xmm0
	jnp	LBB1_16
## %bb.1:
	ucomiss	%xmm5, %xmm5
	jnp	LBB1_16
## %bb.2:
	movaps	LCPI1_0(%rip), %xmm8    ## xmm8 = [NaN,NaN,NaN,NaN]
	movaps	%xmm4, %xmm1
	andps	%xmm8, %xmm1
	ucomiss	LCPI1_1(%rip), %xmm1
	movshdup	%xmm4, %xmm11   ## xmm11 = xmm4[1,1,3,3]
	movaps	%xmm11, %xmm3
	andps	%xmm8, %xmm3
	movshdup	%xmm10, %xmm9   ## xmm9 = xmm10[1,1,3,3]
	jae	LBB1_5
## %bb.3:
	ucomiss	LCPI1_1(%rip), %xmm3
	jae	LBB1_5
## %bb.4:
	xorl	%eax, %eax
	jmp	LBB1_7
LBB1_5:
	movss	LCPI1_1(%rip), %xmm7    ## xmm7 = mem[0],zero,zero,zero
	movaps	%xmm7, %xmm2
	cmpless	%xmm1, %xmm2
	movss	LCPI1_2(%rip), %xmm12   ## xmm12 = mem[0],zero,zero,zero
	andps	%xmm12, %xmm2
	movaps	LCPI1_3(%rip), %xmm1    ## xmm1 = [-0.0E+0,-0.0E+0,-0.0E+0,-0.0E+0]
	andps	%xmm1, %xmm4
	andps	%xmm8, %xmm2
	orps	%xmm2, %xmm4
	cmpless	%xmm3, %xmm7
	andps	%xmm12, %xmm7
	andps	%xmm1, %xmm11
	andps	%xmm8, %xmm7
	orps	%xmm7, %xmm11
	andps	%xmm10, %xmm1
	movaps	%xmm10, %xmm2
	cmpunordss	%xmm10, %xmm2
	movaps	%xmm2, %xmm3
	andnps	%xmm10, %xmm3
	andps	%xmm1, %xmm2
	orps	%xmm3, %xmm2
	movaps	%xmm2, %xmm10
	movb	$1, %al
	ucomiss	%xmm9, %xmm9
	jnp	LBB1_7
## %bb.6:
	movaps	LCPI1_3(%rip), %xmm1    ## xmm1 = [-0.0E+0,-0.0E+0,-0.0E+0,-0.0E+0]
	andps	%xmm1, %xmm9
LBB1_7:
	movaps	%xmm10, %xmm1
	andps	%xmm8, %xmm1
	movaps	%xmm9, %xmm3
	andps	%xmm8, %xmm3
	ucomiss	LCPI1_1(%rip), %xmm3
	jae	LBB1_9
## %bb.8:
	ucomiss	LCPI1_1(%rip), %xmm1
	jae	LBB1_9
## %bb.11:
	testb	%al, %al
	jne	LBB1_15
## %bb.12:
	andps	LCPI1_0(%rip), %xmm6
	cmpeqps	LCPI1_4(%rip), %xmm6
	pshufd	$78, %xmm6, %xmm1       ## xmm1 = xmm6[2,3,0,1]
	por	%xmm6, %xmm1
	pshufd	$229, %xmm1, %xmm2      ## xmm2 = xmm1[1,1,2,3]
	por	%xmm1, %xmm2
	pextrb	$0, %xmm2, %eax
	testb	$1, %al
	je	LBB1_16
## %bb.13:
	movaps	LCPI1_3(%rip), %xmm0    ## xmm0 = [-0.0E+0,-0.0E+0,-0.0E+0,-0.0E+0]
	movaps	%xmm4, %xmm1
	cmpunordss	%xmm4, %xmm1
	movaps	%xmm1, %xmm2
	andnps	%xmm4, %xmm2
	andps	%xmm0, %xmm4
	andps	%xmm4, %xmm1
	orps	%xmm2, %xmm1
	movaps	%xmm1, %xmm4
	movaps	%xmm11, %xmm1
	cmpunordss	%xmm11, %xmm1
	movaps	%xmm1, %xmm2
	andnps	%xmm11, %xmm2
	andps	%xmm0, %xmm11
	andps	%xmm11, %xmm1
	orps	%xmm2, %xmm1
	movaps	%xmm1, %xmm11
	andps	%xmm10, %xmm0
	movaps	%xmm10, %xmm1
	cmpunordss	%xmm10, %xmm1
	movaps	%xmm1, %xmm2
	andnps	%xmm10, %xmm2
	andps	%xmm0, %xmm1
	orps	%xmm2, %xmm1
	movaps	%xmm1, %xmm10
	ucomiss	%xmm9, %xmm9
	jnp	LBB1_15
## %bb.14:
	andps	LCPI1_3(%rip), %xmm9
	jmp	LBB1_15
LBB1_9:
	movss	LCPI1_1(%rip), %xmm0    ## xmm0 = mem[0],zero,zero,zero
	movaps	%xmm0, %xmm2
	cmpless	%xmm1, %xmm2
	movss	LCPI1_2(%rip), %xmm1    ## xmm1 = mem[0],zero,zero,zero
	andps	%xmm1, %xmm2
	movaps	LCPI1_3(%rip), %xmm5    ## xmm5 = [-0.0E+0,-0.0E+0,-0.0E+0,-0.0E+0]
	andps	%xmm5, %xmm10
	andps	%xmm8, %xmm2
	orps	%xmm2, %xmm10
	cmpless	%xmm3, %xmm0
	andps	%xmm1, %xmm0
	andps	%xmm5, %xmm9
	andps	%xmm8, %xmm0
	orps	%xmm0, %xmm9
	andps	%xmm4, %xmm5
	movaps	%xmm4, %xmm0
	cmpunordss	%xmm4, %xmm0
	movaps	%xmm0, %xmm1
	andnps	%xmm4, %xmm1
	andps	%xmm5, %xmm0
	orps	%xmm1, %xmm0
	movaps	%xmm0, %xmm4
	ucomiss	%xmm11, %xmm11
	jnp	LBB1_15
## %bb.10:
	andps	LCPI1_3(%rip), %xmm11
LBB1_15:
	movaps	%xmm4, %xmm0
	mulss	%xmm10, %xmm0
	movaps	%xmm11, %xmm1
	mulss	%xmm9, %xmm1
	subss	%xmm1, %xmm0
	movss	LCPI1_1(%rip), %xmm1    ## xmm1 = mem[0],zero,zero,zero
	mulss	%xmm1, %xmm0
	mulss	%xmm9, %xmm4
	mulss	%xmm10, %xmm11
	addss	%xmm4, %xmm11
	mulss	%xmm1, %xmm11
	movaps	%xmm11, %xmm5
LBB1_16:
	insertps	$16, %xmm5, %xmm0 ## xmm0 = xmm0[0],xmm5[0],xmm0[2,3]
	popq	%rbp
	retq
	.cfi_endproc
                                        ## -- End function

.subsections_via_symbols
```

</div></details>

では、無限大や非数の処理が沢山入っています。

g++ の最適化オプション -ffast-math を付けた

$ g++ -O3 -ffast-math -S sample.cpp

のアセンブリ出力は

```text:sample.s
	.section	__TEXT,__text,regular,pure_instructions
	.build_version macos, 10, 15	sdk_version 10, 15
	.globl	__Z6sampleRNSt3__17complexIfEES2_ ## -- Begin function _Z6sampleRNSt3__17complexIfEES2_
	.p2align	4, 0x90
__Z6sampleRNSt3__17complexIfEES2_:      ## @_Z6sampleRNSt3__17complexIfEES2_
	.cfi_startproc
## %bb.0:
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset %rbp, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register %rbp
	movsd	(%rdi), %xmm1           ## xmm1 = mem[0],zero
	movsd	(%rsi), %xmm0           ## xmm0 = mem[0],zero
	movshdup	%xmm0, %xmm2    ## xmm2 = xmm0[1,1,3,3]
	mulss	%xmm1, %xmm2
	movshdup	%xmm1, %xmm3    ## xmm3 = xmm1[1,1,3,3]
	mulss	%xmm0, %xmm3
	addss	%xmm2, %xmm3
	mulps	%xmm1, %xmm0
	movshdup	%xmm0, %xmm1    ## xmm1 = xmm0[1,1,3,3]
	subss	%xmm1, %xmm0
	insertps	$16, %xmm3, %xmm0 ## xmm0 = xmm0[0],xmm3[0],xmm0[2,3]
	popq	%rbp
	retq
	.cfi_endproc
                                        ## -- End function

.subsections_via_symbols
```

無限大や非数のチェックや分岐などがなくなります。

既定のコンパイル オプションで使う場合は、std::complex の代替を用意したほうがよさそうですが、 -ffast-math で十分な場合も少なくないでしょう。

VisualStudio (VC++) では /fp:mode オプションがありますが、std::complex の実装が違うためか g++ のような大きなプログラムは Release ビルドでは確認できませんでした。