背景
C++では1.0とか書くとdoubleとして解釈され、floatにするには1.0fなどと書かねばならないことはご承知のとおりだ。また、doubleとfloatの演算の結果はdoubleになる。
これはtemplateを使うときに若干面倒になる。例えば、以下のような関数fを実装したとする。
template<typename T> T f(T x, T y) { return x + y; }
さて、浮動小数点数型を受け取るクラステンプレートがあり、中でfを呼ぶとしよう。
template<typename T> struct X { T g(T x) {return f(x, 2.0 * x);} }; int main() { X<float> x; std::cout << x.g(1.0f) << std::endl; return 0; }
何が起きるだろうか。
答：コンパイルエラーになる。
なぜかというと、2.0 * xはdoubleとfloatの積なのでdoubleになり、するとtemplate型推論が上手く行かない。というのも、
template<typename T> T f(T x, T y)
に対して
f(float, double);
が与えられるのだ。Tはどちらとも取れない。
template argument deduction/substitution failed:
note: deduced conflicting types for parameter 'T' ('float' and 'double')
T g(T x) {return f(x, 2.0 * x);}
~^~~~~~~~~~~~
という感じになる。
templateクラスの中で書かれてるリテラルはT型に勝手になってくれよ〜と思うかも知れないが、T型とリテラルが同じ型になるという関係はこのコードを書いた人間の脳みその中にしかないので、コンパイラはそんな勝手な解釈をすることができない。
template<typename T> struct X { T g(T x) {return f(x, 2 * x);} };
これなら通る。2は整数なので、演算結果は浮動小数点数になる。
問題は整数では表現出来ないときで、これは難しい。例えば、倍にするのでなく2.5倍したいときはどうするか。
ナイーブにやると以下のようになる。
template<typename T> struct X { T g(T x) {return f(x, 2.5 * x);} };
これだと2.0の時と同じ問題が発生する。
それを避けるため、明示的にキャストするという手はある。
template<typename T> struct X { T g(T x) {return f(x, static_cast<T>(2.5) * x);} };
単純に長い。一つならまだ良いが、複数個こんなものを書くとしんどい。static_castでなくT(2.5)と書いてもまあ良いのだが……
template引数の推論に失敗するのだから、それを明示的に与えればよい。
template<typename T> struct X { T g(T x) {return f<T>(x, 2.5 * x);} };
これはまあすっきりしている。
他には、定数を集めたクラスを作っておいて、そこから持ってくるという戦略もあり得る。
template<typename T> struct constant { static constexpr T two_and_half = static_cast<T>(2.5); }; template<typename T> struct X { T g(T x) {return f(x, constant<T>::two_and_half * x);} };
果たしてこれでどれほど簡単になるのかという話だが。長くなると面倒な感じになっていくし、非常に冗長だ。piとかroot_piみたいな特殊な値を持っておくなら、この戦略はありなのだが（boost::math::constantsなどは似たようなことをしている）。
この例と同様のことは、例えば以下のような関数でも生じる。
template<typename T> std::complex<T> f(T a, std::complex<T> b);
これをfloatで実体化しつつ、aに2.0、bにstd::complex<float>とかを入れた場合、こうなる。
prog.cc: In function 'int main()':
prog.cc:11:50: error: no matching function for call to 'h(double, std::complex<float>)'
auto c = h(2.0, std::complex<float>(1.0, 1.0));
^
prog.cc:4:17: note: candidate: 'template<class T> std::complex<_Tp> h(T, std::complex<_Tp>)'
std::complex<T> h(T a, std::complex<T> b)
^
prog.cc:4:17: note: template argument deduction/substitution failed:
prog.cc:11:50: note: deduced conflicting types for parameter '_Tp' ('double' and 'float')
auto c = h(2.0, std::complex<float>(1.0, 1.0));
^
エラーの内容は同じで、型推論の結果doubleかfloatかわからなかったよ、という話だ。
速度の話
と、さてここまでで大体めんどくさそうだなあということはわかっていただけたと思う。
面倒くさいだけならまだ良いのだが。以下のリンク先を見て欲しい。
このリンク先はwandboxと同じくオンラインコンパイラなのだが、アセンブリを見せてくれるところが違いだ。アセンブリとC++コードの対応する箇所をハイライトしてくれたり、色々便利だ。
そこで、以下のようなコードをコンパイルした。
template<typename real> real f(real x) { const real cos1 = std::cos(x); const real cos3 = cos1 * (4.0 * cos1 * cos1 - 3.0); return cos1 - cos3; }
同じことをする関数を3つ書き、それぞれ、リテラルを4.0、4.0f、4と変えていって、double、floatで実体化した。ここに載せているのはgcc-8.2.0 -std=c++11 -O3 -ffast-mathの結果だが、clang-6.0.0でも基本的に同じことが起きた。
特筆すべきは、doubleリテラルを使ってfloatで実体化したときだ。コンパイラフラグは-ffast-mathと-O3を立てているので、リテラルはfloatに再解釈されることを期待していた。しかし、現実はこうだ。
float f<float>(float): sub rsp, 8 call cosf pxor xmm2, xmm2 cvtss2sd xmm2, xmm0 movapd xmm1, xmm2 mulsd xmm1, xmm2 mulsd xmm1, QWORD PTR .LC1[rip] subsd xmm1, QWORD PTR .LC2[rip] add rsp, 8 mulsd xmm1, xmm2 cvtsd2ss xmm1, xmm1 subss xmm0, xmm1 ret
ところで、cvtss2sdは"convert single single-precision floating point to single double-precision floating point"の略である。要するにfloatをdoubleに変換するもので（singleが2回書かれているのは、SIMDという複数個一気に計算する命令があるから、それとの区別のためである）、cvsd2ssはその逆だ。
もう何が起きているかはわかるだろうが、もう少し追っておくと、mulsdは"multiply single double-precision floating point"であり、subsdは"subtract single double-precision floating point"の略だ。何が起きているかというと、こいつはfloatで来た引数を一度doubleにしてから計算して後でfloatに戻しているのである。
断っておくが、これはよくある「C言語ではfloatは必ず一度doubleにされてから計算されるからdoubleの方が速い」系の話の亜種ではない。これは現代においては殆どの場合嘘だ。K&R時代とか私は生まれてないぞ。現代ではちゃんと単精度用の命令があるので、floatの計算はmulssやらsubssやらで行われ、実行速度に特に差はない。SIMDとかなら差は出るが、その場合floatの方が速かったりする。
その辺りは詳しくは以下の記事などを参照していただきたい。
さて、対比のために単精度浮動小数点数リテラルを使った場合のアセンブリを見てみよう。
float g<float>(float): sub rsp, 8 call cosf movss xmm1, DWORD PTR .LC3[rip] movaps xmm2, xmm0 mulss xmm2, xmm0 mulss xmm0, DWORD PTR .LC4[rip] add rsp, 8 subss xmm1, xmm2 mulss xmm0, xmm1 ret
変換命令はないし、subssやmulssが使われている。やはり、倍精度リテラルを使ったことで問題が生じたのだと考えるのが自然だろう。
