フィボナッチ数列を再帰関数で計算すると遅い

よく知られた事実ですが、指数的な回数の再帰呼び出しが行われるため、実用的ではありません。


フィボナッチ数列を計算する再帰関数は再帰の仕方がわかりやすいですが、よくわからない再帰の仕方の関数の場合はどうでしょうか?

関数が呼ばれたときにデバッグプリントをする方法では、どのような再帰構造で呼ばれたのかわかりません。

再帰深度だけインデントしてデバッグプリントすればよいのですが、インデント量を保持する必要があります。

一つの方法は、再帰関数に、再帰深度を受け取る引数を追加する方法です。その引数には、外部から呼び出す場合は0を、再帰呼び出しをする場合は、引数の値+1を、それぞれ渡します。 この方法の欠点は、関数原型を変更する必要があるです。また、相互再帰の場合は関係するすべての関数の関数原型を変更する必要があります。

もう一つの方法は、グローバル変数再帰深度を保存する方法です。関数の先頭でその変数をインクリメントし、return文の直前でその変数をデクリメントします。 自己再帰関数であれば、グローバル変数とせず、関数内static変数とすることもできます。 この方法の欠点は、大域脱出されるとデクリメントする機会を失い、値が永久に狂ってしまう点です。

前者は再帰段数に比例する記憶領域を消費するのに対し、後者は定数の記憶領域しか消費していません*1。 これは不思議な感じがありますが、(末尾再帰ではない場合)後者でも再帰段数に比例するコールスタックが消費されていることを合わせて考えると、定数倍の違いになっています。

後者の方法の空間効率が良いのは、「デクリメントすれば元に戻る」、つまり逆操作が存在するという部分が重要です。 木を走査しながら演算を繰り返していき、葉での値を求める、といった場合、一般に逆演算があるとは限りません。 逆演算がない場合、後者のような方法はとれず、前者のような方法を採用することになります。

なお、もともとが末尾再帰だった場合、後者の方法は末尾再帰ではなくしてしまう効果があります。 つまり、先ほどとは逆に、前者の方法は定数の記憶領域のみを消費し、後者の方法は(コールスタック用に)再帰段数に比例する記憶領域を消費することになります。 ただし、コールスタックに積まれる情報は毎回同じものなので、後者の方法の場合でも本質的には圧縮可能です。

*1:実際は、再帰段数の対数桁の情報を覚える必要がありますが、いずれにしても再帰段数の線形よりは少ないです

32bit版RISC-VのLinux syscall 62はlseekではない?

Linuxシステムコール62番は、64bit版の場合、lseekです。

off_t lseek(int fd, off_t offset, int whence);

lseekシステムコールは、第一引数にファイルディスクリプタ、第二引数にどれだけ移動させるか、第三引数に動作の種類を受け取り、最終的な位置を返り値で返します。 第二引数の解釈は、第三引数によって変わります。エラーが発生した場合は、返り値が(off_t)-1になり、errnoによってエラーの種類が通知されます。

Linuxシステムコールの場合、引数や返り値はレジスタで渡されます。そうすると、32bit版では、2GiB程度までのファイルしか扱えないということになりそうです。

調べてみると、Linuxシステムコール62番は、32bit版の場合、_llseekになっています。

int _llseek(unsigned int fd, unsigned long offset_high, unsigned long offset_low, loff_t *result, unsigned int whence);

_llseekシステムコールは、第一引数にファイルディスクリプタ、第二引数と第三引数にどれだけ移動させるか、第五引数に動作の種類を受け取り、最終的な位置を第四引数で指定されたポインタの部分に書き込みます。 第二引数と第三引数は組み合わせて64bitの値として使用されます。エラーが発生した場合は、返り値が-1になり、errnoによってエラーの種類が通知されます。成功した場合、返り値は0になります。

このようにすることで、32bitのレジスタしかもっていないアーキテクチャであっても、2GiBを超えるファイルが取り扱えるようです。

この間ビルドしたspike、qemu3.1.0を使って調べたところ、どちらもlseekの動作をしているのですが、qemu4.1.0を使ってみると_llseekになっています。 gccコンパイル結果を見ても、_llseekを期待しているバイナリになっています。

これはシステムコールの仕様が変わったのでしょうか……? 公式が作っているシミュレータであるspikeが常に正しいと思っていたのですが、古い常識となっているのでしょうか?

ちょっとはまったのでメモ書きでした。

パック展開とバイトニックソート

結構前の、constexpr関数を高速化する方法 - よーるという記事で、低速なconstexpr文脈の中でもある程度の速さで配列操作を行う方法として、パック展開を利用した配列を一気に構築する手法を紹介しました。

これを利用したバイトニックソートは、たとえば私の作った lfl/lfl.hpp at master · lpha-z/lfl · GitHub (171行目から201行目)があります。

