Expression templateとfmaについて

Fused Multiply Addという命令がある。これは、a * b + cという形の演算を1命令で処理するものだ。丸めが一度しか走らないので精度がよくなる。以前、理論値と数値計算を比較してテストするときにこいつが使われるかどうかでテストが通ったり落ちたりして困った思い出がある。

さて、必ずこれを使うようにするにはどうしたらよいか。簡単だ。C99ではmath.hfmafmaffmalが定義されており、この関数を使えばよい。C規格(N1256)の§F.9.10.1 "The fma functions"で、

fma(x,y,z) computes xy + z, correctly rounded once.

とのことなので、CPUがサポートしていなかったらソフトウェア側でエミュレーションするのだろう。実際FP_FAST_FMAというマクロがあり(§7.12 Mathematics <math.h>)、これは積と和を計算するよりもFMAの方が早い場合のみ定義される。FMAFとFMAL版もある。これはほぼ専用命令があれば定義されることと同義なので、これが定義されているかどうかで使い分けもできるだろう。もしCPUがFPUを持っておらず、浮動小数点数演算を全てソフトウェアエミュレーションしていれば(そしてエミュレーションの実装でFMAの方が速ければ)、FMA専用命令がなくとも定義されるかもしれない。いまどきそんなCPUあるか?

さて、C++C++11からstd::fmaがサポートされた。これはCのfmaと同等である。また、SIMDFMAでは以下のようなイントリンシックがある。

__m128 _mm_fmadd_ps(__m128 a, __m128 b, __m128 c);

FMA関係のイントリンシックはもちろんこの一つではなく、幅(__m128, __m256) かける 精度(_ps_pd) 足す 非SIMD_ss_sd) で10通りある。AVX512ではマスクがあるのでもっと色々ある。

さて、これを制御したいとする。ライブラリ側で、ユーザーには気づかれないうちに。

目的からして、ユーザーがどんな演算をするか何も仮定できない。とすると、どういう演算をするかという情報を取り出せる場所は一つしかない。ユーザーのコードだ。

先に書いておいたライブラリがユーザーのコードから情報を取り出す方法としては、まずExpression Template(ET)が思いつく。これは本来はBlitz++という線形代数ライブラリが必要のない要素の計算を避けるために開発した技法で、一種の遅延評価である。実装をかいつまんで書くと、まず以下のようなクラスを作って、演算子が計算結果の一次オブジェクトでなく、両辺への参照を持った特殊なクラスを返すようにする。

template<typename L, typename R>
struct added
{
    using value_type = decltype(
        std::declval<typename L::value_type>() +
        std::declval<typename R::value_type>()
    );
    value_type operator[](std::size_t i) const noexcept
    {
        return l[i] + r[i]};
    }

    L const& l;
    R const& r;
};

このクラスは2つの型(左辺のオブジェクトと右辺のオブジェクトになる)をとって、アクセスしたときに必要な値を計算して返すクラスだ。このクラスはアクセスされるまでは計算を行わない。よって、例えば巨大な行列の積を計算するが対角要素以外いらないとか、そういう時のために使われてきた。

さて、これを使うと、以下の式のdの型はどうなるだろうか。それぞれの変数はmatrix型だとする。

auto d = a * b + c;

まずa*bが行われるので、以下のようになる。

auto d = multiplied<matrix, matrix> + c;

その後、足し算が行われるので、以下のようになる。

auto d = added<multiplied<matrix, matrix>, matrix>;

というわけで、dの型はadded<multiplied<matrix, matrix>, matrix>だ。これは木構造になっている。というか、抽象構文木(AST)になっているのだ。

本来は実装の手間を省くために、std::plusなどを取る型パラメータを付け加えるためよりASTっぽくなるが、今回は例が複雑になるためやめた。

このように、ETを使うことでユーザーコードで出てくる変数の型をASTにすることができる。ここで、templateは部分特殊化、つまりパターンマッチが使えるということを思い出そう。すると、以下のような操作ができる。

template<typename L, typename M, typename R>
struct added<multiplied<L, M>, R>
{
    using value_type = decltype(
        std::declval<typename L::value_type>() *
        std::declval<typename M::value_type>() +
        std::declval<typename R::value_type>()
    );
    value_type operator[](std::size_t i) const noexcept
    {
        return std::fma(l[i], m[i], r[i]};
    }

