基本ベクトルの場合の最適化(3次元幾何)

数カ月触れていなかった自分が書いたコードを見た所、驚くべきものを発見した。

3次元空間内をある方向に向けてある距離だけ動いたオブジェクトが衝突するかどうかを判定し、また衝突するときはどれだけ動いたら衝突するか計算するコードがあった。 このようなロジックは、通常ベクトルを用いて方程式を立てて解く。その結果を実装するわけだが、この場合オブジェクトの動きを表すベクトルが以下のような特殊なものであった場合、計算がかなり簡略化できる。 - 単位ベクトル(単位ベクトルに変換する処理が含まれることが多いため) - 軸に平行なベクトル(複数の軸の成分を考えなくてよいため) よって、この性質を併せ持つ軸に沿った単位ベクトル、つまり基本ベクトルの場合、計算が非常に簡単になるわけだ。

問題は、実行中に来たベクトルが基本ベクトルかどうかを判定するのは若干手間であることだ。 通常、数値計算の結果(1, 0, 0)になる、というような値は、大抵数値誤差などのせいで(0.999999999, 7e-15, 5e-16)のような値になる。 もちろん適当なtoleranceを用いて判定することは可能だが、軸に平行な単位ベクトルが来ることがあまり期待できない場合、これは単にコストになってしまう。 さらに、ゲームやシミュレーションなどを考えた時、軸に沿ってちょうど距離1だけ移動する、というようなイベントはそうそう起きるものではない。

ではそれをする意味はどこにあるのか。それは、ユーザーからの入力がある場合である。 ユーザーからの入力は、まだ偶然生じるよりは基本ベクトルになる可能性がある。特に軸に平行なベクトルになる可能性はある(空間中に突然オブジェクトを作り、そのまま下に降ろしていく、など)。 そう思ったのかは知らないが、数カ月前の自分はそれを想定していたらしい。コード中に、3次元ベクトルクラスの他に基底ベクトルが存在したのだ。

基本ベクトルクラスは通常のインターフェースを持つが変更は許さず、積や和を作ると通常のベクトルクラスを作成して返す、というようになっている。 長さを計算するルーチンのみならずドット積やクロス積もオーバーロードされて高速化されている。 衝突判定のコードはというと、is_pointなどのメタ関数によるSFINAEが飛び交っており、非常に面白くなっていた。

このようにすると、型レベルで分岐が行われるので(コンパイル時に上手く決められるようにすると)基本ベクトルが来ることを保証できる。 なので、基本ベクトルを前提にした最適化が許されるというわけだ。

まあ、手間に見あうとはあまり思えないので、今から消すのだが。単にcollision_zのような関数を定義すれば済んでしまうので、ホンの少しだけ不格好なことを除けば何らデメリットはない。

コンパイル時ルックアップテーブル生成について

(今更感)

目的

非負整数値を取る関数があるとする。そしてそれを実行中に非常に多くの回数呼ぶ(数百兆回とか)、ということは、科学技術計算ではよく見られる光景である。 整数値を取るので、恐らく何度もピッタリ同じ値を計算することになるだろう。 すると、何度も呼びそうな範囲は先に表にしておいて、それを見に行くことで計算時間を短縮する、という発想は自然なものに思える。 配列アクセスが計算より早い場合はこのテクニックは非常に有用だ。ルックアップテーブルという名前まで付いている。

これが実数値を取る関数だと全く同じ浮動小数点数が来る確率は低いので単純なルックアップテーブル作成は難しくなる。 その場合、メモリが足りるなら十分小さい幅ごとの代表値を表にしておくとか、計算が本当に重くて簡単な補間の方が早い場合は、一定値おきの値を使って引数での値を補間するとかいった方法を用いたりする。

さて、C++ではコンパイル時に計算ができる。つまり、コンパイル時にルックアップテーブルが作れるのだ。 ルックアップテーブルを作るとき、普通はビルドスクリプトの中で適当なスクリプト言語から何かライブラリを呼び出してヘッダファイルを生成するのだが、面白いのでコンパイル時に作ってみよう。

例:階乗関数

インターフェース