私が書く前に存在した実装としては、bitonic_sort.hpp · GitHubなどでしょうか。

ところで、パック展開を用いると、シャッフル的なことが簡単にできます。

#include <array>
#include <utility>

template<class T, std::size_t N>
constexpr T comparator( const std::array<T, N>& arr, std::size_t i, std::size_t j ) {
    return i < j ^ arr.at(i) < arr.at(j) ? arr.at(i) : arr.at(j);
}

template<class T, std::size_t N, std::size_t... Indeces>
constexpr std::array<T, N> sorter_impl( const std::array<T, N>& arr, std::index_sequence<Indeces...> ) {
    return std::array { comparator( arr, Indeces, Indeces^1 )... };
}

template<class T, std::size_t N>
constexpr std::array<T, N> sorter( const std::array<T, N>& arr ) {
    return sorter_impl( arr, std::make_index_sequence<N>{} );
}

template<std::size_t... Indeces>
constexpr auto swapper = []( const auto& arr ) { return std::array { arr[Indeces]... }; };
    
constexpr auto sort = []( auto arr ) {
    arr = sorter( arr );
    arr = swapper<0,3,1,2,4,7,5,6>( arr );
    arr = sorter( arr );
    arr = swapper<0,2,3,1,4,6,7,5>( arr );
    arr = sorter( arr );
    arr = swapper<0,7,1,6,2,5,3,4>( arr );
    arr = sorter( arr );
    arr = swapper<0,4,2,6,7,3,5,1>( arr );
    arr = sorter( arr );
    arr = swapper<0,2,1,3,4,6,5,7>( arr );
    arr = sorter( arr );
    return arr;
};    

#include <iostream>
int main()
{
    constexpr auto arr = sort( std::array { 1, 2, 5, 4, 6, 3, 7, 0 } );
    for( auto e : arr ) { std::cout << e << std::endl; }
}

上のコードの、swapperという関数で、配列の要素をシャッフルしています。関数ではなくラムダ式を使っているのは、引数のところをautoと書きたかったからです。こうすることで、template<class T, std::size_t N>のように書かなくて済むようになり、swapper関数のテンプレート引数部分にシャッフルする順番のみを書くことができます。

もちろん、このようなコードは配列の書き換えが多く非効率的です。 ソートした結果をシャッフルする部分はひとまとめにできるので、以下のように改善が可能です。

#include <array>
#include <utility>

template<class T, std::size_t N>
constexpr T comparator( const std::array<T, N>& arr, std::size_t i, std::size_t j ) {
    return i < j ^ arr.at(i) < arr.at(j) ? arr.at(i) : arr.at(j);
}

template<std::size_t... Indeces>
constexpr auto swapper = []( const auto& arr ) { return std::array { comparator( arr, Indeces, Indeces^1 )... }; };
    
constexpr auto sort = []( auto arr ) {
    arr = swapper<0,3,1,2,4,7,5,6>( arr );
    arr = swapper<0,2,3,1,4,6,7,5>( arr );
    arr = swapper<0,7,1,6,2,5,3,4>( arr );
    arr = swapper<0,4,2,6,7,3,5,1>( arr );
    arr = swapper<0,2,1,3,4,6,5,7>( arr );
    arr = swapper<0,1,2,3,4,5,6,7>( arr );

    return arr;
};    

#include <iostream>
int main()
{
    constexpr auto arr = sort( std::array { 1, 2, 5, 4, 6, 3, 7, 0 } );
    for( auto e : arr ) { std::cout << e << std::endl; }
}

このようにしてもアクセスする順序が乱れているため高速ではないでしょう。 このようなバイトニックソートの亜種は大量に存在します。

それにしても、(ハードコードされているとはいえ)こんな単純なコードでバイトニックソート(の亜種)が書けるパック展開は便利ですね。

Brainfuckのインタプリタを書いた

Brainfuckインタプリタの書き方によって、実行速度が変わるのかが気になったため書いてみました。

もちろん、前処理、特にその場でコンパイル (just in time compile) をするなどすれば高速化は容易です。 そうではなく、ごく単純なswitch文を使うような、一命令ずつ見ていくタイプのインタプリタでも、書き方によって速度に大きな差が出るかという疑問です。

(余談)ファイル丸ごと読み込み

求めていたものは以下のページに書かれています。適切な検索キーワードを見つけることが難しいですが、「string fstream 一気に読み込み」で見つかりました。

C++ ファイルを全て読み込む - 備忘録