    L const& l;
    M const& m;
    R const& r;
}:

これは、addedの左辺が何かしらのmultipliedだった場合への特殊化だ。そういうクラスが現れた時は、通常のaddedmultipliedではなくこちらで実体化される。

――という風にすれば、ユーザーが意識せずに書いてもfmaが勝手に生成されるという風にできるなあ、と思いつつも実装していない。面倒なのである。

だが、ライブラリ側でユーザーコードのASTを変換できる可能性というのは面白いのではないかなあと思う。コンパイラが最適化のためにいじっているのも何らかの中間表現なので、理論上は同等のものが実装できるはずだ。最近はETはコンパイラの自動SIMD化を阻害する(解析が難しくなる)ので下火という話も聞くが、こういう方向の発展もありなのではないか。

……それはコンパイラの最適化を阻害しつつ、自分でコンパイラがやってるのと同じ最適化をするってことだから、めっちゃくちゃ特殊な場合を除いて意味ないんじゃないの。

maveの中身について

前回SIMDベクトルライブラリを作った話をしたが、そこでやっていることについてかいつまんで書いておく。

まず、以下のようなクラスがある。

namespace mave
{
template<typename T, std::size_t R, std::size_t C>
class matrix;

template<typename T, std::size_t N>
using vector = matrix<T, N, 1>;
}

これについて、ごく普通の演算子の他に、以下のような演算子オーバーロードを行った。

namespace mave
{
template<typename T, std::size_t N>
std::tuple<vector<T, N>, vector<T, N>>
operator+(std::tuple<const vector<T, N>&, const vector<T, N>&> lhs,
          std::tuple<const vector<T, N>&, const vector<T, N>&> rhs)

template<typename T, std::size_t N>
std::tuple<vector<T, N>, vector<T, N>>
operator+=(std::tuple<      vector<T, N>&,       vector<T, N>&> lhs,
           std::tuple<const vector<T, N>&, const vector<T, N>&> rhs)
}

そう、あまり見る機会はないが、実はoperator+=もクラス外で定義できる。通常この演算子はメンバへの書き込みを伴うのでクラス内でメンバとして定義する(ことが多い)。しかし今回は、std::tupleにはメンバを追加できないのでfreeな関数として実装している。まあこの場合はコンストラクタで与えられるものが全てなので、メンバ関数として実装する必要がそもそもないのだが。

念頭にあるのは以下のようなユーザーコードである。

const mave::vector<float, 3> v1(/*...*/), v2(/*...*/);
const mave::vector<float, 3> w1(/*...*/), w2(/*...*/);
auto [u1, u2] = std::tie(v1, v2) + std::tie(w1, w2);

まず、std::tieは参照を格納したstd::tupleを作って返す。これは本来、タプルのパターンマッチをライブラリレベルでサポートするための機能だ。普通は以下のように使う。

std::tuple<int, double, std::string> f();

int main()
{
    int i;
    double d;
    std::string s;
    std::tie(i, d, s) = f();
}

これは、多値を返す関数の戻り値を分解しつつ受け取るコードだ。何が起きるかというと、まずstd::tieによってstd::tuple<int&, double&, std::string&>が作られ、そこにstd::tuple<int, double, std::string>が代入される。参照への代入なのでもともとのオブジェクトへの書き込みが行われ、idsのそれぞれが書き換わるという寸法だ。その後std::tieによって作成されたオブジェクト(のprvalue)の寿命は人知れず切れる。

ついでにこれはstd::ignoreを使うことで一部の値を無視することができる。

int i;
double d;
std::tie(i, d, std::ignore) = f();

このstd::ignoreは処理系定義なクラスのインスタンスと定義されている。代入演算子template化してどんな型でも代入できるようにしておき、実際には代入演算子は何もしないようにすれば実装できるだろう。こんな感じで。

struct ignore_t
{
    template<typename T>
    ignore_t& operator=(const T&) noexcept {return *this;}
};