例として、単純なので階乗関数を考えてみたい。N=100くらいまでのルックアップテーブルをコンパイル時に作ることを考える。 普通に実装すると以下のようになるfactorialだが、

template<typename Real>
constexpr Real factorial_impl(std::uint32_t i) noexcept
{
    return i==0 ? 1 : static_cast<Real>(i) * factorial_impl(i-1);
}

ルックアップテーブルが完成すれば、factorialは以下のような実装になるだろう。

template<typename Real>
constexpr Real factorial(std::size_t i) noexcept
{
    return factorial_table<Real>::get(i);
}

では中身について考えていきたい。factorial_tableは以下のようになっていると考えられる。

template<typename Real, std::size_t N>
struct factorial_table
{
    constexpr static std::size_t size = N;
    constexpr static std::array<Real, N> value = /* ... ??? ... */;

    constexpr static Real get(std::size_t i) noexcept
    {
        //C++11ではstd::array::operator[]() constはconstexprではない。
        //C++14からstd::array::operator[]() constがconstexprに、
        //C++17からstd::arrayのメンバの殆どがconstexprになる。
        return factorial_table<Real, N>::value[i];
    }
};

template<typename Real, std::size_t N>
constexpr std::size_t factorial_table<Real, N>::size;
template<typename Real, std::size_t N>
constexpr Real factorial_table<Real, N>::value[N];

/* ...???... */の部分をどうするかだ。

配列生成

こういうときのために、C++14にはstd::index_sequenceなるものが存在している。これは整数シーケンスを表現する型、integer_sequencesize_tに対する特殊化である。 cpprefjpによると以下のようなものだ。

integer_sequence - cpprefjp C++日本語リファレンス

namespace std {
  template <class T, T... I>
  struct integer_sequence {
    using value_type = T;
    static constexpr size_t size() noexcept { return sizeof...(I); }
  };
}

これを使うと、variadic templateで渡した整数値を関数内で展開できる。 variadic template argumentの展開は、関数を噛ましつつ行うことができるので、以下のような関数が書ける。

template<typename Real, typename Integer, Integer ... vals>
std::array<Real, sizeof...(vals)>
generate_array(std::integer_sequence<Integer, vals...>)
{
    return {{factorial_impl(vals)...}};
}

パラメータパックの展開により、この関数を以下のようにして呼ぶと

constexpr auto a = generate_array<double>(
    std::integer_sequence<std::uint32_t, 0,1,2,3,4,5>{});

このa{factorial(0), factorial(1), factorial(2), factorial(3), factorial(4), factorial(5)}になるのだ。 これを使うことで、コンパイル時に好きな整数シーケンスを好きな関数で変換した静的配列を生成できる。 しかしこのままだと手でinteger_sequence<std::uint32_t, 0, 1, 2, 3, 4, 5, 6, ..., 100>みたいなのを書く必要があり、無理になってしまう。 もちろんC++がそんなことを許すはずはない。std::make_integer_sequence<Int, N>なるものがあり、これは<Int, 0, 1, ... , N>を一発で作ってくれる。 これにより、

constexpr auto a = generate_array<double>(
    std::make_integer_sequence<std::uint32_t, 100>{});

と書くことができるのだ。

make_integer_sequence - cpprefjp C++日本語リファレンス

C++11縛り

ただ、C++11だとこのinteger_sequenceは存在しない。縛りプレイのために、一応、これをどのように作るか少し考えてみよう。 少々面倒なことが起きるので、まずsize_t決め打ちでindex_sequenceを作る。integer_sequenceは、その後中身をstatic_castすることで作れる。 負の長さを持つ列、というのはないので、深く気にしなくていいのは楽だ。

template<std::size_t ... vs>
struct index_sequence{};

template引数だけが必要なので、中身は別にいらない(sizeとvalue_typeくらいは定義すべきだが、煩雑になるのでやめた)。 0, ... Nみたいなのをつくりたいので、再帰{0, 1, ..., N-1} + {N} -> {0, 1, ..., N-1, N}のようにしていくことを考えた。 そのためには、まず2つのindex_sequenceを結合する必要がある。

template<typename T1, typename T2>
struct index_sequence_concatenator;

