はじめに
既に気付いている人も多いかもしれないが、std::complexの積和の計算を約20%高速化する書き方に気付いたのでメモしておく。今担当している組込リアルタイムDSPの時間制約がキツく、僅かでも高速化したいと思っていた。この方法は手軽さの割に効果が大きい。
\[
% 汎用
\newcommand{\ctext}[1]{\raise0.2ex\hbox{\textcircled{\scriptsize{#1}}}}
% 数学
% 汎用
\DeclarePairedDelimiterX{\parens}[1]{\lparen}{\rparen}{#1}
\DeclarePairedDelimiterX{\braces}[1]{\lbrace}{\rbrace}{#1}
\DeclarePairedDelimiterX{\bracks}[1]{\lbrack}{\rbrack}{#1}
\DeclarePairedDelimiterX{\verts}[1]{|}{|}{#1}
\DeclarePairedDelimiterX{\Verts}[1]{\|}{\|}{#1}
\newcommand{\as}{{\quad\textrm{as}\quad}}
\newcommand{\st}{{\textrm{ s.t. }}}
\DeclarePairedDelimiterX{\setComprehension}[2]{\lbrace}{\rbrace}{#1\,\delimsize\vert\,#2}
\newcommand{\naturalNumbers}{\mathbb{N}}
\newcommand{\integers}{\mathbb{Z}}
\newcommand{\rationalNumbers}{\mathbb{Q}}
\newcommand{\realNumbers}{\mathbb{R}}
\newcommand{\complexNumbers}{\mathbb{C}}
\newcommand{\field}{\mathbb{F}}
\newcommand{\func}[2]{{#1}\parens*{#2}}
\newcommand{\argmax}{\mathop{\textrm{arg~max}}}
\newcommand{\argmin}{\mathop{\textrm{arg~min}}}
% 集合論
\newcommand{\range}[2]{\braces*{#1,\dotsc,#2}}
\providecommand{\complement}{}\renewcommand{\complement}{\mathrm{c}}
\newcommand{\ind}[2]{\mathbbm{1}_{#1}\parens*{#2}}
\newcommand{\indII}[1]{\mathbbm{1}\braces*{#1}}
% 数論
\newcommand{\abs}[1]{\verts*{#1}}
\newcommand{\combi}[2]{{_{#1}\mathrm{C}_{#2}}}
\newcommand{\perm}[2]{{_{#1}\mathrm{P}_{#2}}}
\newcommand{\GaloisField}[1]{\mathrm{GF}\parens*{#1}}
% 解析学
\newcommand{\NapierE}{\mathrm{e}}
\newcommand{\sgn}[1]{\operatorname{sgn}\parens*{#1}}
\newcommand{\cl}[1]{\operatorname{cl}#1}
\newcommand{\Img}[1]{\operatorname{Img}\parens*{#1}}
\newcommand{\dom}[1]{\operatorname{dom}\parens*{#1}}
\newcommand{\norm}[1]{\Verts*{#1}}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\ceil}[1]{\left\lceil#1\right\rceil}
\newcommand{\expo}[1]{\exp\parens*{#1}}
\newcommand{\sinc}{\mathop{\textrm{sinc}}}
\newcommand{\GammaFunc}[1]{\Gamma\parens*{#1}}
\newcommand{\erf}{\mathop{\mathrm{erf}}}
% 逆三角関数
\newcommand{\asin}[1]{\operatorname{Sin}^{-1}{#1}}
\newcommand{\acos}[1]{\operatorname{Cos}^{-1}{#1}}
\newcommand{\atan}[1]{\operatorname{{Tan}^{-1}}{#1}}
\newcommand{\atanEx}[2]{\atan{\parens*{#1,#2}}}
% 畳み込み
\newcommand{\cycConv}[2]{{#1}\underset{\text{cyc}}{*}{#2}}
% 微分
\newcommand{\deriv}[3]{\frac{\operatorname{d}^{#3}#1}{\operatorname{d}{#2}^{#3}}}
\newcommand{\derivLong}[3]{\frac{\operatorname{d}^{#3}}{\operatorname{d}{#2}^{#3}}#1}
\newcommand{\partDeriv}[3]{\frac{\operatorname{\partial}^{#3}#1}{\operatorname{\partial}{#2}^{#3}}}
\newcommand{\partDerivLong}[3]{\frac{\operatorname{\partial}^{#3}}{\operatorname{\partial}{#2}^{#3}}#1}
\newcommand{\partDerivIIHetero}[3]{\frac{\operatorname{\partial}^2#1}{\partial#2\operatorname{\partial}#3}}
\newcommand{\partDerivIIHeteroLong}[3]{{\frac{\operatorname{\partial}^2}{\partial#2\operatorname{\partial}#3}#1}}
% 積分
\newcommand{\integrate}[5]{\int_{#1}^{#2}{#3}{\mathrm{d}^{#4}}#5}
\newcommand{\LebInteg}[4]{\int_{#1} {#2} {#3}\parens*{\mathrm{d}#4}}
% 複素解析
\newcommand{\conj}[1]{\overline{#1}}
\providecommand{\Re}{}\renewcommand{\Re}[1]{{\operatorname{Re}{\parens*{#1}}}}
\providecommand{\Im}{}\renewcommand{\Im}[1]{{\operatorname{Im}{\parens*{#1}}}}
\newcommand{\Arg}[1]{\operatorname{Arg}{\parens*{#1}}}
\newcommand{\Log}[1]{\operatorname{Log}{#1}}
% ラプラス変換
\newcommand{\LPLC}[1]{\operatorname{\mathcal{L}}\parens*{#1}}
\newcommand{\ILPLC}[1]{\operatorname{\mathcal{L}}^{-1}\parens*{#1}}
% 離散Fourier変換
\newcommand{\DFT}[1]{\mathrm{DFT}\parens*{#1}}
% Z変換
\newcommand{\ZTrans}[1]{\operatorname{\mathcal{Z}}\parens*{#1}}
\newcommand{\IZTrans}[1]{\operatorname{\mathcal{Z}}^{-1}\parens*{#1}}
% 線形代数
\newcommand{\bm}[1]{{\boldsymbol{#1}}}
\newcommand{\matEntry}[3]{#1\bracks*{#2}\bracks*{#3}}
\newcommand{\matPart}[5]{\matEntry{#1}{#2:#3}{#4:#5}}
\newcommand{\diag}[1]{\operatorname{diag}\parens*{#1}}
\newcommand{\tr}[1]{\operatorname{tr}{\parens*{#1}}}
\newcommand{\inprod}[2]{\left\langle#1,#2\right\rangle}
\newcommand{\HadamardProd}{\odot}
\newcommand{\HadamardDiv}{\oslash}
\newcommand{\Span}[1]{\operatorname{span}\bracks*{#1}}
\newcommand{\Ker}[1]{\operatorname{Ker}\parens*{#1}}
\newcommand{\rank}[1]{\operatorname{rank}\parens*{#1}}
% ベクトル
% 単位ベクトル
\newcommand{\vix}{\bm{i}_x}
\newcommand{\viy}{\bm{i}_y}
\newcommand{\viz}{\bm{i}_z}
% 確率論
\newcommand{\PDF}[2]{\operatorname{PDF}\bracks*{#1,\;#2}}
\newcommand{\Ber}[1]{\operatorname{Ber}\parens*{#1}}
\newcommand{\Beta}[2]{\operatorname{Beta}\parens*{#1,#2}}
\newcommand{\ExpDist}[1]{\operatorname{ExpDist}\parens*{#1}}
\newcommand{\ErlangDist}[2]{\operatorname{ErlangDist}\parens*{#1,#2}}
\newcommand{\PoissonDist}[1]{\operatorname{PoissonDist}\parens*{#1}}
\newcommand{\GammaDist}[2]{\operatorname{Gamma}\parens*{#1,#2}}
\newcommand{\cind}[2]{\ind{#1\left| #2\right.}} %条件付き指示関数
\providecommand{\Pr}{}\renewcommand{\Pr}[1]{\operatorname{Pr}\parens*{#1}}
\newcommand{\cPr}[2]{\operatorname{Pr}\cPrParens{#1}{#2}}
\newcommand{\E}[2]{\operatorname{E}_{#1}\bracks*{#2}}
\newcommand{\cE}[3]{\E{#1}{\left.#2\right|#3}}
\newcommand{\Var}[2]{\operatorname{Var}_{#1}\bracks*{#2}}
\newcommand{\Cov}[2]{\operatorname{Cov}\bracks*{#1,#2}}
\newcommand{\CovMat}[1]{\operatorname{Cov}\bracks*{#1}}
% グラフ理論
\newcommand{\neighborhood}{\mathcal{N}}
% プログラミング
\newcommand{\plpl}{\mathrel{++}}
\newcommand{\pleq}{\mathrel{+}=}
\newcommand{\asteq}{\mathrel{*}=}
\]
問題と解決法
std::complexの積和(例えばベクトルの内積)を計算するとき、普通は次のように書くだろう(x,y等の変数は適当に定義されているものとする)。
std::complex<float> sum = 0;
for (int i=0; i<N; ++i) {
sum += x[i]*y[i];
}これだとx[i]*y[i]の結果が一時オブジェクトとして生成され、それがsumに加算されるという形をとり、一時オブジェクトの生成でstd::complexのコンストラクタが走るのでコストが掛かる。
しかし、std::complexは実部と虚部がメモリ上で隣接して配置されるという規格を利用し、浮動小数点数の配列(長さ2)として中身を直接扱うことで一時オブジェクトの生成を省略できる。
std::complex<float> sum = 0;
for (int i=0; i<N_sim; ++i) {
float *const sum_ptr = reinterpret_cast<float *>(sum);
const float *const x_ptr = reinterpret_cast<const float *>(&x[i]);
const float *const y_ptr = reinterpret_cast<const float *>(&y[i]);
sum[0] += x_ptr[0]*y_ptr[0] - x_ptr[1]*y_ptr[1]; // real part
sum[1] += x_ptr[0]*y_ptr[1] + x_ptr[1]*y_ptr[0]; // imaginary part
}ループ中でのポインタ定義の実行時コストは次に挙げる理由で無視できる。
reinterpret_castはコンパイル段階で完結する処理なので実行時コストはない。float *const sum_ptr = reinterpret_cast<float *>(sum):sumのアドレスがループ中で不変であることをコンパイラは知っているからsum_ptrはループに入る前に解決され、ループ中での更新はない。sum_ptr自体のアドレスをコード中で使用していないのでRAM上に配置する必要性が無く、おそらくレジスタに割り当てられる。x_ptr, y_ptrはトリップ毎に変わるが、それ以外はsum_ptrと同じ条件が成り立つのでおそらくレジスタに割り当てられる。
関数化
積和の計算の度に煩雑なコードを書くのは生産性・保守性において非効率極まりないので、強制的にインライン展開される関数として定義するのがよい。処理時間計測機能付きのサンプルコードを次に示す。
#include <chrono>
#include <complex>
#include <cstdlib>
#include <iostream>
#include <random>
#include <thread>
#include <type_traits>
/**
* @brief Add the product of two complex numbers "x1" and "x2" to "y"
* This function is expected to work about 20% faster than normal operation such as "y += x1*x2".
*
* @tparam T the number type of real and complex parts
* @param[in] x1 "x1"
* @param[in] x2 "x2"
* @param[inout] y "y"
*/
template <typename T>
inline static void __attribute__((always_inline)) addProd(const std::complex<T> &x1, const std::complex<T> &x2, std::complex<T> &y) {
const T *const x1_vec = reinterpret_cast<const T *>(&x1);
const T *const x2_vec = reinterpret_cast<const T *>(&x2);
T *const y_vec = reinterpret_cast<T *>(&y);
y_vec[0] += x1_vec[0]*x2_vec[0] - x1_vec[1]*x2_vec[1];
y_vec[1] += x1_vec[0]*x2_vec[1] + x1_vec[1]*x2_vec[0];
}
/**
* @brief Add the product of two non-complex numbers "x1" and "x2" to "y".
* Calling this function is completely equivalent to just writing "y += x1*x2".
*
* @tparam T the number type
* @param[in] x1 "x1"
* @param[in] x2 "x2"
* @param[inout] y "y"
*/
template <typename T>
inline static void __attribute__((always_inline)) addProd(const T &x1, const T &x2, T &y) {
y += x1*x2;
}
int main() {
/* Prepares sample points */
constexpr int N_sim = 1e8;
std::random_device seed_gen;
std::default_random_engine rand_engine(seed_gen());
std::uniform_real_distribution<float> dist(-1, 1);
std::vector<std::complex<float>> samplePoints(N_sim);
for (int i=0; i<N_sim+1; ++i) {
samplePoints[i].real(dist(rand_engine));
samplePoints[i].imag(dist(rand_engine));
}
std::complex<float> sum;
std::this_thread::sleep_for(std::chrono::milliseconds(1000)); // short rest
/* normal multiplication */
sum = 0;
const auto t0_normal_mul = std::chrono::system_clock::now();
for (int i=0; i<N_sim; ++i) {
sum += samplePoints[i]*samplePoints[i+1];
}
const auto t1_normal_mul = std::chrono::system_clock::now();
std::cout << "sum" << sum << std::endl; //avoid optimizing out `sum`
std::this_thread::sleep_for(std::chrono::milliseconds(1000)); // short rest
/* special method */
sum = 0;
const auto t0_specialMethod = std::chrono::system_clock::now();
for (int i=0; i<N_sim; ++i) {
addProd(samplePoints[i], samplePoints[i+1], sum);
}
const auto t1_specialMethod = std::chrono::system_clock::now();
std::cout << "sum" << sum << std::endl; // avoid optimizing out `sum`
/* Show result. */
const auto dt_normal_mul = t1_normal_mul - t0_normal_mul;
const auto dt_specialMethod = t1_specialMethod - t0_specialMethod;
const auto dt_ms_normal_mul = std::chrono::duration_cast<std::chrono::milliseconds>(dt_normal_mul).count();
const auto dt_ms_specialMethod = std::chrono::duration_cast<std::chrono::milliseconds>(dt_specialMethod).count();
std::cout << "dt_ms_normal_mul: " << dt_ms_normal_mul << " ms\n";
std::cout << "dt_ms_specialMethod: " << dt_ms_specialMethod << "ms\n";
// CPU: Core 2 Quad Q9650, RAM: DDR2 800MHz 8GiB
// N_sim = 1e8
// dt_normal_mul: 411 ms
// dt_specialMethod: 330 ms
}性能比較
次の表に示すのは、前掲のコードで処理時間を測定した結果である。繰り返し回数は結果を安定させるのに十分な値とし、10回以上実行した結果を平均した。最適化レベルは最高に設定している。
表の3行目は今使っている組込DSPであり、恐ろしいほどに効果が認められる。このプロセッサにはHWによるループ支援機能があるため、ループ内での一時オブジェクトの生成処理を取り去ったことで極めて高効率な機械語が生成されたものと考えられる。
| HW | OS | コンパイラ | 繰り返し回数 | 処理時間 [ms] | 処理時間削減量[%] |
| CPU: Core 2 Quad Q9650 RAM: DDR2 800MHz 8GiB | Windows 10 64bit | clang 11.0.1-2 | $10^8$ | normal method: 411 special method: 330 | 19.7 |
| CPU: Core-i5 1035G7 RAM: DDR4 8GiB | Windows 10 64bit | clang 11.0.1-2 | $10^8$ | normal method: 161 special method: 127 | 21.0 |
| CPU: TMS320C6748 RAM: DDR2 256MiB | TI-RTOS 6.83.0.18 | C6000 CGT v8.3.11 | $10^6$ | normal method: 763 special method: 191 | 75.0 |