これらの機能は便利だったのだが、std::tieだと初期化にコストのかかるオブジェクトを受け取りたい時に不必要なコストがかかってしまうことがある(先にデフォルトコンストラクタが走り、その後std::tieで受け取るときムーブコンストラクタが走る)。また、そもそもデフォルトコンストラクタを持っていないとstd::tieは使いにくい。適当な引数を与えてあとで代入することで使えないわけではないが、明らかに意味不明だろう。というわけで、結局コア言語側で構造化束縛が導入されることになった。

閑話休題std::tieは参照のタプルを作るので、参照のタプルを受け取って使うことにすればよい、というのは使い方を決めると自動的に決まる。要素数が2個の場合はstd::pairを使えるようにしようかと一瞬思ったが、微妙に複雑になってしまう(結局tupleもサポートしないといけないので転送するだけの関数が増えてしまう)のでやめた。

さて、これで大丈夫だろう、と思ったのだが、少し困ったことが起きた。先ほどの状態でこれは通ったのだが、

const mave::vector<float, 3> v1(/*...*/), v2(/*...*/);
const mave::vector<float, 3> w1(/*...*/), w2(/*...*/);
auto [u1, u2] = std::tie(v1, v2) + std::tie(w1, w2);

以下は通らなかった。

mave::vector<float, 3> v1(/*...*/), v2(/*...*/);
mave::vector<float, 3> w1(/*...*/), w2(/*...*/);
auto [u1, u2] = std::tie(v1, v2) + std::tie(w1, w2);

例えば、以下のようなtemplate関数があったとする。

template<typename T1, typename T2>
void f(std::tuple<T1 const&, T2 const&> tied)
{
    std::cout << std::get<0>(tied) << ", "
              << std::get<1>(tied) << std::endl;
}

これにstd::tuple<int&, double&>などを渡すと、const T1intがマッチしないのでこれは候補にならない。なので、以下は通らない。

int i = 42;
double pi = 3.14;
f(std::tie(i, pi)); // no matching function!

std::tie(T1, T2)std::tuple<T1&, T2&>を返すからだ。だが、もちろんT&const T&に変換できるので、templateをやめると通る。

void f(std::tuple<int const&, double const&> tied)
{
    std::cout << std::get<0>(tied) << ", "
              << std::get<1>(tied) << std::endl;
}

int main()
{
    int i = 42;
    double pi = 3.14;
    f(std::tie(i, pi)); // conversion from int& to int const&
    return 0;
}

これが通ってtemplateにすると通らないのはハマりどころな気がする。実際、書いてエラーを見るまでtemplateでも通ると思っていた。一度、呼び出す関数の解決順序を見直した方が良さそうだ。

さて、これをどうするかだ。多分templateをやめてmave::vector<float, 3><double,3>matrix<float, 3, 3>……とそれぞれに演算子オーバーロードを用意すれば良いのだが、その場合用意していないサイズのmatrixは、作れるが演算だけできないというよくわからないことになる。全てのサイズ、型についてオーバーロードを用意することはできない。そういう問題を回避する為にtemplateがあるのだ。では非templatetemplate演算子を共存させれば、と考えるわけだが、当然template版と非template版の両方がマッチするので、conflictしてオーバーロード解決が失敗する。

解決策として、こうした。

#define MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(FUNC_NAME, L_MODIFICATION, R_MODIFICATION)\
  //...
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, MAVE_EMPTY_ARGUMENT, MAVE_EMPTY_ARGUMENT)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, const&             , MAVE_EMPTY_ARGUMENT)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, &                  , MAVE_EMPTY_ARGUMENT)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, MAVE_EMPTY_ARGUMENT, &)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, const&             , &)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, &                  , &)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, MAVE_EMPTY_ARGUMENT, const&)
MAVE_GENERATE_FORWARDING_BINARY_FUNCTIONS_MATRIX_MATRIX_MATRIX(operator+, &                  , const&)

グワーッ!

渡されうる引数はconst&&(値渡し)の3通りがありえて、2引数取る場合はその組み合わせで増えるので大変なことになる。 これ、何かいい解決法がないだろうか。考えているのだが、浮かばなくて困っている。

SIMDベクトル・行列ライブラリを作った

TL;DR

高水準なAPISIMD命令を使って小さなベクトル・行列計算するためのライブラリを作った。

github.com

名前は安直だが MAtrix と VEctor からとった。