int main( int argc, char* argv[] ) {
    std::ifstream ifs( argv[1] );
    std::string str( (std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>() );
}

(std::istreambuf_iterator<char>(ifs))とかっこがついているのには意味があります。これがないと、以下のパターンにマッチしてしまうため、関数宣言だと解釈されてしまいます。

ReturnType func_name( Arg1Type (arg1), Type2() );

(arg1)の部分は不自然ですが、このようなかっこをつけることは許されています。

Type2()の部分は、Type2 (*)()、つまり関数ポインタ型(引数なし、返り値型がType2)が第二引数の型であるというように解釈されます。 第二引数の名前は省略されています。

丸かっこを使ったコンストラクタ呼び出しが関数宣言と解釈される現象があることは知っていたのですが、引数なしの場合だけだと思っていました。第一引数のように無駄なかっこ扱いされる場合、第二引数のように引数なし関数のかっこ扱いされる場合はどちらも知りませんでした。教育的な例ですね……。

このように丸かっこを使ったコンストラクタ呼び出しは問題が発生しやすいです。 こういう場合は統一初期化構文を使うのが良い解決策です。

例えば、std::istreambuf_iterator<char>()ではなく、統一初期化構文を使ってstd::istreambuf_iterator<char>{}とすれば問題が解決します。

自然な実装

普通にインタプリタを書くと以下のようになります。データメモリをスタック上にとるとスタックオーバーフローが発生しうるので避けたほうが賢明です グローバル領域に確保するのが最も簡単です(よりきれいな書き方は、動的確保する方法でしょう)。

#include <cstdint>
#include <array>

std::array<std::uint8_t, 1000000> data;

#include <iostream>
#include <fstream>
#include <string>
int main( int argc, char* argv[] ) {
    std::ifstream ifs( argv[1] );
    std::string insn( (std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>() );
    auto ip = insn.begin();
    auto dp = data.begin();
    while (ip != insn.end()) {
        switch(*ip) {
        case '+': ++*dp; break;
        case '-': --*dp; break;
        case '>': ++dp; break;
        case '<': --dp; break;
        case '.': std::cout << *dp; break;
        case ',': std::cin >> *dp; break;
        case ']': if ( *dp) for( int loop = 0; loop += *ip == ']', loop -= *ip == '['; ) --ip; break;
        case '[': if (!*dp) for( int loop = 0; loop += *ip == '[', loop -= *ip == ']'; ) ++ip; break;
        default:;
        }
       ++ip;
    }
}

もっと速い実装

実は、以下のような意味不明なswitch文を書くと、先の実装に比べて1.37倍ほど高速になります。 速度の測定は、Intel Core i7-4771 (Haswell)、g++-9.1 -O3 -march=native、動作プログラムはErik Bosman氏作mandelbrot.bです。

#include <cstdint>
#include <array>
std::array<std::uint8_t, 1000000> data;
#include <iostream>
#include <fstream>
#include <string>
#include <cstdio>
int main( int argc, char* argv[] ) {
    std::ifstream i( argv[1] );
    std::string insn( (std::istreambuf_iterator<char>(i)), std::istreambuf_iterator<char>() );
    auto ip = insn.begin();
    auto dp = data.begin();
    while (ip != insn.end()) {
        switch(*ip) {
        case '>': ++dp; break;
        case '<': --dp; break;
        default:;
        switch(*ip) {
        case '+': ++*dp; break;
        case '-': --*dp; break;
        case '.': std::cout << *dp; break;
        case ',': std::cin >> *dp; break;
        case ']': if ( *dp) for( int loop = 0; loop += *ip == ']', loop -= *ip == '['; ) --ip; break;
        case '[': if (!*dp) for( int loop = 0; loop += *ip == '[', loop -= *ip == ']'; ) ++ip; break;
        }
        }
        ++ip;
    }
}

なぜ速度に違いが出るのか

観察

case文を"上側"または"下側"のswitch文に移動した場合のコンパイル結果を観察します。すると、以下のような違いがあることがわかります。

  • case '['case ']'に対応する文は、どちらに入れても大きく変わらない
    • case '['case ']'に対応する文はループのある長いコードであり、コンパイラはこれを別立てで考えている?
  • その他のcase文6つのうち、5つ以上が片方に入ると、ジャンプテーブルを使った実装になる
    • case式の値が狭い範囲に分布しているかとは無関係にジャンプテーブルになる
    • case式の値が狭い範囲に分布していると、ジャンプテーブルのサイズは小さくなる
  • ジャンプテーブルにならない場合、高速に実行される
  • case '<'case '>'のみを"上側"のswitch文に入れると最高速
    • <>は出現頻度が最も高い二つ
    • 一般に、出現頻度の高い事象から検討するのが最適

速度に差が出る原因

一般にジャンプテーブルを使った間接分岐でswitch文を構成する場合、分岐予測が当たりにくくなります。

分岐予測器は条件分岐(if文)が成立するか否かを当てることには強いです。一方、その分岐が飛ぶ先を当てることは苦手であり、前に飛んだアドレスにもう一度飛ぶだろうくらいの予測しかしてくれません(最近は間接分岐用の分岐予測器も研究されているようです。この前、パーセプトロンを使ったような間接分岐予測器が発表されていました)。

ジャンプテーブルを使う場合、そのswitch文を処理するためには分岐予測ミスは高々一回しか起こりません。 一方、ジャンプテーブルを使わない場合、(おそらく二分探索的な)if文のツリーになるため、分岐予測ミスは複数回起こりえます。

こう聞くとジャンプテーブルのほうがよさそうに思えますが、今回作ったBrainfuckインタプリタは以下の条件が効いてくるため、if文のツリーの分岐予測ミス率は意外と低いものとなり、ジャンプテーブルを使った場合の性能を上回ります。

  • インタプリタを動かす場合、同じような命令列が流れてくることが多い。よってswitch文の引数に来る値の系列には法則性があり、過去の分岐方向履歴を活用した分岐予測は当たりやすい
  • 一見256方向分岐だが、実行時間のほとんどを占めているのは+-<>[]の6文字だけである(,は現れず、.もほとんど現れない)。しかもほとんどの部分が<>の連続である。つまり、「前回と同じ方向に行く」という雑な分岐予測器でさえ、かなりの精度で当てることができる

g++-9ではコード生成エンジンが大幅に変わっており、特にswitch文の最適化を強化したとうたっていました。 今回のような結果になったのは、g++-9を使ったことも原因の一つでしょう。

g++-8との違い

g++-8を使った場合、case文を4つずつに分けたときだけはif文のツリーになりました。つまり、「4つ以下ならif文のツリーにする」というルールは変わっていないようです。一方、g++-9では、case '['case ']'に対応する文のような長いコードを特別扱いするようです。長いコードの場合は分岐予測ミスによる損失が相対的に小さいとみなしているのでしょうか……?

まとめ

g++-9ではswitch文の最適化が強化されており、小型のswitch文をif文のツリーに分解する最適化をしてくれることがあります。 よって、switch文を二つに分解するという一見意味の分からないコード変形を行うことで、if文のツリーへの変形を行ってほしいということをコンパイラに伝えることができます。

インタプリタを実装するときは、ジャンプテーブルではなくif文のツリーを使ったほうが高速な場合もあります。 よって、このようなコード変形は、インタプリタの高速実装に役立てられるかもしれません。

32bit版RISC-V用のspikeを動かす

オープンなISAであるRISC-Vには、32bit版も定義されています。 spikeは、RISC-Vの公式が開発している命令セットシミュレーターです。

以下の手順で必要なソースコードを持ってきます。

$ git clone https://github.com/riscv/riscv-tools.git
$ cd riscv-tools
$ git checkout bce7b5e
$ git submodule update -i --recursive

最新のコミットでは、riscv-gnu-toolchainが含まれなくなってしまいました。 RISC-V用のクロスコンパイラをすでに持っている場合は問題ありませんが、そうでない場合は手に入れる必要があります*1riscv-gnu-toolchainが含まれていたころのbce7b5eをチェックアウトすれば自動でやってくれるため楽です。

このあと./build-rv32ima.shを動かすのですが、名前の通り、RV32IMA用になってしまいます。 RV32IMAだと浮動小数点数命令拡張が含まれないため、スクリプトの中のrv32imarv32gに変更します。

インストールしたいパスをシェルのRISCV変数に入れた後、先ほど書き換えたビルドスクリプトを走らせます。

$ export RISCV=~/rv32
$ ./build-rv32ima.sh

70分ほど待つとビルドが完了します。

使い方は、以下のような感じです。前に作ったRV64G用の場合は~/rv64/bin/spike pk a.outのようにpkをそのまま打てばよかったような気がしますが、~/rv32/bin/spike ~/rv32/riscv32-unknown-elf/bin/pk a.outのようにpkのパスを打つ必要があるようです。

$ echo -e '#include<stdio.h>\nint main(){puts("Hello, RISC-V!");}' > main.c
$ ~/rv32/bin/riscv32-unknown-elf-gcc main.c
$ ~/rv32/bin/spike ~/rv32/riscv32-unknown-elf/bin/pk a.out
bbl loader
Hello, RISC-V!

*1:OSの代わりであるproxy kernelをRISC-V用にビルドする必要があるため、クロスコンパイラが必要になります

忙しく感じる

なんか最近忙しくないはずなのに、忙しいような気がしています。

たぶん原因は、塊として認識できない、こまごましたやらないといけないことが多いことのような気がします。 あとは、それをやることのできる時間が限られている(持ち時間が少ないという意味ではなく、昼限定イベント的な感じ)のが原因のような気がします。

坂井先生もおっしゃっていましたが、やはりOSのスケジューリング問題だと思ってやるのがよいようです。