template<std::size_t ... v1, std;;size_t ... v2>
struct index_sequence_concatenator<
    index_sequence<v1...>, index_sequence<v2...>>
{
    typedef index_sequence<v1 ..., v2 ...> type;
};

これを使って{0, ..., N}を生成しよう。

template<std::size_t N>
struct index_sequence_generator
{
    typedef typename index_sequence_concatenator<
            typename index_sequence_generator<N-1>::type,
            index_sequence<N>
        >::type type
};

template<>
struct index_sequence_generator<0>
{
    typedef index_sequence<0> type
};

というわけでできた。使いやすいようにエイリアスを定義しよう。 長さがNで0始まりなので{0, 1, ..., N-1}になることをすっかり忘れていたのでエイリアスでこっそり修正しておこう。

template<std::size_t N>
using make_index_sequence =
    typename index_sequence_generator<N-1>::type;

というわけで、C++11でもindex_sequenceが作れて、コンパイル時ルックアップテーブルができた。

記事のライセンスについて

ローカルに翻訳記事が溜まってきている。ブログやチュートリアルの翻訳が多いが、専門書の要約や一章まるまるの翻訳などもある。英語のまま読む方が一読する場合の効率はよいが、読んでいる途中でこんがらがる可能性があるときや見返す可能性が高い場合に翻訳メモを残したりする。そのメモを補完したのち全体を手直しすればおおよそ翻訳記事が出来上がる。

自分で使うだけでも十分役に立っているが、結構な時間を注ぎ込んでいるのでローカルに置いて埋もれてしまうのは勿体無い気もする。公開したいという気持ちはあるが、一部は明らかにライセンス違反になるだろう。ブログもそうだが、専門書はそもそもの著作権が切れていなかったり、翻訳権が日本の出版社によって購入されていたりする。そういうものの翻訳を勝手に、例えばこのブログで公開すれば、おそらく何らかの法に触れるだろう。ブログはそもそもライセンスが明記されていないことが多い。

CC BYやCC BY-SAなどでライセンスされていれば非常にわかりやすいので公開できるのだが、何も書かれていないとどうしようもない。プライベートなストレージに記事が溜まっていくばかりである。

とか言いつつこのブログもライセンスについて何も言っていない気がする。どうせ自分のブログを引用したり再頒布したい奴などおるまいと思っているせいだ。しかしこういう文句を垂れておきながら自分はライセンスについて何も言わないというのはあまり筋が良くなかろう。決めるとしたらCC BYにしたいところだが、はてな全体にライセンス条項があったりするのだろうか(はてなのサービスを使うためにはCC BY-NC-SAでなければならない、など)。あとで調べてみなければなるまい。

エラー時はメッセージを返すが、成功しても返すものがない関数

今、Cライブラリをラップしているのだが、副作用を起こす関数の戻り値をどうするか考えている。

前提:全体のエラーハンドリング

前提として、エラーハンドリングにはここまで例外送出でなくboost::optionalexpectedを用いている。 基本的に、value_or()でなく無理やりunwrap()的なことをしたり、std::stringstd::vectorなどが内部で例外をぶん投げたりしない限り例外は送出されないようにしている。

ところで、例外はどこからどのタイミングで投げられるかとか、どこでキャッチするのかを考えるのが非常にしんどい。 コードが長くなるにつれて飛んでくるかもしれないものと考慮しないといけない条件がどんどん増えていくので(リソース確保途中で飛んできたら解放を保証できるか?)、最近はもう諦め気味で、例外が飛んできたらもうstd::terminate()で良いかな……という状態になった。 それならoptionalを使って毎回受け取る度にチェックする方がまだマシである。

エラーメッセージなしに失敗する可能性のある関数は、以下のようなシグネチャになる。

optional<result_type> may_fail(args...);

ラップする対象のライブラリでエラーメッセージが取得できる場合は、以下のようにしている。

expected<result_type, boost::string_view> may_fail(args...);

Cライブラリなので、get_error()がライブラリ内部に取られたC文字列へのポインタを返す。これをstring_viewにくるんで返している。 std::stringは(Short String Optimizationが働かない場合)動的メモリ確保を行うが、boost::string_viewはそれを行わないため速度面で優位であると予想できる。

