satysfi.vimを作った

というわけで中身を読もうとやっていっていたわけだが、SATySFiにはVim向けのシンタックスハイライターがない。当たり前だ。新興の言語なのだから。だが読むに当たってシンタックスハイライトがされないというのは非常にしんどい。なので最低限のハイライトがされるようなスクリプトを書くことにした。 気持ちとしては、SATySFiを使えるようになりたいがシンタックスハイライトがないと読みにくい、しかしSATySFiについて詳しくないのでどこをハイライトしたらいいかわからない、というデッドロックを破るための、ブートストラッピングの最初の一手というところである。

まだとても完成と言える状態ではないが、最低限度の質素なハイライトはされるので、自分でやるのが面倒だという方は使ってみて欲しい。そしてあなたがハイライトに物申したいときは遠慮なくissueなりPRを送っていただきたい。

github.com

これが突貫工事であることを除いても、私はまだvimシンタックスファイルにもOCamlにもSATySFiの構文にも詳しくないので、ハイライトされるべき箇所がされていない可能性は高いし、それを見つけられる可能性も高い。

まあ、vimの構文ファイルもSATySFiも知らない人間が1時間で作ったということに思いを馳せれば、たいていのことは許していただけると思う。


ところで私は今までVim用のシンタックスファイルを書いたことがない。作り方も導入の仕方も知らない。今まで自作言語でも作っていれば経験があっただろうが、仕方のないことだ。調べなければならない。

vimには素晴らしいドキュメントがある。しかも日本語だ。

syntax - Vim日本語ドキュメント

あるいは、こちらでもよい。

Creating your own syntax files | Vim Tips Wiki | FANDOM powered by Wikia

まず、vimの普段のシンタックスハイライトは、c.vimのようなファイルがあることがわかる。単純なコードならこれを見ることで色々とわかるかもしれないが、とても読めるような代物ではない。私と同様の環境なら、/usr/share/vim/vim??/syntax/*.vimのような場所にシンタックスファイルが死ぬほど転がっているが、これを材料に独学するのは非常に難しい。

素直にドキュメントを読み進めることにしよう。ドキュメントの構文ハイライトファイルのセクションを見れば、ある程度のことがわかる。

しかし、書き始める前に確認できる状態を整えねばならない。最終的にはvimのパッケージマネージャで管理できるようにするべきだが、書いている途中に毎度updateするのはしんどい。それ以上に、パッケージマネージャの仕組みを今から調べる時間はない。

よって、一番原始的な方法で実行する。まず、~/.vim/ディレクトリにftdetectsyntaxというディレクトリを作る。まずftdetectsatysfi.vimというファイルを置き、その中に以下を書き込む。

autocmd BufNewFile,BufRead *.saty,*.satyh set filetype=satysfi

ftdetectは file type detect の略で、これによって*.satyファイルや*.satyhファイルを開いた時にファイルタイプが設定されるらしい。 ファイルタイプがわかったら、シンタックスファイルを選べる。syntaxディレクトリにsatysfi.vimなるファイルを設置する。で、とりあえず以下を書こう。

if exists("b:current_syntax")
  finish
endif

これは、既にシンタックスが定義されていた場合に終了するためのものだ。やっておいた方が行儀がよいだろう。

では、まず手始めにコメントがコメントと認識されるようにしよう。

syn match satysfiComment /%.*/
hi def link satysfiComment Comment

一行目がsatysfiCommentという構文要素の定義だ。名前に続く正規表現、つまり%から始まる文字列にマッチする場合、それはsatysfiCommentである。 続いて、hi def linkは構文要素をCommentというハイライトグループに紐付ける。 カラースキームはこれらの構文要素に色を設定していくものだろう、と思っている(作ったことがないのでこれは予想である)。

グループについてはドキュメントの以下の部分が、

syntax - Vim日本語ドキュメント

vimで使える正規表現は以下の記事が参考になる。