数年前にもExpression Templateを使って似たようなものを書いたことはあったが、これは的を絞った分ちゃんと狙いが定まった感がある。

Motivation

サーベイ不足は自覚しているので、以下は雑な感想である。正直に言うと、理由の半分くらいは「面白そうだったからやってみたかった」だ。

AVXをサポートしているCPUを持っているなら、floatが8個まで同時に加算できる。これは便利だが、ごく普通に3次元ベクトルに対して演算子オーバーロードをするとAVX命令を使うかどうかはコンパイラ任せになる。

ここで絶対に8個同時に計算したい、という場合はあるだろう。というわけで、以下のような演算子をサポートした。

mave::vector<float, 3> v1(/*...*/), v2(/*...*/);
mave::vector<float, 3> w1(/*...*/), w2(/*...*/);
auto [u1, u2] = std::tie(v1, v2) + std::tie(w1, w2);

C++17の構造化束縛が使えない場合は、std::tieを使うと良い。

maveは3次元ベクトルのためのライブラリなので、長さを求めるなどの計算も同様に出来て、

// l1 = mave::length(v1 - w1);
auto [l1, l2] = mave::length(std::tie(v1, v2) - std::tie(w1, w2));

と書ける。maveでは、4つまでのtieをサポートしている。

これらの演算子は、対応するSIMD命令が存在しなかった場合はできる量に分けて実行する。なので、常に最大個数で書いておいても大丈夫だ(レジスタの量が足りなかったり、命令の順番が悪くてパイプラインの効率が悪くなる可能性はあるし、端数の処理で減速する可能性もある)。

よくある線形代数ライブラリの目標は大きな行列やベクトルの複雑な計算を高速にすることなので、こんなケースのサポートは後回しだろう(と思う)。 また、maveは小さなベクトルの扱いに特化しているので、簡単なlengthなどの関数を用意している。このあたりが汎用のSIMDライブラリと異なる点だ。まあ、汎用のゼロオーバーヘッドSIMDライブラリを作ってからその上にこれ実装すれば良いのでは、というのはあるが、まあそれは置いておいて。

仕組みの概要

中で何をしているかに関しては、意外と大変なことがいくつかあったので回を改めてちゃんと書きたいが、大雑把に言うと要するにtieによって作られるタプルのそれぞれについて演算子オーバーロードをしているわけだ。 クラスの外でもoperator+operator+=は実装できるので、特定の中身を持ったstd::tupleについて演算子オーバーロードできる。ここに注目して、上記のようなインターフェースを作っている。

boost::compressed_pairの使いドコロ

使いドコロがわかりにくいと噂のboost::compressed_pairというものがある。今日はこれの具体的な使いドコロを紹介してみたい。

ところでなぜ私の使っているIMEは「つかいどころ」を「使いドコロ」と変換するのだろう。まあそんなに違和感もないのでこのまま行くことにする。

背景

さて、そもそも存在を知られていない気がするので、boost::compressed_pairの説明を少ししていこう。これは、std::pairと基本的には同じなのだが、どちらかの型引数が空クラスだった場合にその領域を消し去るというものだ。

そもそも、通常はstd::pair<int, void>などは作ることができない。以下はコンパイルエラーになる。

#include <utility>
#include <iostream>

void f(){return;}

int main()
{
    std::pair<int, void> p(42, f());
    std::cout << p.first << std::endl;
    return 0;
}