ただし何の代償もないわけではなく、このやり方だとライブラリ側がエラー文字列をfreeしたり上書きするとstring_viewが不正になってしまう。 まあ次のAPI呼び出しよりも先にエラーチェックをしろという話だが、これに関しては、expected<T, E>Tから変換可能な型UEから変換可能な型Fを持ったexpected<U, F>に変換可能にしておけば、boost::string_viewstd::stringへ変換することでエラーメッセージを保持できる。

expected<result_type, std::string> result = may_fail(args...);

とりあえずはこれでいいかな、という感じでここまで来た。

本題:エラーが発生するかもしれないが、成功した場合返すものがない関数

さて、副作用がある関数がある。例えば画面に四角形を描画するとしよう。APIは以下のようなものだろう。

int draw_rectangle(window* dest, const rect* rectangle);

返却値が0なら成功、-1なら失敗、というように決まっている。失敗したらget_error()を使ってconst char*を受け取る。 この関数をこれまでと一貫性を保ってラップするとしたらどうすればいいか。

状況を整理しよう。この関数は既知のWindowに既知の四角形を描画することが目的だ。 なので、成功した場合返すものが存在しない。ただし、失敗した場合はその理由を書き込んだ文字列を取得できる。

もちろんこれはWindowの状態を変えるので状態変更後のWindow構造体を返しても良いわけだが、このCライブラリはそのようにできてはいないし、毎度のコピーは速度を下げる。 純粋関数型言語コンパイラならそういう状況を考慮して作られていそうで、最適化が十分仕事をするかもしれないが、C++コンパイラにそのコピーを消し去るような最適化を要求できるかというと微妙な気持ちになる。 まあムーブをし続ければまだましなのかな。

なので、ここは「成功した場合返すものはないが、失敗した場合はエラーメッセージを返したい」としよう。 すぐさま思いつくのは、「無いかも知れないもの」を包むoptionalだろう。だがこれには、ここで使うには致命的な問題がある。

使ってみるとしよう。上の例を継続して使う。

optional<boost::string_view> draw_rectangle(window& win, const rect& rectangle);

// ...
auto result = draw_rectangle(win, box);
if(result) {
    /* 1: error? */
} else {
    /* 2: error? */
}

どちらがエラーっぽいだろうか。当然、else節だろう。普通、boolへの変換は成功した場合にtrueだ。他の場合はどれもそうなっている。 だが、optionalはあくまでも「値がある場合にtrue」なので、この使い方だと「エラーメッセージがある場合(失敗した場合)にtrue」になってしまう。

例えばウィンドウを作成して四角形を描画する一連の流れを書くと、以下のようになるわけだ。

auto win = make_window({640, 480});
if(win) {
    auto result = draw_rectangle(*win, box);
    if(result) {
        std::cerr << "an error occured!: " << *result << std::endl;
    } else {
        // successfully drawn. go next step.
    }
    // ...
} else {
    std::cerr << "an error occured!: " << win.unwrap_error();
}

このコードのどこに一貫性があると言うのか。一瞬でifelseの関係が逆転するこの設計は、明らかにミスを誘発する。

返すものがないならエラーコードや真偽値を返せばいいのでは? とも思わなくはないが、ユーザーにget_error()呼び出しをさせるのは避けたい。 ライブラリが提供するエラーメッセージは上書きされる可能性もあるし、適切なタイミングで呼ばせるというのは多少負担になる。 微々たるものとはいえ、消しされる負担は消し去るべきだ。

また、他にエラーメッセージ書き込み用のstd::string&を受け取っておく、というのも策としてあり得る。 ただ、私は関数の引数は少ないほうが良いと思っており、必要なステップも短いほうがいいと思っているので、これはあまり個人的には好ましくない。 それにこのやり方だと成功する場合でも関係なくstd::stringの初期化コストがかかってしまう。

解決策

で、どうするか。

要するに以下のように書ければ良いわけだ。

auto result = draw_rectangle(win, box);
if(result) {
    /* go next step */
} else {
    std::cerr << "an error occured!: " << result.unwrap_error();
}

または、