vim正規表現リファレンス - Qiita

基本的には、ハイライトしたいと感じる部分を見出し、一応Lexerを見つつどういったものか確認して、マッチする正規表現を書く、という形でやっていった。 括弧ごと何かしたい場合などは、syn region start="" end=""が、単なるキーワードにはsyn keywordが使える。


困ったのはキーワードのinで、これを単にkeywordにすると、例えばプリミティブのtext-in-mathなどの真ん中のinが誤ってハイライトされる。 普通にやるとケバブ・ケースの命名規則はサポートされないのだ。 ケバブ・ケースといえばLispなのでLispの構文ファイルを少し見に行ったが、難しすぎたので逃げ帰ってきた。 だがこの位なら自力でなんとかできるだろう、と思いしばらくsatysfiパッケージのソースなどを読んでいたところ、in直後に改行もしくは空白のマッチでよいのではという気がしてきたのでそのようにした。

とまあ、いくつかの困ったところに関しては死ぬほどアドホックなことをしているので、まだ色々変わると思われる。 順調にSATySFi(とvimの構文ファイル)に詳しくなれればより改善されていくと思うので、暖かく見守ってほしい。

SATySFiのインストール

SATySFiは静的型検査の恩恵を受けられる組版処理システムである。

github.com

Ubuntu16.04を使っているのだが、SATySFiを使ってみようとして少し困ったので。ちなみにそれまでのOCaml経験はHello, Worldのみ、つまりコンパイラを入れたことがある以外はゼロ。むしろデファクトな環境構築法を知らないまま適当な環境を導入した分マイナスかも知れない。

先に結論から行くと、以下のようにすると上手く行った。シェルはfish-2.7。opamがfishにも対応していて少しうれしくなった。普段使う大抵のソフトウェアはfishに対応していないので。

$ sudo apt-get install opam
$ opam init # 少し時間がかかる。configを弄るのはokした
$ eval (opam config env)
$ opam switch list
... # 確認のため
$ opam switch 4.05.0 # 時間がかかる
$ eval (opam config env)
$ git clone https://github.com/gfngfn/SATySFi.git
$ git submodule update -i
$ opam pin add satysfi .
$ opam install satysfi
$ satysfi --version # 確認
  SATySFi version 0.0.1

最初、出ないはずのシンタックスエラーに遭遇してわけがわからなくなっていた。opamのインストールから始めると上手く行ったので恐らく最初の環境構築時点で何かやらかしていたのだとは思うが、よくわからない。最初のopamのインストールは公式のinstallスクリプトwgetしてシェルに流し込むやつでやったのだが、素直にaptで入れると上手く行った。もちろんopam initとかを忘れていただけの可能性はある。

ちなみにopamをソースからビルドしようとしたところ、opamが依存しているライブラリをダウンロードしようとして、md5sumが違うと言って終了するという事件があった。もちろん確認せずに勝手に解凍されては困るのでこの処理は疑いようもなく正しいのだが、こちらとしてはどうにもならないので辛かった。素直にaptで入れよう。

ちなみに、デモやドキュメントをコンパイルしようとするとフォントがないので落ちる。これに関してはフォントを持っていないこちらの問題なわけだが、フリーなフォントを使用するようにしたブランチとそれを入れていくDockerfileがあるので、それを参考にやっていく。

github.com Comparing gfngfn:master...pandaman64:use-free-font · gfngfn/SATySFi · GitHub

ちなみに、インストール前なら上の比較を見ながら変えればよいが、インストール後は~/.satysfiを見に行くのでそちらの同名ファイルを書き換えなければならない。

というわけでPDFが作成できるようになる。

しばらくデモとドキュメントのPDFとソースを比べながら手さぐりしていけばそのうちなんとかなるだろうと考えているが、果たして。

基本ベクトルの場合の最適化(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を取って割るコードになる。恐らく対応する命令がないのだろう、とは思うが、どうなのだろう。