最大公約数をもっと高速に求める(その3)

先週の記事の続きです。

最大公約数をもっと高速に求める(その2)【cmova命令は遅い】 - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = t < m;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

このコードの機械語命令列

gcd_impl関数のループ一回分の機械語命令列は、以下のようになっています(実際にはループアンローリングされており、レジスタ割り付けだけが異なる命令列が3回周期で現れます)。

sub    rdi, rsi    # t = n - m;
cmp    rdi, rsi    # q = t < m;
mov    rax, rdi
cmovb  rax, rsi    # n = q ? m : t;
cmovnb rdi, rsi    # m = q ? t : m;
test   rdi, rdi    # m == 0
jz     Label       # if( m == 0 )

ループ一回分は7命令になっているわけですが、だからといって7 cycleかかるわけではありません。近年のCPUでは、これくらいの命令列は並列に実行されます。

今回の命令列で、並列実行のボトルネックになるのは依存の連鎖です。 具体的には、t = n - mを計算し、その結果を使ってq = t < mを計算し、その結果によって次ループのnmが算出される部分です。

qtが算出されなければ計算できませんし、nmqが算出されなければ計算できません。また、次のループのtnmが算出されなければ計算できません。

つまり、どんなに並列に実行できるCPUを持ってきても、この部分は各々1cycleかかってしまいます。

よってこのコードでは、ループ一回分に3cycleかかることになります。

なお、mov命令が気になるかもしれませんが、最近のCPUでは実質的に0cycleで実行されると考えてよいようです。

互減法をさらに高速化する

互減法の一ステップでは、以下のように計算が進みます。

GCD(n, m) = if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

これを拡張し、以下のようにしてみます。

GCD(n, m) = if n-Km >= m then GCD(n-Km, m) else if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

もしn-Km >= mであれば、互減法のKステップを一気に進めることができるため、互減法の効率が上がることが予想されます。 また、それにより、アルゴリズム全体に占める除算の割合が減少し、高速化すると予想されます。 ただし、以下の点に注意する必要があります。

  1. 途中結果n-Kmが負になる可能性がある。
  2. Kmがオーバーフローする可能性がある。

1.については、移項することで回避可能です。 2.については、前処理でmを小さくしておけば回避可能です。

改良されたメインループ

メインループの内部を以下のようにしてみます。メインループの前に前処理を置くことでうまくmは小さくしているとして、m*Kはオーバーフローしないと仮定します。

uint64_t t = n - m;
uint64_t s = n - m * K;
bool q = t < m;
bool p = t < m * K;
n = q ? m : t;
m = q ? t : m;
if( m == 0 ) { return n; }
n = p ? n : s;

n = p ? n : sが追加されたため、さらに依存の連鎖が伸び、このループ一回には最低でも4cycleかかるようになりました。 もし、m*Kの計算に4cycle以上かかる場合、そのせいで依存の連鎖が伸び、ループ一回にかかるサイクル数はもっと長くなります。 あるいは、m*Kの計算が他の命令の実行の裏で並列に実行できない場合も、ループ一回にかかるサイクル数はもっと長くなります。

最近のIntel社製CPUであれば、MUL命令は3cycleで完了するため、m*Kの計算のせいで依存の連鎖が伸びるということはありません(Kの値によってはコンパイラ最適化によりLEA命令を駆使した命令列となり、もっと速く計算されます)。 また、m*Kの計算は、他の命令の実行の裏で行うことができるため、そのせいでループ一回にかかるサイクル数が長くなるということはありません。

元のコードのループ一回には3cycleかかっていたところ、4cycleかけることで、Kステップ一気に実行できる可能性が追加されたことになります。 nmの大きさに大きな差がなかった場合でも、少なくとも1ステップは進むため、それほど損したことにはなりません。

Kとして何をとるか

Kとして2や3などの小さな数をとれば、n-Km >= mである可能性がそれなりに高いため、改良により互減法のステップを効率よく進めることができるようになります。 Kとして大きな数をとる場合、n-Km >= mである可能性が低く、意味のある仕事をできる確率(pfalseになり、n = sが実行される確率)が下がります。しかし、意味のある仕事をできた場合、互減法のステップを一気に進めることができます。

このトレードオフを考慮に入れ、いくつかの実験を行った結果、以下にようにすることでstd::gcdより二倍以上高速なgcd関数を作れることがわかりました。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

gcd_pre関数は、m*Kがオーバーフローしないようにするための前処理です。互減法を4回繰り返すことで、mは大きくとも0xffff'ffff'ffff'ffff / 5であると保証できます。

Kとして5をとったことにより、32のような大きな商が立つ場合でも、6回程度ループを回れば互除法の一ステップを行ったことになります。これ以上大きな商が立つ場合は、除算を行った方がよさそうです。 そのようなことの発生頻度が低いため、単純な互減法であれば10回程度で除算に切り替えたほうがよかったところ、80回まで増やすことが可能になりました。 なお、gcd_implのループ回数80回は目安であり、必ずしも最適な値とは限りません。

最大公約数をもっと高速に求める(その2)【cmova命令は遅い】

先週の記事の続きです。

最大公約数をもっと高速に求める(その1) - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = m > t;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

コンパイラの気まぐれにつきあう

先ほど示したコードの

                bool q = m > t;

                bool q = t < m;

に変えてみましょう。私の環境では、25%ほど高速化しました。わけがわかりません。こんな高速化は誰も思いつけません。

アセンブリを確認しても、レジスタの割り付けは一切変わっておらず、単にcmp命令のオペランドが逆転し、後続のcmov命令の条件が変わるだけです。

他の環境でも試してみましたが、25%などという劇的な高速化ではないものの、10%程度は性能が向上するようです。

Intel社製CPUとつきあう

Intel Architecture Code Analyzerを使って調べてみたところ、Broadwell以降の場合、

  • cmovb命令やcmovnb命令は1つのμOPに変換される。実行は一サイクルで完了する。
  • cmovbe命令やcmovnbe命令は2つのμOPに変換される。実行には合計で二サイクルかかる。
    • cmova命令はcmovnbe命令の別名称。
    • cmovna命令はcmovbe命令の別名称。

というようになっているようです*1

Code Analyzerなんて持ち出さなくても、インテルの最適化リファレンスマニュアルに書いてあったようです(最後の方、C-17ページ。pdfとしては759ページ目)。

https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf

bool q = m > tの場合、cmovbe命令とcmovnbe命令が使われ、レイテンシが二サイクルとなります。一方、bool q = t < mの場合、cmovb命令とcmovnb命令が使われ、レイテンシは一サイクルとなります。 これにより高速化したということのようでした。

なお、Haswell以前の場合、前者は2μOP、後者は3μOPに分解され、それらが直列に実行されるようなので、その差は小さくなるとはいえ、やはり性能に差が出ます。

まとめ

cmov命令は、使う条件によっては遅いことがあります。(はじめて知りました!) そのため、不等号の方向を変えるだけという一見等価な変換で高速化したり、遅くなったりすることがありえるようです。

今回の記事はもっとアルゴリズム的なところを掘り下げようと思っていたのですが、Intel社製CPUに最適化するだけの回になってしまいました。(でも25%高速化!)

まだまだ続きます。

*1:似たような命令なのに動作原理が違うあたり、x86はやはりCISCだなぁという感じです

最大公約数をもっと高速に求める(その1)

最大公約数を求める非常に有名なアルゴリズムとして、ユークリッドの互除法が知られています。 このアルゴリズムは、二つの入力を十進法で表したとき、小さいほうの桁数の高々五倍の回数の剰余演算で最大公約数が求まることが知られています(ラメの定理)*1。二つの入力が符号なし64bit整数であれば、高々92回です。これは最悪ケースであり、多くの場合はかなり少ない回数で済みます。実験してみると、平均的には37回程度のようです。

しかし、剰余演算は非常に遅い演算です。Intelの最近のCPUでは、21cycle~83cycleかかるとされています。ここから計算すると、二つの符号なし64bit整数の最大公約数を求めるには、1μ秒弱の時間がかかることになります。

実際、以下のコードを実行してみると、乱数生成を除いて7秒ほどかかります。つまり一回当たり0.42μ秒ということで、だいたい計算と合います。

#include <iostream>
#include <random>
#include <chrono>