auto result = draw_rectangle(win, box);
if(result.is_error()) {
    std::cerr << "an error occured!: " << result.unwrap_error();
}

このインターフェースなら、expectedでいいのではないか。 variantを空にする場合、boost::blankstd::monostateなどを入れると良い。その気持ちで行けば、単純に

expected<boost::blank, boost::string_view>
draw_rectangle(window& win, const rect& rectangle);

のようにすれば良いのではないだろうか。

成功した場合は無を返し、失敗した場合はエラーメッセージを返す。普通にexpectedを使うことを考えた場合、単純で最もわかりやすい方法な気がする。

コンパイラの最適化と専用命令について

Fast inverse square rootなるものをご存知だろうか。これは、1.0 / sqrt(x)を高速に近似計算する命令である。 Fast inverse square root - Wikipedia

この計算は非常によく使う。例えばベクトルがあるとする。このベクトルを規格化したい。つまり長さを1にしたい。どうすればよいか。 当然、各要素をベクトルの長さで割れば良い。これ以上簡単なことはない。

inline std::array<float, 3> normalize(std::array<float, 3> const& v)
{
    const auto l = length(v);
    return {{v[0] / l, v[1] / l, v[2] / l}};
}

だが、一般に除算は他の演算に比べて遅い。依然、「アセンブリ解読」においてコンパイラが必死で割り算を消し去っていたことは記憶に新しい。 なので、(即値割り算をあれだけいじくり回すなら、多分コンパイラはオプションによっては上のコードでも勝手にやるだろうが)割り算をなくそう。

inline std::array<float, 3> normalize(std::array<float, 3> const& v)
{
    const auto inv_l = 1.0 / length(v);
    return {{v[0] * inv_l, v[1] * inv_l, v[2] * inv_l}};
}

ところで、length()の実装はどうなっているだろうか。

inline float length(std::array<float, 3> const& v)
{
    return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}

まあこんな感じだろう。ということはだ。

const auto inv_l = 1.0 / std::sqrt(...);

ここに、fast inverse square rootが出てきている。 この手の演算は、CGやゲーム・科学技術計算で死ぬほど出てくる。精度が要求される場合は微妙だが、近似計算で間に合う場合はこの計算はできるだけ早く行いたい。

という背景を思えば、以下の謎の関数が存在する理由もわかってもらえるのではないだろうか。

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

Wikipediaからの孫引きで、それによると引用元は"quake3-1.32b/code/game/q_math.c". Quake III Arena. id Software. Retrieved 2017-01-21."だ。 途中で謎のbit演算とマジックナンバーが出現しており、WTFどころではない。

とまあこんな手法まで登場するほどこの処理は本当に使用頻度が高くそれにしてはコストも高いので、intel製のCPUには専用命令を搭載している(rsqrtss)。 この命令の中身は恐らく上記のような形になっているのだろう。これが使える場合は普通にこれを使ったほうがよいので、マクロによる分岐が無限に生じていく。可読性は消え去る。

で、思ったことがある。素人が適当に最適化するよりも、無言でコンパイルオプションを-Ofastにしたら十分速くなるとよく言われる。 さらにClangはLLVMの恩恵を十全に活かして、まるでコードの意味を理解しているかのようなアグレッシブな最適化をするということで有名だ。 Clangがこの専用命令を知っている場合、1/sqrt(x)はその命令を使う形に勝手に最適化できないだろうか。 使うにしてもこれは近似計算なので、オプションで何も言わなければ使わないことも考えられる。どこまで教えれば使うのだろう。

というわけで試してみる。手元のマシンはhaswell世代のCPUで、Ubuntu 16.04、clang-4.0だ。以下のようなコードを書く。

#include <math.h>
float rsqrt(float x)
{
    return 1.0 / sqrtf(x);
}

とりあえず上記の最適化がなされるか確かめるため、いつものやつでいく。

$ clang-4.0 -std=c99 -ffast-math -march=native -mtune=native -Ofast rsqrt.c -c -o rsqrt_clang.o
$ objdump -d rsqrt_clang.o > rsqrt_clang.asm

というわけで出たのが以下のコードだ。