[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ

ではこういう場合どうしようか。例えば、空のクラスを作り、それをいれてみる。

#include <utility>
#include <iostream>

struct empty_t{};

empty_t f(){return empty_t{};}

int main()
{
    std::pair<int, empty_t> p(42, f());
    std::cout << p.first << std::endl;
    return 0;
}

ところで、このペアのサイズはどれくらいになるだろう。よくある環境なら、これは8バイトになる。だが、intは4バイトだ。 何故8バイトも取るのか? その理由はアライメント要求にある。intは4バイトアライメントを要求する(大抵の場合は)。これは何を意味するかというと、intは4の倍数アドレスにしか置けないということだ。 なので、pairを例えば連続したメモリ領域に置こうとした時(std::vectorにいれるなど)、intempty_tintと来たこの2つのintは、両方共4の倍数に置かれなければならない。 もしempty_tがゼロでないサイズのオブジェクトなら、間隔は4バイト以上にならざるを得ない。なので、4バイトの倍数になるようにパディングが入ることになる。そして実際、empty_tは1バイト要求する。

なぜempty_tが1バイト要求するのか? C++ではアドレスがオブジェクトの一意性に関わっているので、空のクラスでも一意的なアドレスを取れるようにしないといけない。例えば空のクラスを構造体のメンバとして使った場合、もしサイズがゼロだったら、空クラスへのポインタはその構造体中の次のメンバ変数の領域を指してしまうだろう。その要素が最後に配置されたとしても、メモリ空間中で隣接する別のオブジェクト(配列の次の要素など)を指してしまうかも知れない。そういうことを避けるため、最低限の1バイト(多分sizeof(T) == 1になるようにしていると思うが、これは大抵の場合1バイトだ)を取ることになる。中身は不定だと思う。

そういうわけで、このペアは本質的にint一個分のサイズで済むはずなのだが、実際には2個分要求することになってしまう。そして中身の半分はガラクタだ。いくらメモリが潤沢な現代とはいえ、CPUキャッシュまで潤沢なわけではない。むしろメモリ上でのデータの局所性が重要になってきている現代において、ガラクタによって必要なデータが遠くに行ってしまうのは辛いだろう。

というわけで、何とかしてこういう場合にサイズをint一個分で済ませるというのがboost::compressed_pairの目的である。

使いドコロ

そもそも空のクラスとのペアって何だよ、ペアじゃなくていいじゃないか、というのは妥当なツッコミである。一体全体どういう時にそんな問題が発生するのか?

この使いドコロは、メタプログラミングにある。例えば、重い処理をするにあたって先に計算を済ませておける部分があるとき、済ませておきたいというのは当然の欲求だ。そこで、値のインデックスと先に計算しておいた値のペアを置いておきたいとする。

template<typename Calculation>
struct cache
{
    std::vector<std::pair<std::size_t, typename Calculation::parameter_type>> values_;
};

だが、先に計算しておくべきパラメータが存在しないCalculationというのがあり得る。この場合、parameter_typevoidにすると上記の構造体はコンパイルできないし、empty_typeにするとメモリを無意味に消費してしまう。そいう場合のためにboost::compressed_pairが存在する。

少し話が抽象的すぎてよくわからないので、具体例を出してみる。読者は分子動力学シミュレーションをご存知だろうか。

分子動力学法 - Wikipedia

なんだかゴチャゴチャ書いているが、要するに原子間にかかる力(分子力場)をうまく設定して、ニュートン運動方程式を数値積分して原子が動くのを見て楽しむことだ。よく材料力学や生物物理学などの分野で活用されている。例えば、イオン同士の間に静電相互作用のポテンシャルを入れるとか、化学結合についてバネ(調和振動子)ポテンシャルを入れるとか、そういうことをする。

さて、この中でも最も基本的だと思われる、Lennard-Jonesポテンシャルを例に上げることにしよう。これは単純な式ながら分子間力のかなり良い近似になっているので、殆どの分子力場に何らかの形で含まれている。分子間力は2つの要素からなり、瞬間的な電子分布のゆらぎによる双極子間相互作用、つまりロンドン分散力による引力と、電子同士が近接した時のパウリの排他率による斥力である。いろいろおもしろい話はあるが、そういう話をしにきたわけではない。

さて、分子シミュレーションの計算量は多い。Lennard-Jonesのような分子間力は全ての粒子の間に作用するので、計算には全てのあり得る粒子のペアを考える必要があり、これにはO(n2)の時間を要する。これで粒子数が1万とかだと(この数はタンパク質のシミュレーションでは標準的な数だ)、愚直に計算していると統計力学的に意味のある話ができるくらいの時間シミュレーションを実行するには信じられないくらいの時間がかかる。

で、効率化したいというのは当然の欲求ということになる。一般に、最も効率のよい最適化は、計算しないことだ。分子間力にはカットオフという概念がある。多くのポテンシャルは無限遠点で値がゼロになるようにされていて(これは強い相互作用によるクオークの閉じ込めみたいな例外的な場合を除いて妥当な話だろう)、そのため、ある有限の距離でも無視できる程度の弱さまで減衰する。そのため、ある程度以上遠い粒子の間の力は計算しなくていい。これがカットオフだ。

これだけでもある程度速くはなるが、どうせ計算しないペアを見てスキップするというのは少しアホらしい。このため、先に「相互作用しそうなペア」をリストアップしておいて、そのペアだけについて計算するという方法が取られている。もちろん粒子はシミュレーション中に動くので、このリストを時々アップデートする必要があるわけだが、そこらへんのトリックは各自調べてほしい。簡単に説明しておくと、先にカットオフよりも長い距離でペアを探索しておいて安全マージンを作っておき(このためリストの一部は無駄になるが、まあそれはそんな数ではないので)、粒子が動いた分だけマージンを減らしていくという感じだ。

さて、話が長くなってきたが、もう少し説明を追加しないといけない。このLennard-Jonesポテンシャルは2つのパラメータを取る。粒子の大きさと、力の強さだ。もちろん、これらのパラメータは原子の種類によって変わる。特に、これは原子のペアに依存するので、原子の種類が増えると組み合わせ的に増大する。これはやっていられないので、各原子についてのパラメータからペアについてのパラメータを決める方法が考えられている。ローレンツ・ベルテロ則を使えば、各原子のパラメータだけを覚えておいて、ペアについてのパラメータは相加平均と相乗平均でオンデマンドに計算することができる。

やっと本題だが、相乗平均はsqrtを必要とする。だがsqrtは重い処理だとして知られている。そのため、力の計算のときに何度もしたくはない。ところで、我々は相互作用しそうなペアのリストを作って使いまわしているのだった。なら、このリストについでに相乗平均の値を持たせてしまえばよい!

というわけで、以下のようなクラスが発生するだろう。

template<typename Potential>
struct NeighborList
{
    // list of { pair of {
    //    neighboring particle index, pre-calculated parameter
    // } }
    std::vector<std::pair<
        std::size_t, typename Potential::value_type
        >> pairs;

    // list of { range of {neighbor of i-th particle} }
    std::vector<std::pair<std::size_t, std::size_t>> neighbors;
};

ああ、ようやく戻ってきた。だが、このPotentialの中には、例えばパラメータが一様なのでそもそも相乗平均なんかとらなくていいとか、事前計算が掛け算一個とかで格納するメリットがほぼ無いとか、そういうものが後々出てくるだろう。そうすると、これはまさしく上記の問題に遭遇することになる。こういう場合にcompressed_pairを使えば、何も気にすることなく全ての場合に対応できる。そうでなければ、意味不明なメタプログラミングか(constexpr ifを使えば簡単になるかもしれない)、メモリの無駄遣いに耐えるしか選択肢がなくなる。

長くなりすぎてしまった。少し珍しい場合のことを考えるためには、前提を共有する必要があって説明が不必要に長くなってしまう。抽象的な例を挙げているだけの方がマシだったかもしれない。

仕組み

ところで、ではcompressed_pairは何をどうやって空クラスのサイズをゼロにしているのだろうか。これを理解するには、empty base optimization (EBO)というものを説明せねばなるまい。基底クラスには、その派生クラスとポインタの値が違わなければならないという保証はないので、基底クラスが空の場合、本当にサイズをゼロにすることができる。このとき、基底クラスへのポインタは派生クラスのポインタと一致する。さらに特定の場合、最初のメンバ変数へのポインタとも一致するが、それはさておき。

EBOをうまく使っている例はちょいちょい見受けられる。例えば、short string optimization(SSO)だ。これは文字列が十分短い時は動的確保をせずにクラス内に持ってしまうというもので、できるだけ小さいクラスにできるだけ多くのchar buf[N]を詰め込むために、いろいろなテクニックを駆使している。その一つがEBOで、アロケータはたいてい状態を持たない(単にnew/deleteするだけ)ので、空のクラスだ。同じサイズのクラスに1文字でも長い文字列を格納したいので、これに費やす1バイトをcharに割り当てたい。なので、他のメンバを一度構造体のメンバとして持ち、その構造体がアロケータを継承することでEBOを有効にしている。SSOは以下の記事で詳しく解説されている。

qiita.com

さて、もうお分かりの通り、compressed_pairもEBOを活用している。ペアの片方が空だった場合、それを継承してしまうのだ。するとEBOによって空のクラスのサイズが消え去る。使う側は何も気にしなくてよい、というわけだ。

まとめ

というわけでboost::compressed_pairは、template引数に応じて必要だったり不要だったりする値があるとき、メタプログラミングなしでメモリの使い方を最適化できる便利なクラスだ。 これを期に検討されてはいかがだろうか。

追記:リテラル数値について

ところで、昨日の記事で、以下のコードはコンパイルエラーになって面倒だという話をした。

template<typename T>
T f(T x, T y)
{
    return x + y;
}

auto z = f(1.0f, 2.0); // f(float, double)

一応追記しておくと、これはtemplateだから生じる問題とかではなくて、オーバーロードでも普通に生じる。

double f(double x, double y)
{
    return x + y;
}
float  f(float  x, float  y)
{
    return x + y;
}

auto z = f(1.0f, 2.0); // f(float, double)

もちろん、理由はambiguousだからだ。

prog.cc: In function 'int main()':
prog.cc:14:25: error: call of overloaded 'f(float, double)' is ambiguous
     auto z = f(1.0f, 2.0);
                         ^
prog.cc:3:8: note: candidate: 'double f(double, double)'
 double f(double x, double y)
        ^
prog.cc:7:7: note: candidate: 'float f(float, float)'
 float f(float x, float y)
       ^

一般に、浮動小数点数リテラル値を書くときは色々と面倒が生じるということになる。

一度明示的に定義した型で受けてしまえば、この手の問題は発生しない。

float x = 1.0 / w;
float y = 2.0 * w;
auto z = f(x, y);

冗長ではあるが。

リテラル数値と精度について

背景

C++では1.0とか書くとdoubleとして解釈され、floatにするには1.0fなどと書かねばならないことはご承知のとおりだ。また、doublefloatの演算の結果は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 * xdoublefloatの積なので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で実体化しつつ、a2.0bstd::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));
                                                  ^