// 16M
std::pair<uint64_t, uint64_t> arr[0x100'0000];

int main() {
        std::mt19937 mt;
        std::uniform_int_distribution<uint64_t> dist( 1, 0xffff'ffff'ffff'ffff );
        for( auto& [e1, e2] : arr ) {
                e1 = dist( mt );
                e2 = dist( mt );
        }

        const auto start = std::chrono::system_clock::now();

        for( auto& [e1, e2] : arr ) {
                std::gcd( e1, e2 );
        }

        const auto finish = std::chrono::system_clock::now();
        std::cout << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() / 1e9 << " sec" << std::endl;
}

コンパイルclang++ -std=c++17 -march=native -O3

環境:Windows subsystem for Linux、clang5.0.1、Intel(R) Core(TM) i5-7200U@2.50GHz(ターボブースト時3.1GHz)

最大公約数を求めるだけの問題に出会うことはそれほどありませんが、アルゴリズムの内部で最大公約数を求める必要があり、ユークリッドの互除法を部品的に使うアルゴリズムはそれなりにあります。そういったアルゴリズム群の(定数倍)高速化のために、gcd関数を高速化することを考えます。

互除法の仕組み

まずは、ユークリッドの互除法で最大公約数(greatest common divisor: GCD)が求まる理由を復習します。以下では、GCD(n, m)と書いたらnmの(数学的事実としての)最大公約数を意味することとし、gcd(n, m)と書いた場合は、(帰納的な定義を持つ)最大公約数を求める関数の定義か、それにnmを代入したものを意味することとします。

まず、GCD(n, m) = GCD(n-m, m)であることが知られています(これの証明は大学入試問題として有名ですが、ここでは割愛します*2)。そのため、gcd(n, m) = if n < m then gcd(m-n, n) else gcd(n-m, m); gcd(0, m) = mといった形で最大公約数を求める関数を、帰納的に定義することができます。これは互減法というアルゴリズムであり、これを使うとコードサイズを小さくできる可能性があります(下記参照。golfer向けであり、実際にプログラムサイズを小さくするのに役に立つとは限りません)。しかし、最悪の入力、つまりgcd(99999,1)などに対して99999回の引き算が発生するため、低速です。

http://golf.shinh.org/reveal.rb?Modular+Inverse/tails_1484205185&c


2019/03/04 01:50 追記:(注意)上記コードは、互減法そのものではありません。モジュラ逆数を求めよという問題に対する解答であり、いわゆる拡張ユークリッドの互除法を使うのが一般的であるところ、互除法ではなく互減法を用いてコードを縮めた例です。

このコードを提出したtailsさんから、混乱を招く恐れがあるとの指摘をいただいたため、追記しました。


GCD(n, m) = GCD(n-m, m)を繰り返し適用すると、GCD(n, m) = GCD(n-m, m) = GCD(n-2m, m) = ... = GCD(mod(n, m), m)となることがわかります。nmは対称なので、n≧mと仮定しても一般性を失いません。この条件のもと、gcd関数は、gcd(n, m) = gcd(m, mod(n, m)); gcd(n, 0) = nのように定義することができます。これがユークリッドの互除法です。

互除法の各ステップの解析

いくつか手で実験をしてみればわかりますが、[1, N]の区間から一様ランダムに取り出した二つの整数n, mに対する互除法の場合、商は1が立ち剰余はn-mになる、というパターンが多いことがわかります。

これの理由は単純で、商として2以上が立つには、片方がもう一方の二倍以上である必要があるからです。そのようになる確率は、二数の選び方が一様ランダムであれば1/2です。 互除法を進めていくと確率は片寄っていきますが、商として1が立つパターンが頻発することは変わりません。また、商が1でなくとも、商は小さいことが多く、大きな商の出現は稀です。一つ一つの数に対して確率が小さいのはもちろんのこと、それらすべての確率の総和も小さいです。

N=264-1として実験してみると、1が立つのは41%、2が立つのは17%、3が立つのは9%、4が立つのは6%、5が立つのは4%、6が立つのは3%、までで八割を占めています。そのあとを追ってみると、13以下が90%、26以下が95%、128以下が99%です。一様ランダムな場合、商はその逆数に比例した確率で出現します。互除法を進めていくにつれ、一様ランダム性が崩れる(区間が減少するのに加え、偶然小さい数が出現する確率が大きくなる)ため解析は複雑になります。しかし、小さい数が出現すると互除法が速やかに終了するため、全体的な確率には大きな影響を及ぼさないようにも思われます。結局のところ、商の出現確率はその逆数に比例していると考えても、大きく外れてはいないでしょう。

剰余演算は遅いので避ける

剰余演算を使うことで単純な引き算と比べて高速化できる、というのが互除法でした。しかし、剰余演算は非常に遅い演算であり、できる限り避けたいところです。

互除法の各ステップにおいて、商として1が立つ場合、数十サイクルもかけて除算を行うのは非効率的であり、1cycleで終わる引き算を行えばよいことになります。 先ほどの解析によれば、商として1が立つことは多い割合で発生するため、数十サイクルのロスを頻繁に行っていることがうかがえます。

商として何が立つかは実際に割り算を行わないといけないため、そのような効率化はできないように思えるかもしれません。

しかし、m != 0の時、`n >= m && m > n - mn / m == 1と等価なため、以下のようにgcd関数を作ることができます。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        if( m == 0 ) { return n; }
        uint64_t t = n - m;
        if( m > t ) {
                return gcd_impl( m, t );
        } else {
                return gcd_impl( m, n % m );
        }
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

しかし、これを使っても速度が数%改善される程度にとどまります。これは、商として1が出るパターンの出現頻度は四割であり、かつそれは予測不可能に出現するため、分岐予測が全く当たらないのが原因です。

当てられない分岐は遅いので避ける

高速化のテクニックに、当てられない分岐を取り除くというものがあります。

m == 0も予測不可能に出現しますが、これなしというわけにはいかないこと、出現は最後に一回だけであり影響は軽微であることから、この分岐を取り除くことは考えません。 n > mは、入力の最初に50%の確率で発生し、これは絶対に当てられません。

よって、m > tの分岐を取り除くことを考えます。その前にちょっと準備をします。

予測できない分岐なしに互減法を実装する

予測できない分岐を取り除くため、x86の命令セットには、cmov (conditional move)命令が用意されています。これは分岐なしにcond ? a : b;のような計算をすることができる命令です。昔のIntelはさぼっていたのかこの命令のレイテンシが2cycleでしたが、最近では1cycleとなっているため、さらに気軽に使えるようになりました。cmov命令を使う欠点は、分岐すれば計算する必要がなかったかもしれないのに、abも計算しないといけなくなる点です。

互減法は以下のように書けます。

while( m != 0 ) {
        uint64_t t = n - m;
        bool q = m > t;
        n = q ? m : t;
        m = q ? t : m;
}

このコードを書いたとしても、必ずしもcmov命令が使われるとは限りません。実際g++でコンパイルすると条件分岐を使ったコードにコンパイルされてしまいます。一方clang++でコンパイルするとcmov命令が使われました。

互除法と互減法を組み合わせる

先ほどの実験では、商として1が現れることが多く、特に128以下の商が立つ状況は99%の確率で発生することを明らかにしました。そこで、互除法ではなく互減法の方が速いのではないか、と思いたくなりますが、これは誤りです。一様ランダムに選ばれた符号なし64bit整数同士の互除法は実験してみると平均37ステップかかることがわかります。つまり、三回に一回は128を超える商がどこかのステップで出現することとなり、そのような場合に引き算を続けるのは非常に非効率的です。

一方、多くの場合に商が1であり、そのような場合は除算より引き算の方が一桁速いことを考えると、互減法が優れている局面は確かに存在することもわかります。

よって、うまく互除法と互減法を切り替える、または組み合わせるなどの手法が有効そうであることがわかります。

ここで、切り替えるために条件分岐を使ってしまうのは、やはり当てられない分岐となるためによくなさそうです。実際、先ほどの失敗は当てられない分岐を作ってしまったことによるものでした。

分岐せずにうまく切り替える方法はあるのでしょうか?

キーとなるアイデアは、「分岐予測に失敗して時間を浪費するくらいなら、その時間に無駄になるかもしれない仕事をしておけば、その仕事が無駄ではなかった場合に高速化する」というものです。このアイデアニュートン・ラフソン法を用いている平方根関数を高速化するときにも使いました。これは、「収束判定という予測できない分岐をするかわりに、どんな入力に対しても収束する回数だけ回す」というものです。すでに収束している場合は無駄な加速を行っていることになりますが、それでも分岐予測に失敗するよりは速いのです。

このアイデアを使うと、以下のようなコードになります。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = m > t;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

ここで、ループの回数10回は適当であり、最適化の余地があります。このコードは、互減法を10回行った後、互除法的な除算を一回やる、というのを繰り返すものになっています。商が小さければ互減法の部分でnmが小さくなってくれるため、トータルの除算回数が減ることになり、高速化します。一方、商が大きくなる場合でも、10回しか引き算せずに除算を行うようになっているため、延々と引き算をし続ける、といった非効率的な状況にはまることを避けられます。

なお、10回のループに関してはループアンローリングが行われることを期待しています。もしそうでなくても、少なくとも当てられる分岐にはなっています。よって、この関数の中には、if( m == 0 ) { return n; }という避けようのない分岐を除いて、予測不可能な分岐はありません。

先ほどのアイデアがどう使われているかを別の言葉でも説明してみます。

もし商が大きい場合、10回の引き算しても意味はなく、なるべく早く除算モードに切り替えることが高速化につながります。しかし、商が大きいという予測不可能な分岐をしてモードを切り替えるくらいなら、10回程度の引き算という無駄な仕事をした方が高速です。

逆に(運よく)商が小さい場合、引き算は無駄な仕事とならず、高速化に貢献します。

まとめ

最大公約数を求める有名なアルゴリズムユークリッドの互除法には除算が含まれており、一般的なCPUで計算すると遅くなりがちです。互減法を部分的に組み合わせて高速化を図ろうとしましたが、予測不可能な分岐が追加されたため、それほど高速化はしませんでした。そこで、無駄かもしれない仕事と引き換えに予測不可能な分岐を除去する方法を考えました。その手法を使うと、元のナイーブな互除法に比べて1.5倍ほど高速に動作するプログラムができました。

*1:実際の比例係数は、黄金比phiを底とする対数を用いてlog_phi(10)とあらわせて、これは4.8ほどです。

*2:というより忘れて思い出せません……。

完全精度 expf 関数の作り方

C言語では、expf関数などの超越関数の精度は一切保証されないことになっています*1。 そうは言っても、数学的な結果に近い結果が返ってくることが、常識的に信じられています。

数学的な結果に最も近いfloatdoubleが返ってくることを、完全精度(丸め誤差0.5ulp保証*2)といいます。また、数学的な結果に最も近い、または二番目に近い値が返ってくることを、丸め誤差1ulp保証といいます。

一般に、超越関数の結果は、丸め境界に非常に近い値になり得るため、完全精度を実現するのは困難です(テーブルメーカーのジレンマ)。一方、1ulp保証ならば現実的な計算量で実現可能です。

ちょっと調べてみたところ、手元でビルドしたg++8.2では、double std::exp(double)は完全精度っぽいですが、float std::exp(float)は単にdouble std::exp(double)関数の結果を丸めた結果*3を返しているような挙動を示しており、完全精度になっていません。 そこで、自分で作ってみます。

完全精度expf関数の作り方

ここでは、高速化については考えず、とにかく完全精度を追求するものとします。また、テーブルに頼ることはできるだけ避けます。

正しい方針

実は、expfであれば完全精度を達成することは意外と簡単です。

最も丸め境界に近いのはexpf(-0x1.d2259ap+3) = 0x1.fa6636p-22 であり、真の値は0x1.fa66350000001p-22付近にあります。この場合であっても、double精度で正確に計算し最後に一回だけ丸めればぎりぎり有効数字が足ります。

しかし、そのためには簡易的なdouble-double演算*4が必要になってくるため、やや計算量が高くつきます。

どのくらい雑にやってよいのか

実際には丸め境界に近い値はそれほど多くありません。double精度の下2bitを正確に求める必要があるのは、上記の例に加えて

  • expf(-0x1.e1dbe2p-8) = 0x1.fc3fd2p-1 真の値は0x1.fc3fd10000002p-1付近
  • expf(-0x1.c1c4b8p-10) = 0x1.ff1f4ep-1 真の値は0x1.ff1f4effffffcp-1付近
  • expf(0x1.fdff02p-17) = 0x1.000100p+0 真の値は0x1.0000ff0000003p+0付近
  • expf(0x1.cd3982p-14) = 0x1.000734p+0 真の値は0x1.000734ffffffcp+0付近
  • expf(0x1.8d7cb6p-12) = 0x1.0018dap+0 真の値は0x1.0018d90000004p+0付近

くらいしか存在しないため、運が悪くなければ慎重にdoubleで計算するだけで丸め誤差の問題を発生させずに済むかもしれません。

形式的な証拠とはなりえませんが、定義域のすべての値について正しく丸められていることを確認すれば、完全精度を達成していること自体は確かめられます。 expfの定義域は-104.0fから88.0fほどであり、22.4億通り(231よりちょっと多いくらい)で正しく丸められているかを確認すれば目的を達成できます。 これは、最近の*5パソコンであれば、特に工夫をしなくても1時間程度で終わります。

実際にとった方針

基本的にはテーラー展開を使って求めることになりますが、引数の値が大きい場合、収束に必要な項数が多くなります。これは誤差の増大も意味するため、避けるべきです。

方針としては、 \exp(s+t) = \exp(s) \times \exp(t) をうまく活用していきます。 まず、 \exp(s)が簡単に計算できるように、かつtの絶対値が小さくなるように、うまくsを選びます。 こうすると、 \exp(t)テーラー展開で求めるときに必要な項数が少なくなります。

 \exp(s) の計算法

 \exp(s) を簡単に計算する手法として、メジャーなのは以下の二つです。

  • s \ln(2) の整数倍になるようにとる。すると単に二のべき乗となり、std::ldexp関数で効率よく計算できる。tの絶対値は0.5以下となる。
  •  \exp(s) のテーブルを持つ。tの絶対値を小さくすることができるため、上の方法より収束が速い。収束の速さはテーブルサイズとのトレードオフとなる。

しかし、今回は難しいことを考えず、sとして整数をとりました。 \exp(s)は繰り返し二乗法により求めました。これによりいくらか誤差が乗りますが、最終的な結果への影響はないようです。

 \exp(t) の計算法

tの絶対値が0.5以下の時、| t15 / 15! | <  2^{-55} であるため、double精度を出すためには14次までのテーラー近似を利用すれば十分です。

なお、テーラー展開を13次までとした場合は最終的な結果に誤差が発生しました。

テーラー近似は展開中心点から離れるほど近似度が悪くなることが知られています。この問題を解決する方法として、minimax近似というものあります。 minimax近似は、特定の区間内での最大誤差が最も小さくなる近似多項式です。しかしその係数は複雑になります。

今回はややこしいことを考えず、単純なテーラー近似を利用しました。

ここで、テーラー近似を求める段階で注意が必要です。

一般に多項式の値を求める場合、乗算回数が少なくなる*6ホーナー法が用いられます。

ホーナー法とは、  a_0 + x \times ( a_1 + x \times ( a_2 + x \times a_3 ) ) のように計算する多項式計算法のことです。

この方法は、係数の誤差が蓄積する性質があるため、限られた精度で正確な値を求めるのに向いていないようです(むしろ乗算回数が少なくなることから、多倍長演算向けだと思われます)。

今回は愚直に、 \sum_{i=0}^{14} c_i t^{i} を求めることにしました。ここで、絶対値の小さいほうから足していくことで計算誤差を最小化します。

実際のコード

constexpr std::array<double, 15> make_exp_coeff() {
        std::array<double, 15> arr {};
        double fact = 1.0;
        for( int i = 1; i < 15; ++i ) {
                fact *= i;
                arr[i] = 1. / fact;
        }
        return arr;
}

double power( double x, unsigned int n ) {
        double ret = 1.0;
        for( ; n; n >>= 1 ) {
                if( n & 1 ) { ret *= x; }
                x *= x;
        }
        return ret;
}

// |t| <= 0.5  --> | t^15/15! | < 2^-55
double exp_taylor14( double t1 ) {
        constexpr auto coeff = make_exp_coeff();
        const double t2 = t1 * t1;
        const double t4 = t2 * t2;
        const double t8 = t4 * t4;
        double ret = 0.0;
        ret += t8 * t4 * t2      * coeff[14];
        ret += t8 * t4      * t1 * coeff[13];
        ret += t8 * t4           * coeff[12];
        ret += t8      * t2 * t1 * coeff[11];
        ret += t8      * t2      * coeff[10];
        ret += t8           * t1 * coeff[9];
        ret += t8                * coeff[8];
        ret +=      t4 * t2 * t1 * coeff[7];
        ret +=      t4 * t2      * coeff[6];
        ret +=      t4      * t1 * coeff[5];
        ret +=      t4           * coeff[4];
        ret +=           t2 * t1 * coeff[3];
        ret +=           t2      * coeff[2];
        ret +=                t1 * coeff[1];
        return ret + 1.0;
}

float myExp( float x ) {
        constexpr double exp_plus_one  = 0x1.5bf0a8b145769p+1;
        constexpr double exp_minus_one = 0x1.78b56362cef38p-2;
        const int s = std::round( x ); // round to nearest(ties to max magnitude)
        const double t = x - s; // |t| <= 0.5

        const double exp_s = s > 0 ? power( exp_plus_one, s ) : power( exp_minus_one, -s );
        const double exp_t = exp_taylor14( t );
        return static_cast<float>( exp_s * exp_t );
}

検証コード

#include <iostream>
#include <cstdint>
#include <cstring>
#include <cassert>

#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

using idd = kv::interval<kv::dd>;

template<class To, class From>
To bit_cast(const From& from) noexcept {
        To to;
        static_assert( sizeof to == sizeof from );
        std::memcpy(&to, &from, sizeof to);
        return to;
}


float next( float x ) {
        assert( !std::signbit(x) );
        return bit_cast<float>( bit_cast<uint32_t>(x) + 1 );
}

float before( float x ) {
        assert( !std::signbit(x) );
        if( x == 0 ) { return -0.0f; }
        return bit_cast<float>( bit_cast<uint32_t>(x) - 1 );
}

idd intervalFloat( float x ) {
        double y = next(x);
        double z = before(x);
        return idd( (x+z)/2, (x+y)/2 );
}

float nearestFloat( idd ia ) {
        // 二重丸め
        float x = (float)(double)mid(ia);
        float y = next(x);
        float z = before(x);
        idd ix = intervalFloat(x);
        idd iy = intervalFloat(y);
        idd iz = intervalFloat(z);
        if( subset(ia, ix) ) { return x; }
        if( subset(ia, iy) ) { return y; }
        if( subset(ia, iz) ) { return z; }
        std::cerr.precision(50);
        std::cerr << x << std::endl;
        std::cerr << y << std::endl;
        std::cerr << z << std::endl;
        std::cerr << ia.lower() << std::endl;
        std::cerr << ia.upper() << std::endl;
        throw nullptr;
}

#include "my_exact_expf.hpp"

void check( float x ) {
        idd ix(x);
        auto exact_exp_range = exp(ix);
        float nearest = nearestFloat(exact_exp_range);
        float my = myExp(x);

        if( nearest != my ) {
                std::cerr << "x     " << x << std::endl;
                std::cerr << "x     " << std::hexfloat << x << std::defaultfloat << std::endl;
                std::cerr << "my    " << my << std::endl;
                std::cerr << "exact " << exact_exp_range << std::endl;
                std::cerr << abs(mid(exact_exp_range) - my) / (nearest - before(nearest)) << " ulp" << std::endl;
        }
}

int main() {
        std::cerr.precision(50);
        const uint32_t minimum = bit_cast<uint32_t>(-104.0f);
        const uint32_t maximum = bit_cast<uint32_t>(88.0f);
        for( uint32_t i = bit_cast<uint32_t>(-0.0f); i < minimum; ++i ) {
                check( bit_cast<float>(i) );
        }
        for( uint32_t i = bit_cast<uint32_t>(0.0f); i < maximum; ++i ) {
                check( bit_cast<float>(i) );
        }
        return 0;
}

この検証コードのコンパイルには、区間演算ライブラリkvが必要です。

github.com

なお、そのままコンパイルするにはboostが必要ですが、このライブラリはC++03対応を維持するためにC++11で導入された標準関数を使っていないようです。 そのため、以下の強引な置換により、boostへの依存は解消されます。

sed -i s/boost/std/g
sed -i s/enable_if_c/enable_if/g

コンパイルオプションは、高速化のため、g++-8 -std=c++17 -march=native -O3 -DKV_FASTROUND -DKV_USE_TPFMA -I ../../kv/を指定しました。使ったCPUはIntel(R) Core(TM) i5-7200U CPU @ 2.50GHzです。かかった時間はほぼ一時間ぴったりでした。

区間演算によりexp(x)の真値を含む区間が求まります。正しくfloatに丸めた場合の可能性が二つ以上ある場合はthrow nullptrされるはずですが、そのようなことは起こりませんでした。 また、(float)(double)mid(ia)の部分はdoubleに丸めてからfloatに丸めているので正しく丸められていない可能性がありますが、依然1ulp保証ではあるので前後のfloat値を確認しています。 本当は正しく求める方法もあるのですが、手を抜きました。

このプログラムが何も出力せずに終了したことにより、実装したmyExp関数は完全精度であることを確かめました。0.0付近の非正規化数であっても0.5ulp保証されています。

ただし、inf付近でどちらに丸めるべきかについては、よくわからなかったので、その部分について0.5ulp(??)保証されているかはわかりません。

g++のよくわからない挙動

g++では、なぜかstd::expf関数が使えません。::expf関数はあるらしいです。::expf関数はおしいですが完全精度になっていません(この記事で書いたコードが意味のないものとならなくてよかったです)。また、type detectionイディオムを使った確認により、std::exp(float)オーバーロード自体は存在し、返り値の型はfloatであることを確認しました。

#include <cmath>

template<class T>
struct TD;

int main() {
  TD<decltype(std::exp(1.0f))> td;
}

なぜstd::exp(double)::expf(float)は完全精度であるのにstd::exp(float)だけ精度が低いのでしょうか?

::expf 関数の精度

以下の98個の入力について0.5ulpを超える誤差があることを確認しました。最大の誤差は0.5000031997 ulp 未満でした。

  • -0x1.c1c4b8p-10 (上で紹介した入力のうちの一つ)
  • -0x1.abb574p-9
  • -0x1.f5723ep-9
  • -0x1.f9c67ep-9
  • -0x1.fd2ddcp-9
  • -0x1.021d74p-8
  • -0x1.2e4ab4p-8
  • -0x1.400432p-8
  • -0x1.4b0f34p-8
  • -0x1.4be042p-8
  • -0x1.4d8260p-8
  • -0x1.4dcdc2p-8
  • -0x1.58ef7ap-8
  • -0x1.db8546p-7
  • -0x1.f691dap-7
  • -0x1.f99ab6p-7
  • -0x1.02f812p-6
  • -0x1.074b54p-6 (最大誤差となる入力。真値は0x1.f7d67affff94ap-1付近ですが、0x1.f7d67cp-1が返ってきます)
  • -0x1.a59e4ep-6
  • -0x1.b17ee0p-6
  • -0x1.b1db4ap-6
  • -0x1.8e02c8p-5
  • -0x1.cd9502p-4
  • -0x1.d0a5b6p-4
  • -0x1.fdcbf6p-4
  • -0x1.2b6a50p-3
  • -0x1.0faa34p-2
  • -0x1.19662ep-2
  • -0x1.d66d74p-2
  • -0x1.074b34p-1
  • -0x1.c2aaeap-1
  • -0x1.df1c00p-1
  • -0x1.efdceap-1
  • -0x1.08abe2p+0
  • -0x1.74a396p+0
  • -0x1.7f4296p+0
  • -0x1.b75242p+0
  • -0x1.ff2602p+0
  • -0x1.7dd9aap+1
  • -0x1.d2445ap+1
  • -0x1.f66118p+1
  • -0x1.2b0b04p+2
  • -0x1.90ced0p+2
  • -0x1.bc933ep+2
  • -0x1.e4c6c0p+2
  • -0x1.e8ebf0p+2
  • -0x1.1089b0p+3
  • -0x1.287d08p+3
  • -0x1.93e788p+3
  • -0x1.d2259ap+3 (上で紹介した、最も丸め境界に近くなる入力)
  • -0x1.de6f7cp+3
  • -0x1.1f8dc4p+4
  • -0x1.433aaep+4
  • -0x1.b43edcp+4
  • -0x1.c95efap+4
  • -0x1.23134ap+5
  • -0x1.444328p+5
  • -0x1.2219dep+6
  • -0x1.2c4726p+6
  • -0x1.5ce26ap+6
  • 0x1.6438b6p-8
  • 0x1.9785f6p-8
  • 0x1.1c2c32p-6
  • 0x1.2dc3f0p-6
  • 0x1.2e3554p-6
  • 0x1.c649eap-6
  • 0x1.3a7784p-5
  • 0x1.980654p-5
  • 0x1.eaddb4p-5
  • 0x1.eff8cep-5
  • 0x1.4feb94p-4
  • 0x1.599d64p-3
  • 0x1.72793ap-3
  • 0x1.b28da4p-3
  • 0x1.f3518cp-3
  • 0x1.54f076p-2
  • 0x1.0d5f92p-1
  • 0x1.39a300p-1
  • 0x1.70f4e4p-1
  • 0x1.b94234p-1
  • 0x1.bf4b50p-1
  • 0x1.13f37ap+0
  • 0x1.c8264ap+0
  • 0x1.4c1e18p+1
  • 0x1.69512ep+1
  • 0x1.172096p+2
  • 0x1.24b172p+3
  • 0x1.3983a8p+3
  • 0x1.56f768p+3
  • 0x1.ddc228p+3
  • 0x1.5d4bd4p+4
  • 0x1.6255f4p+5
  • 0x1.91fceap+5
  • 0x1.d8afbep+5
  • 0x1.097f6ap+6
  • 0x1.15298ap+6
  • 0x1.1ddcb8p+6
  • 0x1.3829eep+6

*1:超越関数だけでなく、sqrt関数以外の無理関数は精度が保証されません。

*2:ulpはunit in the last placeの略で、最後の桁の重みのことです。数学的な結果に最も近い値に丸めた場合、その丸め誤差は0.5ulpになります。

*3:これでは二回丸めが発生し、正しくありません。なぜ二回丸めると正しくないかというと、0.495→0.5→1.0のような間違いが発生するためです。

*4:doubleを二個使う多倍長表現のことをdouble-doubleと呼びます。疑似四倍精度とも呼ばれます。

*5:FMA命令に対応しているHaswell以降のCPUでは、double-double演算が高速化します。FMA命令は計算結果を一回だけ丸めるため、演算誤差の補正が簡単になるためです。

*6:必ずしも乗算回数が最小となるわけではありません。さらに乗算回数が少ない方法として、Knuthのadaptation of coefficientsという方法があるようです。

C++におけるinline指定子の誤解

C++にはinline指定子があります。しかし、このinline指定子は、その名前から誤解されがちです。

(誤解)inline指定子を付けると、その関数はインライン展開される

これは誤解であり、正しくありません。正しくは、以下のような意味になります。

  • inline指定子を付けた場合、複数回の定義が許される。また、複数の翻訳単位*1で定義されたとしても実体は一つになることが保証される*2
  • (C++14まで)inline指定子を付けた場合、コンパイラはその関数をインライン展開するヒントとすることができる。

「inline指定子を付けた場合、コンパイラはその関数をインライン展開するヒントとすることができる。」という記述はC++17以降で無くなりました。 ただし、コンパイラの最適化は見える結果が変わらない範囲で自由に行えるため、inline指定子をヒントにすることは問題ありません。

ちなみに、テンプレート関数、constexpr関数、クラス定義内で定義される関数、はいずれも暗黙にinline指定されたものとみなされます。default定義されたメンバ関数も暗黙にinline指定されたものとみなされます。

ところで、inline指定されている場合(これには上述の暗黙指定の場合を含みます)、複数回の定義が許されますが、その場合でも全ての定義は字面上同一かつ意味上同一でないといけません(one definition rule, ODR)。

外部リンケージを持つ(static指定されておらず、無名名前空間に囲まれていない)関数がODRを守れていない場合、コンパイラは特に指摘してくれません(診断不要:no diagnostics required)。このようなことが発生すると、リンクの順番によってプログラムの実行結果が変わるなどの意味不明なバグにつながり、非常にデバッグが困難です。つまり、うっかり同名の関数や変数を書いてしまってもコンパイラが指摘してくれず、バグが埋め込まれるということになります。そのため、特にグローバル名前空間inlineなものを置くことは忌避されるべきです。

うっかりを防ぐためには、static inlineとすべきです。ただし、実行効率まで考えると常に最善とは限りません(後述します)。

(誤解)inline指定子を付けると、その関数はインライン展開されやすくなる

先ほどの誤解よりは少しましになりましたが、依然誤解です。「コンパイラはその関数をインライン展開するヒントとすることができる」をうがって読むと、「コンパイラはその関数をインライン展開 しない ことの根拠としても使える」ことがわかります。そんな例は実在するのかという感じですが、以下に実例を示します。

inc.h

#include <array>
#include <cstdint>

#define M1(i) s += (i)*(i), arr[s%SIZE] = i;
#define M2(i) M1(i) M1(i+1)
#define M4(i) M2(i) M2(i+2)
#define M8(i) M4(i) M4(i+4)
#define M16(i) M8(i) M8(i+8)
#define M32(i) M16(i) M16(i+16)
#define M64(i) M32(i) M32(i+32)
#define M128(i) M64(i) M64(i+64)


DECL_SPEC std::array<std::size_t, 128> LargeFunction(std::size_t s) {
    std::array<std::size_t, 128> arr {};
    for( int t = 0; t < 128; ++t ) {
        M128(0);
    }
    return arr;
}

source1.cpp

#include "inc.h"

std::array<std::size_t, 128> callLargeFunction(std::size_t s) { return LargeFunction(s); }

source2.cpp

#include "inc.h"

std::array<std::size_t, 128> anotherCallLargeFunction(std::size_t s) { return LargeFunction(s); }

main.cpp

#include <iostream>
#include <array>
#include <cstdint>

std::array<std::size_t, 128> callLargeFunction(std::size_t s);
std::array<std::size_t, 128> anotherCallLargeFunction(std::size_t s);

int main() {
    std::cout << callLargeFunction(1)[3] << std::endl;
    std::cout << anotherCallLargeFunction(2)[4] << std::endl;
}

以下のコンパイルオプションで分割コンパイルし、リンクした場合、バイナリのサイズは13120バイトとなります。コンパイルオプションとして-Osの代わりに-O1-O2-O3を使う場合、17176バイトとなります。-fltoオプションを付ける場合、四つの最適化オプション-O1-O2-O3-Osのいずれでも12936バイトとなります。

g++-8 -std=c++17 -Os -DDECL_SPEC=inline source1.cpp -c -o source1_inline.o
g++-8 -std=c++17 -Os -DDECL_SPEC=inline source2.cpp -c -o source2_inline.o
g++-8 -std=c++17 -Os main.cpp -c -o main_inline.o

g++-8  source1_inline.o source2_inline.o main_inline.o -o inline

一方、-DDECL_SPEC=inlineの代わりに-DDECL_SPEC=static-DDECL_SPEC='static inline'を指定した場合、最適化オプションによらず17176バイトとなります。また、-fltoオプションを付けた場合、-O1の場合のみ12896バイト、その他の三つの最適化オプション-O2-O3-Osのいずれでも12936バイトとなります。

つまり、 staticstatic inlineの代わりにinlineとすることでインライン展開を抑制することができる 場合がある(具体的には-Os指定かつ-flto指定なし)ということです。

なぜこのようなことが起こったのでしょうか?

-fltoオプションを付けた場合、他の翻訳単位で使われるか否かを見てからインライン展開をするかどうかを判断することができるようにコンパイルすることになります。よって、基本的にはインライン展開せずにコンパイルし、リンク時にインライン展開を行うか判断することになります。そのため約13KBの小さなバイナリとなります。

一方、-fltoオプションを付けない場合、他の翻訳単位で使われるか否かはわからないままインライン展開するべきかの判断を下す必要があります。インライン展開をすべきと判断した場合、二つの関数にインライン展開されるためバイナリサイズが増大し、約17KBのバイナリとなります。

ここで、static指定子がついている場合(-DDECL_SPEC=static-DDECL_SPEC='static inline'とした場合)を考えてみると、他の翻訳単位からその関数が呼ばれることはない*3ため、実体を作る必要はありません。そのため、その翻訳単位で一回しか呼ばれないのならインライン展開してしまうのが最善となります。

一方、単独でinlineとしている場合、「複数の翻訳単位で定義されたとしても実体は一つになることが保証される」ことから、インライン展開しないで実体を作っておけば、他の翻訳単位で作った実体とリンク時にマージできる可能性があり、バイナリサイズを小さくする役に立てることができます。そのため、むしろインライン展開を抑制する働きを持つことがあります。

まとめ

inline指定子はもはやコンパイラにインライン展開を指示するものではありません。また、インライン展開を促進する効果があるとも限らず、むしろインライン展開を抑制する働きを持つことさえあります。この効果はstatic inlineとした場合には働かないため、ODR違反に細心の注意を払いつつstatic inlineとしないことが有効な局面もあり得るということを示しました。

補遺

inlineはむしろインライン展開を抑制する働きを持つことさえある、と書きましたが、これはあくまで特殊例です。現にコンパイラinlineキーワードがついている関数についてインライン展開閾値を下げるということを行っているようです。

また、関数のマージはinlineを付けない限り発生しないわけではなく、コンパイラの最適化によって勝手に行われることがあるようです。しかも、関数ポインタをとっているのに勝手にマージされて同じアドレスになって困るなどの現象も発生することがあったらしいです(VC9で発生した例→RTTI を使わずに Boost.Any を実装する方法 - melpon日記 - HaskellもC++もまともに扱えないへたれのページ)。

*1:難しい用語ですが、分割コンパイルするときの一つのcppファイルだと理解すれば大体あっています。

*2:平たく言うと、ヘッダファイルで定義してソースファイルにインクルードする場合、inlineを付けないと実体が複数になって困りますが、inlineを付ければそういうことが起こらないということです。つまり、ヘッダファイルに定義を書くことができるようになるという利点があります。

*3:"その関数"自体が呼ばれることはあり得ますが、staticがついている以上、実体としては別の関数として扱われます。

いろいろなボードゲームの複雑性クラスとその直感的説明

いろいろなボードゲームの複雑性クラスは語られることはありますが、多くの場合直観的説明がなく、論文を読むしかない状況にありました。調べたことをまとめておきます。

チューリングマシンの機能

決定性チューリングマシン

普通のコンピュータのことだと思えばよいです。

非決定性チューリングマシン

決定性チューリングマシンに、『ひとり用のゲーム*1を遊んでいるとき、「ここはどっちに進むとクリアできるか?」という質問に答えてくれる機能』が追加されたものだと思えばよいです。これは何回でも使えます。 「なぜそうなのか?」という質問には答えてくれないので、本当にクリアできることを確かめるには自分でやってみるしかありません。

この機能は使わなくてもよいので、決定性チューリングマシンより強力です。

交代性チューリングマシン

決定性チューリングマシンに、『ふたり用のゲーム*2を遊んでいるとき、「この局面は必勝局面か?」という質問に答えてくれる機能』が追加されたものだと思えばよいです。 あるいは、『絶対間違わないプロ棋士次の一手を聞ける権利』と思っても同じです。これも何回でも使えます。 やはり「なぜそうなのか?」という質問には答えてくれないので、本当に必勝であることを確かめるのは難しいです。

この機能は、ひとり用のゲームを解くのにも使えるため、非決定性チューリングマシンより強力です。

複雑性クラス

多項式時間

P

決定性チューリングマシンを用いると多項式時間で解ける決定問題のクラスです。 ようするに、普通のコンピュータで多項式時間で解ける、Yes, Noを答える問題です。

NP

非決定性チューリングマシンを用いると多項式時間で解ける決定問題のクラスです。 "証拠"を見せられた時に多項式時間で納得できるものであるとも解釈されます。この"証拠"としては、"ここはどっちに進むとクリアできるか?"という質問への答えの列が挙げられます。

当然Pを含んでいますが、真に包含しているかは不明です(P≠NP予想)。

AP

交代性チューリングマシンを用いると多項式時間で解ける決定問題のクラスです。 交代性チューリングマシンは強すぎてなかなか直感的説明が難しいです。

当然NPを含んでいますが、真に包含しているかは不明です。

多項式空間

PSPACE

決定性チューリングマシン多項式量のメモリを与えると解けるけ決定問題のクラスです。 多項式時間では自明に多項式空間しか使用できない(時間計算量と空間計算量の関係)ため、P、NP、APを全て包含します。

実はAP=PSPACEであることがわかっています(並列計算原理)。P, NPを真に包含するかは不明です。

NPSPCAE

非決定性チューリングマシン多項式量のメモリを与えると解ける決定問題のクラスです。

実は、非決定性チューリングマシン多項式空間を使って解ける問題は、その二乗程度の空間を使ってよければ決定性チューリングマシンでも解けることが知られています(サヴィッチの定理)。 多項式を二乗しても多項式であるため、PSPACE=NPSPACEです*3

APSPACE

交代性チューリングマシン多項式量のメモリを与えると解ける決定問題のクラスです。 当然NPSPACEを含んでいますが、真に包含しているかは不明です。

指数時間

EXPTIME

決定性チューリングマシンを用いると指数時間で解ける決定問題のクラスです。 PSPACEやAPSPACEを含んでいます。なぜなら、多項式量のメモリで表される状態数は高々指数個しかないため、それを全探索するのには指数時間しかかからないからです(空間計算量と時間計算量の関係)。

実は、APSPACE=EXPTIMEであることがわかっています(並列計算原理)。 また、Pを真に包含することもわかっています(時間階層定理)。 NP, PSPACEを真に包含しているかは不明です。

NEXPTIME

非決定性チューリングマシンを用いると指数時間で解ける決定問題です。

当然EXPTIMEを含んでいますが、真に包含しているかは不明です。

AEXPTIME

交代性チューリングマシンを用いると指数時間で解ける決定問題です。

当然NEXPTIMEを含んでいますが、真に包含しているかは不明です。

指数空間

EXPSPACE

決定性チューリングマシンに指数量のメモリを与えると解ける決定問題です。

PSPACEの時の議論と全く同じ議論により、NEXPTIMEを包含すること(時間計算量と空間計算量の関係)、NEXPSPACEと同じであること(サヴィッチの定理)、AEXPTIMEと同じであること(並列計算原理)がわかっています。 また、PSPACEを真に包含することがわかっています(空間階層定理)。

いろいろなボードゲーム

スリザーリンク*4をはじめとする多くのペンシルパズル、充足可能性問題

その問題に解があるかを判定する問題はNP完全であることがわかっています。

  • 指数的な分岐があるため、多項式時間では解けなさそう
    • Pではなさそう
  • 答えを見れば「確かにこの問題の答えである」と多項式時間で納得できる(唯一解であるかはまた別の問題)
  • 非決定性チューリングマシンの機能に『ここに線を引くべきか?』などを聞きながら解いていくことができる
    • NPに含まれる

オセロ*5をはじめとする単調に終局に近づく二人零和完全情報ゲーム

ボードサイズを一般化し、(初期配置から到達できるとは限らない)特定の局面に必勝手順があるかを判定する問題はPSPACE完全であることがわかっています。 語彙の集合が与えられる"しりとり"(Generalized geography)も全く同様です。

  • 指数的な分岐があるため、多項式時間では解けなさそう
    • Pではなさそう
  • 最善手順のみ見せられても多項式時間では納得できない(他の手を打ったらどうなるの?となってしまう)
    • NPではなさそう
  • 多項式手数で終わるため、深さ優先探索をすれば多項式空間で解くことができる
    • PSPACEに含まれる
  • 多項式手数で終わるため、交代性チューリングマシンの機能を用いて、『絶対間違わないプロ棋士』同士を戦わせれば多項式時間で解くことができる
    • APに含まれる

倉庫番をはじめとする可逆なギミックがある一人用コンピュータゲーム

その問題に解があるかを判定する問題はPSPACE完全であることがわかっています。 マリオブラザーズ*6などもここに分類されます。

  • 指数的な分岐があるため、多項式時間では解けなさそう
    • Pではなさそう
  • 最短手順が指数的に長い問題もあるため、それを見せられても多項式時間では納得できなさそう
    • NPではなさそう
  • 盤面を覚えておくのは多項式空間でできるため、非決定性チューリングマシンの機能に『次の一手』を聞きながら手順を進めていけば多項式空間で解くことができる
    • NPSPACEに含まれる

ところで、ちょっと考えただけでは

  • 最短手順が指数的に長い問題もあるため、反復深化深さ優先探索を行ってもコールスタックのために指数的な空間が必要そう
    • PSPACEでなさそう……?

となってしまいそうですが、PSPACE=NPSPACEなので、これはおかしいです。実は、時間を犠牲にしてメモリ使用量を削るテクニック(これこそがサヴィッチの定理です)が存在し、それを使うとPSPACEに収まることがわかります。 反復深化深さ優先探索の時間計算量は、Cを適当な定数、Nを入力長として(たぶん)Ο(CN)ですが、メモリ使用量削減手法では最悪時Ο(CN2)になるはずです。

このテクニックは一人用ゲームにしか適用できない手法です。二人用ゲームの深さ優先探索に適用する場合、特定の手を打つと必勝手順に持ち込まれてしまうか、という情報をコールスタックに記録する必要がある点が異なります。一人用ゲームであれば、全探索したことさえわかればよいので、コールスタックを愚直に記録する必要はない、という発想によりメモリ使用量を削ることができます。

時間を犠牲にしてメモリ使用量を削るテクニック(サヴィッチの定理)の適用

具体的なテクニックは驚くほど簡単です。盤面は指数通り(例えば2T通り、ここでTはNの多項式)しかなく、これを枚挙することは多項式空間でできます。また、その盤面が解局面であるかは簡単に判定できます。そこで、あり得る解局面全てを枚挙することができるので、それらすべてについて、「初期盤面からその解局面に到達できるか(到達する手順が存在するか)?」という部分問題を解けば元問題が解けることになります。他の部分問題の解はその部分問題の解に影響を及ぼさず、どの部分問題まで解いたかは多項式空間(T bit)で覚えておくことができます。

ここで、「局面Aから局面Bにk手で到達できるか?」という問題を解くには、「局面Aから局面Cにk/2手で到達できる、かつ、局面Cから局面Bにk/2手で到達できる」かをあり得るすべての局面Cについて調べれば十分です。 これで問題サイズが半分になりました。以下再帰的に問題サイズを小さくしていけば、この再帰はlog k段で終了します。kが1である基底ケースは自明に解けます。

全く同じ局面は最短手順の中に出てこないはずなので、その手順の長さは盤面の通り数2Tで抑えられます。よって、「初期盤面からこの解局面に到達できるか?」という問題は「初期盤面からこの解局面に2T手以内で到達できるか?」という問題と同値になります。0から2T手までの各手数について、先の手法を使うと再帰の段数は高々T段になります。また、各再帰のコールスタックでは、局面Cとしてどこまで検索したかというやはり多項式に収まる情報(T bit)を持っておけば十分です。合わせてT2程度の空間計算量で解けることになります。

この手法はメモリ使用量は削減できますが、同じ仕事を非常に多数回行うため非常に非効率的であり、実用的な手法ではありません。

詰将棋*7、日本ルールの囲碁*8など

ボードサイズを一般化し、(初期配置から到達できるとは限らない、任意の位置に任意の駒を配置できる)特定の局面に必勝手順があるかを判定する問題はEXPTIME完全であることがわかっています。

  • 指数的な分岐があり、多項式時間では解けなさそう
    • Pではなさそう(実際、Pではないことが示されています)
  • 最善手順だけ見せられても多項式時間では納得できない(他の手を指したら/打ったらどうなるの?となってしまう)
    • NPではなさそう
  • 最善手順は指数的に長い可能性があるため、交代性チューリングマシンの機能を用いても多項式時間で解くことはできなさそう(『絶対間違わないプロ棋士』同士の対戦は指数時間かかる)
    • APではなさそう
  • 最善手順は指数的に長い可能性があるため、深さ優先探索を行うと指数的な空間が必要になる
    • PSPACEではなさそう
  • 局面を覚えておくのは多項式空間でできるため、交代性チューリングマシンの機能を用いて『絶対間違わないプロ棋士』同士を戦わせることにより多項式空間で解ける
    • APSPACEに含まれる
  • 局面の情報のみから必勝手順があるかが決まり、あり得る局面は指数個しかないので、メモ化深さ優先探索すれば指数時間かけて解ける
    • EXPTIMEに含まれる

補足

日本ルールの囲碁はアゲハマの数を覚えておく必要がありますが、指数的な長さの最善手順の対局では指数的な個数のアゲハマしか発生せず、それは多項式bitの情報量しか持ちません。

中国ルールの囲碁は、指数的に長い手順が必要なパターンがあることが示されておらず、EXPTIME完全かはわかっていないようです。それでも少なくともPSPACE困難であることはわかっています。コウの取り扱いに依存しない囲碁の部分問題として、「このシチョウ(珍瓏)はどちらが勝つか?」のような決定問題が考えられます。この問題はオセロ同様に単調に終局に向かう性質があり、実際PSPACE完全です。

将棋*9をはじめとする同一局面の再現を制限するルールのある二人零和有限確定完全情報ゲーム

ボードサイズを一般化し、(任意の位置に任意の駒を配置できる)特定の局面に必勝手順があるかを判定する問題はEXPSPACE完全であることがわかっています。 同一局面を再現する手を合法手としないルールは、スーパー劫ルールとも呼ばれます。将棋の千日手ルールは、すでに三回出現した局面を再現する手を合法手としないルールですが、似たようなものです。 いずれも相手の有用な着手を妨害する効果を持つ点が話をややこしくする原因です。

  • 指数的な分岐があり、多項式時間では解けなさそう
    • Pではなさそう(実際、Pではないことが示されています)
  • 最善手順だけ見せられても多項式時間では納得できない(他の手を指したら/打ったらどうなるの?となってしまう)
    • NPではなさそう(実際、NPではないことが示されています)
  • 最善手順は指数的に長い可能性があるため、交代性チューリングマシンの機能を用いても多項式時間で解くことはできなさそう(『絶対間違わないプロ棋士』同士の対戦は指数時間かかる)
    • APではなさそう(実際、APではないことが示されています)
  • 最善手順は指数的に長い可能性があるため、深さ優先探索を行うと指数的な空間が必要になる
    • PSPACEではなさそう(実際、PSPACEではないことが示されています)
  • どの局面が既に出現したかも含めたあり得るゲームの状態の数は指数個に収まらないので、深さ優先探索は指数時間では終わらない
    • EXPTIMEではなさそう
  • どの局面が既に出現したかという指数的な情報量を覚えておく必要があるため、交代性チューリングマシンの機能を用いて『絶対間違わないプロ棋士』同士を戦わせる方法は多項式空間では足りない
    • APSPACEではなさそう
  • 最善手順は指数的な長さに収まる*10ので、深さ優先探索をすれば指数的な空間を使って解くことができる
    • EXPSPACEに含まれる
  • 最善手順は指数的な長さに収まるので、交代性チューリングマシンの機能を用いて『絶対間違わないプロ棋士』同士を戦わせる方法を使えば指数時間で解ける
    • AEXPTIMEに含まれる

補足

ところで、詰将棋にも千日手の問題がありそうに思われますが、実際はEXPTIMEに収まります。連続王手の千日手は攻方の負け=不詰であることから、攻方に不利にしか働きません。また、千日手にならない程度に同じ手順を繰り返しても手順が長くなるだけであり、やはり攻方に不利にしか働きません。ただし、双玉図式では攻方に有利に働く場合がありうる*11ため必ずしもEXPTIMEに収まるとは限らなくなります。最短路問題に置き換えてみると直感的な理解ができるかもしれません。つまり、普通の詰将棋では単に距離が正の閉路が追加されるだけで、問題が複雑になるわけではありません。一方双玉図式の場合、使える回数に制限のある負閉路が増えた状況にあたり、問題が複雑化します。

日本ルールの囲碁のコウ(劫)ルールは千日手問題と似ていますが、問題を複雑化させるに至らず、EXPTIMEに収まります。囲碁では着手によりその手番の石が盤上に必ず一つ増えることから、元の盤面に戻る着手はその増えた石をとる着手しかありえません。つまり、その増えた石はどの石なのか=直前の着手はなにか、を覚えておく必要がありますが、それ以前の情報を覚えておく必要はありません。チェスの50手ルールも、直前50手のみ覚えておけばいい点で同様です。なお、この50手というのは理論的に導かれたものではなく、実際、一般化される前の普通のチェスでは、50手ルールが悪影響を及ぼす盤面の存在が知られています。そういった問題を避けるためにこの50手という部分がボードゲームサイズに依存して決められる場合、EXPTIMEに収まるとは限らなくなります。

中国ルールの囲碁には、任意の既出局面に戻る着手は禁止される、というスーパー劫ルールが存在します。そのためEXPSPACE完全な雰囲気がありますが、そもそも指数的に長い手順が必要なパターンが知られておらず、スーパー劫ルールを持つ中国ルールの囲碁がEXPSPACE完全かはわかっていないようです。日本ルールの囲碁がEXPTIME完全であるという証明は通常のコウ(劫)ルールに依存しており、スーパー劫ルールを導入すると証明が破壊されるそうです。

余談

将棋の古い千日手ルール「仕掛けたほうが打開する」は、千日手はほとんどの場合取って取り返され打って相手が守りのために打つという一連の手順を繰り返す場合として発生することから定められたルールのようですが、序盤の駒組み(駒の取り合いが発生せず、どちらが仕掛けたともつかない状況)でも発生したことからルールが改正されました。

しかし、その改正された千日手ルール「同一手順で同一局面に戻ることが三回発生した場合」も、同一局面に戻る手順が複数存在する場合に問題となりました。

そこでさらに改正された千日手ルールでは「同一局面四回目」というルールになりました(同一局面に戻る手順が唯一の場合、等価になります)。

他のボードゲームではたいてい「同一局面三回目」というルールになっていますが、これは、同一局面に戻る手順があることは二回目には認識できるため、三回目の発生は打開しなかったことが明白だという考え方に基づくのでしょうか?

余談の余談:同一局面に戻る手順が複数存在する場合にどう問題になるのか?

そのような手順A, Bがあるとします。すると、A B BA BAAB BAABABBA BAABABBAABBABAAB ... のような手順で繰り返す場合、この手順の中にはAAAなどはもちろんのこと、AABAABAABのような一連の手順で同一局面に戻ることが三連続するパターンも存在せず、無限に続くことになります。この手順中に同一手順が三連続するパターンがない(非三連である)ことは直感的には明らかですが、具体的に示すのはちょっと難しいです。

もう少し簡単に示せる構成もあります。

  1. iが三で割って一余る数なら、i回目にはAを指す。
  2. iが三で割って二余る数なら、i回目にはBを指す。
  3. iが三で割り切れる数なら、i回目にはi/3回目と同じ手順を指す。

まず、AAAとBBBのような三手順が全て同じということは絶対にありません。ABABABのような、六手順で、二手順を三回繰り返すようなものが存在しないことも構成から明らかです。問題はABXABXABXのような九手順で、三手順を三回繰り返すようなものがないかですが、そういうのがあるとするとXが全部同じだということになり、三手順すべて同じものはありえない話と矛盾します。以下同様にして、3(3n+1)手順や3(3n+2)手順で同一手順が三連続することがないのは構成から明らかであり、9n手順で同一手順が三連続するパターンがないことは3n手順でそういったものがないことから示されます。

メモ:読んでいない(というか手に入れることができておらず読めていない)論文

  • ローグライクゲームで、「残り体力がXであるとき、脱出できるか?」はNP完全(武永、2005)
    • どういう定式化をしたのか知りたいです。敵は動いたり発生したりしない(ので指数的に長い手順は発生しない)・武器や防具が存在する(のでどの敵から倒すべきかは簡単には決まらない)、あたりでしょうか?World Wide Adventureを思い出します。
  • ローグライクゲームで、空腹度も考えるとPSPACE完全(武永ら、2007)
    • PSPACE完全ということは何かしらの可逆なギミックがあるはずで、どういう定式化をしたのかが知りたいです。無限に敵がわく・その敵の肉を食べることで空腹度を回復できる・ボス敵を倒すためには時間がかかるので途中で雑魚敵を狩って空腹度を回復する手順が挟まる(ので指数的長さの手順が発生しうる)、とかなのかな?
  • 日本ルールの囲碁はEXPTIME完全(Robson, 1983)
    • 日本ルールのコウ(劫)ルールが重要な役割を果たしているらしいのですが、あまり自明ではないため読みたいです。
  • 倉庫番問題のPSPACE完全性をGeneralized geographyの帰着により示した論文(?)
    • なんか昔見た記憶があるのですが、見つかりません(知っている方がいらっしゃいましたら教えてください)。日本語論文だったような気がしますが、うろ覚えです。NP困難性を示したものは、(略証ですが)上原氏によるもの(LA 1999)があり、これの発展形だったように思います。なお、倉庫番はPSPACE完全だと示した最初の論文(Cullberson 1998)は線形拘束オートマトンの帰着です。

*1:自分の選択のみがクリアできるかに関わる、有限確定完全情報ゲーム

*2:二人零和確定有限完全情報ゲーム

*3:余談:対数を二乗すると対数ではなくなるので、対数領域の時はこの議論は使えず、L≠NLかはわかっていません。少なくとも、L⊆NL⊆L2であることはわかっています。また、空間計算量と時間計算量の関係によりL⊆Pであることがわかります。一方、2log x多項式ですが、2(log x)2は超多項式になるため、L2はPに含まれるとは限りません。それとは無関係に、NL完全問題(例えば有向グラフの二点間到達可能性問題)を多項式時間で解く方法(例えば深さ優先探索多項式空間を使うがPなのでこれは許される)は知られているので、NL⊆Pであることがわかっています。

*4:二コリの登録商標です。

*5:メガハウス登録商標です。

*6:Demaineらにより、2016年に示されました。

*7:双玉図式では、最後の審判に代表される、ルールの不備が指摘されています。

*8:終局判定に"両者の同意"なるものが必要になるなど、定式化の障害となる点が残ります。

*9:将棋はEXPTIME完全とする論文では千日手を考慮していません。

*10:あらゆる局面は指数個しかなく、これを(三回)網羅すれば必ず終局するので、手順の長さ自体は指数的な長さに収まります。

*11:実際、最後の審判の作意手順はそれを意図したものです。

DAGのトポロジカルソートのうち最適なものを見つけたい

DAG (Directed acyclic graph) のトポロジカルソートは一般にΟ(V!)通りありますが、その中で特定の評価関数を最小にするものを見つけたい、という問題群を考えています。

最適化コンパイラを作るときに出てきた問題やコードゴルフに関連する問題も含まれているので、NP困難な問題もありそうです。

解けていない問題の概要は以下のような感じです。形式的な定義は後述しますので、用語(方言)の意味が分からなくても大丈夫です。

  1. そのトポロジカルソートの"幅"の最小化(gitのコミットロググラフであまり横幅が広がらないようにしたい。あるいは、いつかやらないといけない仕事を覚えておく脳容量を減らせるスケジューリングがしたい)
  2. そのトポロジカルソートにおける最長の辺の長さの最小化(最適化コンパイラで命令の順番をどのようにすべきかに関連)
  3. そのトポロジカルソートにおける辺の長さの総和の最小化(関数型難解言語Grassのgolfをする時に解く必要のある問題)

 

  • 1.と3.の問題はbitDPを行うことで最適解を求めることできますが、グラフのサイズはそれなりに大きい(たとえばV~1000)ことを想定しているので現実的ではありません。
  • 2.の問題は、無向グラフの場合はbandwidth問題と呼ばれる問題(指数時間アルゴリズム入門の87ページ)で、Ο(5^|V|)時間アルゴリズムを構成したという論文が2012年に出たそうです(アイデアは2008年頃が初出?同著者からもっと底が小さいアルゴリズムも提案されている)。このアイデアは有向グラフでも適用できて、有向であるという制約のおかげで指数の底が小さくなりますが、2よりは大きそうです。

いろいろ調べてみたところ、無向グラフを扱っているものは計算量理論的な議論が多く、有向グラフを扱っているものはいかに見やすいDAGの可視化をするかというヒューリスティクス的アプローチの論文が多そうに感じました。

形式的定義

グラフとトポロジカルソートについて

頂点の集合をVとし、有向辺の集合をEとします。また、有向辺から有向辺の始点と終点の順序対への写像をfとします。 この三つ組G=(f,V,E)を有向グラフと呼びます。

このグラフに閉路が存在しないとき、それは有向非巡回グラフ (DAG) と呼ばれます。

その場合、二頂点間の到達可能性を二項関係とみなすと、これは半順序となります。

トポロジカルソートとは、その半順序関係を全順序に拡張したもののことを言います。

具体的には、Vの並べ替え(全単射 s:V→\{1, 2, 3, ..., |V| \}であって、次の条件を満たすもののこととします。

 \forall v, w \ (vからwに到達可能) \rightarrow s(v) \le s(w)

一般にトポロジカルソートsは複数種類あり得ます。

参考

  • トポロジカルソートを何か一つ求める問題は、時間計算量Ο(|V|+|E|)で解くことができます。
  • トポロジカルソートが何通りあるかを調べる問題は#P完全です。いわゆるbitDPの典型問題として、Ο(|V|2^|V|)時間*1アルゴリズムが存在します。

トポロジカルソートの"幅"の最小化

【問題】DAGであるGが与えられる。 \mathop{max}\limits_{i \in \{1, 2, 3, ..., |V|-1 \}} |\{ e | e \in E, g(v) \le i \land i \lt g(w) \ where \ (v,w) = f(e)\}|を最小化するトポロジカルソートsを求めよ。

無向グラフなら、(方言ではない)グラフの言葉を使うと、Cutwidthを最小化するLayoutを求めてください、という判定版がNP完全な問題になります。

bitDPで解ける、について

トポロジカルソートを前から順番に確定させていくとします。この時、すでに順番を確定させた頂点の集合をSと呼ぶことにします。

グラフを現在の位置で切った時の"幅"とはつまりSに含まれる頂点からSの補集合に含まれる頂点への辺の数であり、それはSのみによる情報です。そこにはSの中での並べ方の情報は一切含まれていません。

そこで、Sに含まれる頂点をうまく並べた時の最大"幅"は、「現在の位置で切った時の"幅"」と「Sから一つ頂点を除いた集合に含まれる頂点ををうまい順番で並べた時の最大幅」(Sから一つ頂点を除いた集合は一般に複数あります)を全て計算して、その最大値を取ることで計算できます。

この性質を利用すると、トポロジカルソートをΟ(|V|!)通り全探索するのではなく、Vの部分集合の包含関係に関する動的計画法(いわゆるbitDP)を行うことでΟ(|V|2^|V|)時間のアルゴリズムを構成できます(Sの部分集合が2^|V|通り、各SについてSから一つ頂点を取り除いた集合はΟ(|V|)通り)。

トポロジカルソートの最長の辺の長さの最小化

トポロジカルソートにおける辺の長さ は、 E→\mathbb{N}の関数lであって、 l(e) = s(w) - s(v) \ where \ (v, w) = f(e)で与えられるものとします(方言)。

【問題】DAGであるGが与えられる。 \mathop{max}\limits_{e \in E} l(e)を最小化するトポロジカルソートsを求めよ。

無向グラフなら、(方言ではない)グラフの言葉を使うと、Bandwidthを最小化するLayoutを求めてください、という判定版がNP完全な問題になります。

無向グラフ版を指数時間で解くアルゴリズム M. Cygan, M. Pilipczuk (2012)

スカラー最適化問題を解くアルゴリズムを構成するうえで重要なテクニックに、「K以下の解があるかという決定問題を解くアルゴリズムを作る。するとKについて二分探索することで最適値を求めることができるようになる」というものがあります。この変換を行う利点は、「 K以下の解があると仮定する 。すると解をこのように構築できて…………あれ?うまくいかなかったから答えはNoです」みたいな動作をする、解の性質を用いた貪欲アルゴリズムが許されるようになる点です。普通、解の性質は解を求めるまで分かりませんから、それを使えるようになることはアルゴリズム的に有利に働くことがあります。

Bandwidth問題を解くアルゴリズムは、このテクニックを使っています。まず、最長の辺の長さをKだと仮定します。この時、トポロジカルソートの結果を格納する配列に、K+1要素ごとに仕切りを入れ、"部屋"に分けることを考えます。次に、頂点をどの"部屋"に入れるかであり得るパターンを全列挙することを考えます。何も制約がなければΟ( (|V|/K)^|V|)通りとなってしまって階乗通りと大差ありませんが、無向辺で結ばれている二つの頂点は、同じ"部屋"に入れるか二つある隣の”部屋”のどちらかに入れるかしかできません。よってパターン数はΟ(|V|/K 3^|V|)通りとなります(グラフが連結でなければ各連結成分の解を合わせれば解が求まるため、連結なグラフについての問題に帰着されます。よって連結だと仮定しても一般性を失いません。連結なら適当な全域木を考え、適当な頂点を根にとります。根は|V|/K個の部屋のどこにでも入れられます。他の頂点は全域木に含まれる無向辺の制約により3つの部屋のどれかにしか入れられません)。

部屋の割り当てが決まったからと言って、自動的に最長の辺の長さがK以下になるとは限りません。よって、頂点のその部屋内での位置を決める必要もあります。

部屋に入れるとき、インデックスの小さい方から入れていくことにします。頂点を入れていく順番はΟ(|V|!)通りありますが、実際には頂点を入れていった順番は何の情報も持っていなく、部屋に入れられた頂点の集合の区別のみが意味のある情報になっています(これは入れる順番と部屋の区切り方がうまいからです)。よって状態数はΟ(2^|V|)通りであり、特定の"部屋割り"について時間計算量Ο(|V|2^|V|)となります。

これを組み合わせれば、時間計算量Ο(|V|^2/K 6^|V|)のアルゴリズムが得られたことになります。

実際には最長の辺の長さがKを超えない状態はΟ(5^|V|)通りしかないため、bitDPではなくメモ化深さ優先探索を行うことでΟ*(5^|V|)時間で完了します。

なお、この方法は空間計算量がΟ*(2^|V|)ですが、同じ計算を何回も行うので少し無駄です。空間計算量が大きい代わりに時間計算量が小さいアルゴリズムがいくつか提案されています。

有向グラフ版を指数時間で解くアルゴリズムの案

同様の方法を使うと、部屋の割り当ては二通り(同じ部屋に入れるか次の部屋に入れるか)なのでΟ(|V|/K 2^|V|)通りになり、指数の底は4以下(3.5?)になります。

トポロジカルソートにおける辺の長さの総和の最小化

【問題】DAGであるGが与えられる。 \mathop{\sum}\limits_{e \in E} l(e)を最小化するトポロジカルソートsを求めよ。

無向グラフなら、(方言ではない)グラフの言葉を使うと、Minimum Linear Arrangement (MinLA) を最小化するLayoutを求めてください、という判定版がNP完全な問題になります。

bitDPで解ける、について

先ほどと同様にトポロジカルソートを前から順番に確定させていくとします。この時、すでに順番を確定させた頂点の集合をSと呼ぶことにします。

ここに一つ頂点を追加した時に増える辺の長さの計算には、やはりS内部の順番は無関係であることが重要な性質です(完全に自明というわけではありませんが、どのSに含まれる頂点からSの補集合に含まれる頂点への辺も区別する意味がないという性質があります)。

よってやはり先ほどと同様にbitDPを行うことができて、Ο(|V|2^|V|)時間のアルゴリズムを構成できます。

限定的な状況における解法

無向グラフ版の各種結果(参考文献1)

近似アルゴリズムについての結果

入力グラフに仮定を置いた場合の結果

  • MinLAは幾分簡単であり、木なら時間計算量Ο(|V|^(log3/log2))で、根付き木*2なら時間計算量Ο(|V|log|V|)で、それぞれ解けるアルゴリズムがあります。
  • Cutwidth問題は、最大次数3のグラフであっても判定版がNP完全となります。木なら時間計算量Ο(|V|log|V|)で解けるアルゴリズムがあります。
  • Bandwidth問題はさらに進んで二分木(最大次数3かつ木)であってすら判定版がNP完全となります。
    • 区間グラフの場合、時間計算量Ο(|V|+|E|)で解けるアルゴリズムがあります(対応する区間表現順に並べた場合が最適になります)。

固定パラメータ容易性 (Fixed Parameter Tractable, FPT)

  • Bandwidthが2のLayoutがあるか?という問題は時間計算量Ο(|V|)で解くことができます。
  • Cutwidthが2のLayoutがあるか?という問題は時間計算量Ο(|V|)で解くことができます。
  • 固定のkについて、Cutwidthがkのトポロジカルソートがあるか?という問題は時間計算量Ο(|V|^2)で解くことができます。kについてのオーダーがどのような恐ろしいものになっているかはわからないという点に注意してください。

有向グラフの場合、それとは対照的に、もしグラフの各頂点の入力辺が1本(つまり根付き木になっている)であれば、2-opt最適化を貪欲に繰り返していくことで解決できることが分かっています。

また、現実的な入力の多くでは、各頂点の入力辺が多くの場合1本、それほど少なくない割合で2本、という問題設定であり*3、入力辺が3本以上あることはないと仮定したアルゴリズムも実用上重要だと思われます。もっとも、最小有向閉路除去問題は入力辺を2本に限定しても判定版がNP完全だったはずなので、入力辺が2本という制約は本質的ではなくて、依然NP完全問題を帰着するのには十分な表現力があるのかもしれませんが……。

なお、この制約のもと、bitDPの時間計算量の多項式部分から|V|が消えますが、依然として2^|V|が立ちはだかります。

NP困難だった場合は適当な近似率保証のある近似アルゴリズムを見つけたいところです。

いずれの問題も、最悪のトポロジカルソートを持ってくると近似率がΘ(N)倍になってしまう例が簡単に作れます。

具体的には、2000頂点のグラフで、{ 1→2, 2→3, 3→4, ..., 999→1000 } ∪ { 1001→1002, 1002→1003, 1003→1004, ..., 1999→2000 } ∪ { 1→1001, 2→1002, 3→1003, ..., 1000→2000 }を辺に持つようなグラフを考えた時、1, 2, 3, 4, ..., 999, 1000, 1001, 1002, 1003, 1004, ..., 1999, 2000というトポロジカルソートは最悪であり、1, 1001, 2, 1002, 3, 1003, 4, 1004, ..., 999, 1999, 1000, 2000というトポロジカルソートは最良です。

参考文献

    1. Diaz, J.Petit, M. Serna, ACM Comput. Surv. 2002, 34, 313-356.

*1:トポロジカルソートは最大で|V|!個あるため、それを数えるには|V|log(|V|) bitの多倍長整数演算が必要になり、本当の計算量のオーダーはさらに|V|log(|V|)倍となります。ただし、適当な数で割った余りを求める問題として定式化する場合はこの問題はなくなります。あるいは、現実的な時間で解ける問題しか出題されないという制約、たとえば |V| \le 20であるとすれば64bitに収まるのでこの問題は隠されます。本文中の他の部分では、そこまで大きな整数の演算は必要ない(log(|V|) は64を超えない)という妥当な仮定のもと、log(|V|)項を無視した表記になっています。そもそも指数時間アルゴリズムなので非常に些細な違いであり、あまり気にする必要はありません。

*2:唐突に有向グラフが出てきたように思えますが、本当にあっているんでしょうか?

*3:gitのコミットグラフ、最適化コンパイラ、Grassのgolfという全く出自の違う問題群なのに、出来上がるDAGの特徴は似ている、というのはなかなか面白いことだと思います。