0000000000000000 <rsqrt>:
 0:  vrsqrtss    %xmm0,    %xmm0,%xmm1
 4:  vmulss      %xmm1,    %xmm0,%xmm0
 8:  vfmadd213ss 0x0(%rip),%xmm1,%xmm0
11:  vmulss      0x0(%rip),%xmm1,%xmm1
19:  vmulss      %xmm0,    %xmm1,%xmm0
1d:  retq

少なくとも、専用命令が使われている。では-march=native-mtune=nativeを外してみよう。

0000000000000000 <rsqrt>:
 0:  rsqrtss    %xmm0,    %xmm1
 4:  mulss      %xmm1,    %xmm0,
 8:  mulss      %xmm1,    %xmm0,
 c:  addss      0x0(%rip),%xmm0
14:  mulss      0x0(%rip),%xmm1,
1c:  mulss      %xmm1,    %xmm0
20:  retq

どうやら効いていたのはFMAだけだったっぽい。-ffast-mathを外してみる。この辺りでこの命令は使われなくなるのでは……と思ったが、先のアセンブリコードと完全に一致した。

で、-Ofast-O3にしたところ、ついにrsqrtssが消えた。

0000000000000000 <rsqrt>:
 0:  movaps %xmm0,%xmm1
 3:  xorps  %xmm0,%xmm0
 6:  sqrtss %xmm1,%xmm0
 a:  ucomiss %xmm0,%xmm0
 d:  jnp    1c <rsqrt+0x1c>
 f:  push   %rax
10:  movaps %xmm1,%xmm0
13:  callq  18 <rsqrt+0x18>
18:  add    $0x8,%rsp
1c:  movss  0x0(%rip),%xmm1
23:  
24:  divss  %xmm0,%xmm1
28:  movaps %xmm1,%xmm0
2b:  retq   

もしかしてOfastだと-ffast-mathが勝手に追加されてるのかな(うろ覚え)と思い、-Ofast -ffast-mathでやってみると、復活した。man clangをチラ見したもののどのフラグが自動で立つかは見当たらなかった。 やはり-ffast-mathによる最適化を行えば、rsqrtssが使われるらしい。恐らく理由は、これが近似計算だからだろうと思う。

コンパイラも結構やってくれるんだな、ということがわかった。ちなみにgcc-7.2.0-ffast-math -O3で同等のバイナリを吐いた。 ここまでずっと32bit浮動小数点数でやってきたが、64bitでやると普通にsqrtを取って割るコードになる。恐らく対応する命令がないのだろう、とは思うが、どうなのだろう。

コードの美しさについて

コードが綺麗とはどういうことか、の基準は人によって異なると思う。しかし例えば以下のような関数があまり綺麗ではないという気持ちは、共有できるのではないだろうか。

void func(int N, double*** v1, double*** v2, double** a, double** b, int* c, int* d, int* e, int f, int g, int h);

コードの場合、美しいとは、シンプルであることだ。見ただけですぐに何をしたいかがわかること、要素ごとに纏められていること、名前や型が適切に付けられて意図が伝わりやすいこと、スタイルが一貫していること、そして、可能なら、短いこと。 間違えようのない最小限のことを正しくこなすコードが相互作用しつつ大きく複雑なものになっているというのが理想的に思える。

構成要素をシンプルに保っていれば(そして、その構成要素に適切な名前がついていれば)、可読性、保守性、開発効率はどれも上がる。 小さく単純なものは理解しやすいので可読性が高い。可読性が高いものは多くの人が理解しやすいので修正もききやすい。当然理解しやすいものは書き足しやすいし、それ一つで完結しているものは切り出して他のところへ組み込みやすい。 非常に大雑把に言ってしまえば、そうなっているもののことを美しいコードと呼ぶのだ。

続きを読む

数値積分のインターフェース

子細合ってヤバい形の積分を解く羽目になったので数値積分をすることにした。 別に数値積分の複雑なアルゴリズムの解説や異常なまでの最適化などを紹介するわけではなく、非常に基本的な台形積分を使い、そこまで最適化も考えずにやる。 だがそれだけだと何の面白みもないので、使い回すことを考えた時に数値積分を実行する関数のインターフェースをどうするか少し考えてみたい。

続きを読む