[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ

エラーの内容は同じで、型推論の結果doublefloatかわからなかったよ、という話だ。

速度の話

と、さてここまでで大体めんどくさそうだなあということはわかっていただけたと思う。

面倒くさいだけならまだ良いのだが。以下のリンク先を見て欲しい。

godbolt.org

このリンク先は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.04.0f4と変えていって、doublefloatで実体化した。ここに載せているのは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"の略である。要するにfloatdoubleに変換するもので(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の方が速かったりする。

その辺りは詳しくは以下の記事などを参照していただきたい。

qiita.com

さて、対比のために単精度浮動小数点数リテラルを使った場合のアセンブリを見てみよう。

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

変換命令はないし、subssmulssが使われている。やはり、倍精度リテラルを使ったことで問題が生じたのだと考えるのが自然だろう。

まとめ

中に入れる型を決めないままコンテナだけを決める

背景

これはある実験機器から出てくる特殊なフォーマットの実験データを読むために作ったライブラリのために考えたやり方である。

そのデータは画像データとヘッダ情報から出来ており、基本的にはそれらのペアを返すことになる。だが数百KB〜数百MB程度と大きさが色々なので、メモリアロケーションの戦略を色々と変えたいことがあるだろうと考えた。

というわけで、read関数にオプションtemplate引数としてコンテナ+アロケータの種類をカスタマイズできるものを渡したい。だが、画像データとして生データを取得するか、ヘッダ情報から単位換算などを済ませたデータを倍精度・単精度浮動小数点数として読み込むかは変更できるようにもしたかった。すると、データクラスはそもそもデータ型が不定ということになる。その状況でどのようにしたらコンテナの種類だけ変えられるだろうか?

template<typename T, typename ???>
class ImageData
{
    using value_type = T;
    using container_type = ???;
    // 第二template引数で`container_type`を
    // `std::vector<T>`や`std::deque<T>`に変えられるようにしたいが……
};

解決策

みなさんはstd::allocatorのメンバ構造体rebindをご存知だろうか。細かいことを言うとC++17でdeprecatedなので知らなくても良いのだが、その場合はstd::allocator_traits::rebind_allocを見ていただきたい。これは異なる型を同じ戦略でアロケートするクラスを取り出すための型である。

例えば、連結リストは実際にはノードをアロケートする必要があるが、std::listではユーザーにノードの実装は見せていない。見せる必要もない。だが、アロケートの戦略は変えられるようになっている。すなわち、第二引数にAllocatorを取る。

namespace std
{
template<typename T, typename Alloc = std::allocator<T>>
class list;
}

アロケータには本当はstd::allocator<_Node<T>>などを渡したいわけだが、このノード型がどうなっているかはユーザーには不明である。なのでどうあれユーザーはstd::allocator<T>を渡すことになる。すると、std::listはどうにかして、ユーザー指定のアロケート戦略を尊重しつつ、_Node型をアロケートしなければならない。このために、allocatorrebindメンバを持っており、異なる型を同じ戦略でアロケートするクラスを取り出せるようになっているのだ。

何言ってんだお前という人のためにサンプルを載せておく。

template<typename T>
struct allocator
{
    template<typename U>
    struct rebind {typedef allocator<U> other;};
};

allocatorの中に構造体が定義されており、その構造体がtemplateになっている。これによって、typedef allocator<T>::template rebind<U>::other allocator_UのようにしてUのためのアロケータを取り出せるのだ。

この発想が転用できる。

まず、以下のようなクラスを用意しておく。

struct vec
{
    template<typename T>
    struct rebind
    {
        typedef std::vector<T, std::allocator<T>> other;
    };
};

struct deq
{
    template<typename T>
    struct rebind
    {
        typedef std::deque<T, std::allocator<T>> other;
    };
};

これが渡されてくる前提で、最初のクラスで以下のようにしておく。

template<typename T, typename tagT>
class ImageData
{
    using value_type = T;
    using container_type =
        typename tagT::template rebind<value_type>::other:
};

ユーザーが呼ぶための関数を用意する。

template<typename T, typename tagT = vec>
ImageData<T, tagT> read(std::string const& file_name);

template<typename T, typename tagT>
ImageData<T, tagT> read(std::string const& file_name, tagT);

カスタムする気のないユーザーは、普通にreadを呼べば良い。カスタムする気のあるユーザーは、テンプレート引数に渡すか、第二引数としてそのクラスを渡すとよいことになる。

const auto data = read<double>("sample.dat");

const auto data1 = read<double, deq>("sample.dat");
const auto data2 = read<double>("sample.dat", deq{});

さらに、先にユーザーにこのクラスが持っておくべきものを伝えておけば、ユーザーが作ったアロケータやコンテナすら射程に入る。

struct user_alloc
{
    template<typename T>
    struct rebind
    {
        typedef std::vector<T, user::allocator<T>> other;
    };
};

struct user_container
{
    template<typename T>
    struct rebind {typedef user::container<T> other;};
};

ユーザー指定のコンテナのインターフェースがstd::vectorなどと同じなら、ライブラリの実装は買える必要はない。

違う場合も、メンバ関数呼び出しをフックできるようにしておけば最小限の変更で対応できる。

例えばsize関数呼び出しには、以下のような関数を用意し、ライブラリ内部ではそれだけを使っておけば、自分でコンテナを作るレベルのユーザーはそれをフックして対応してくれるだろう。

// size() メンバ関数があればそれを呼ぶ
template<typename T, typename std::enable_if<
    has_mem_func_size<T>::value, std::nullptr_t>::type = nullptr>
std::size_t size(const T& container);

// sizeメンバがないコンテナを使うユーザーは、オーバーロードを追加する
template<typename T>
std::size_t size(const user::container<T>& c)
{
    return c.get_size();
}

他の解決策

ところで、ImageDataクラスはtemplate template引数をとるような実装もできる。

template<typename T, template<typename> class containerT>
class ImageData
{
    using value_type = T;
    using container_type = containerT<T>;
};

template<typename T>
using deq = std::deque<T>:

const auto data = read<double, deq>("sample");

この実装はシンプルだし、template引数だけで完結するなら何も問題はない。どうせstd::vectorより凄いコンテナを作ってくる人間は限界突破しているはずなので、この程度にはひるまないだろう。

唯一なにかあるとしたら、

const auto data = read<double>("sample.dat", deq{});

のような呼び出しが出来ないことだろうか。ディスパッチ用の何かは今回エイリアスなので、実体化出来ない。まあしかし、それだけだ。