Rで確率分布を扱う

　勉強会で、Rで確率分布を扱う方法を教えてもらった。かなり便利だなと思った。
自分なりに問題を設定して、解いてみた。

問題１

日本では１年間（※1）に交通事故が平均80万件起こる。

1. 今日（※2）、日本で交通事故がちょうど2200件起こる確率を求めよ。
2. 今日、日本で交通事故が2200件以下起こる確率を求めよ。
3. 今日、日本で交通事故が2200件より多く起こる確率を求めよ。
4. 「１日の間で交通事故が起こる件数はx以下である。」という主張が99.9%の確率で信用できるxの値を求めよ。
5. １週間にわたる１日の交通事故件数を乱数で生成せよ。

(※1) １年＝365日とすること。

(※2)  休日か平日で状況が違うのではとかは考えなくてよいものとする。

解答１

ある期間内に事象が平均何回起こるか分かっているので、ポアソン分布を使って解きます。

> # problem 1.1
> dpois(2200, 800000/365)
[1] 0.008375249
>
> # problem 1.2
> ppois(2200, 800000/365)
[1] 0.5752186
>
> # problem 1.3
> ppois(2200, 800000/365, lower.tail=FALSE)
[1] 0.4247814
>
> # problem 1.4
> qpois(0.999, 800000/365)
[1] 2338
>
> # problem 1.5
> rpois(7, 800000/365)
[1] 2151 2202 2146 2110 2169 2197 2196

問題２

あるプロジェクトでは、バグが発見されてから次のバグが発見されるまでの期間が平均１カ月（※3）であった。不運なことに、今日バグが発見された。

1. プロジェクトマネージャはお客様に謝罪にした。「今後半年間は決してバグを出しません。もしバグが出たら罰金をお支払いします。」罰金を支払わなければならない確率を求めよ。
2. プロジェクトマネージャは無謀な宣言をしないように細心の注意を払った。「今後x日間は決してバグを出しません。もしバグが出たら罰金をお支払いします。」これは90%以上の確率で実現可能だという。xの値を求めよ。

（※3）この問題では一律で１カ月=30日とすること。

解答２

繰り返し起こる事象の到着時間が分かっているので、指数分布を使って解きます。

> # problem 2.1
> pexp(30*6, rate=1/30)
[1] 0.9975212
> # problem 2.2
> floor(qexp(0.9, rate=1/30, lower.tail=F))
[1] 3

BIT(Fenwick Tree)でRMQ

　値の更新ありのRMQは普通セグメント木で解くが、実はBITでも解ける。こんな感じ。query関数内の[l, r]がセグメントに含まれていないときの処理の仕方がかっこいい。
　計算量は、queryはO(logN)、updateはO(log²N)。

#include <iostream>
#include <cstdlib>
#include <cassert>

using namespace std;

const int N = 100000;    // # of elements in array
const int Q = 100000;    // # of queries
int a[N+1];
int bit[N+1];

int query( int l, int r ) {
    int m = 0;

    while (l <= r) {
        int p = r - (r & -r);       // p+1 is the left end of this segment
        if (p+1 >= l) {             // when [l, r] is completely included in this segment
            if(a[bit[r]] > a[m])
                m = bit[r];
            r -= r & -r;            // move to the next segment
        } else {                    // when [l, r] is not completely included in this segment
            if (a[r] > a[m])
                m = r;
            --r;                    // shrink the segment
        }
    }

    return m;
}

void update(int x) {
    for(int y = x; y <= N; y += y & -y) {
        if(bit[y] == x) {           // the maximum value was at the x-th position, but its value has changed
            int z = query(y-(y&-y)+1, y-1);
            if (a[y] > a[z])
                bit[y] = y;
            else
                bit[y] = z;
        } else {
            if(a[bit[y]] < a[x])
                bit[y] = x;
        }
    }
}

int main(int argc, char **argv) {
    srand(time(NULL));
    
    for (int i = 1; i <= N; i++) {
        a[i] = rand();
        update(i);
    }

    for (int i = 0; i < Q; i++) {
        int type = rand() % 2;

        if (type) {
            int x = rand() % N + 1;
            int v = rand();
            a[x] = v;
            update(x);
        } else {
            int l = rand() % N + 1;
            int r = rand() % N + 1;
            int x = query(l, r);
            
            for (int j = l; j <= r; j++)
                assert(a[x] >= a[j]);
        }
    }

    cout << "success" << endl;
    
    return 0;
}

Cコードのフォームの数は

数え上げの問題を考えてみました。

問題

十分に長く、十分に細い６本の指を持つ宇宙人が24fギターを弾いています。Cコードの異なるフォームの数を求めなさい。
ただし、Cコードとはドミソの３音からなる和音です。これら意外の音を鳴らしてはいけませんし、これらの３音はすべて鳴っていないといけません。
それぞれの弦は解放弦でもいいし、ミュートしてもいいこととします。

全探索解

24フレット（＝２音階分）あるのである音をならすポジションは2個（解放弦が目的の音に等しいときは3個）あります。よってそれぞれの弦について、ドミソのいずれかの音を鳴らすことのできるポジションは高々3 * 3 = 9個です。これにミュートする場合を加えると、一つの弦の有効な押え方は10個程度です。全体で10^6通りを考えればよいため、全探索でも十分計算できます。
また、音階を12を法とした正数群と捉えるとプログラミングしやすいです。

答えは110362とおりです。僕が使っているフォームは5個くらいなので、聞いたことのない響きがたくさん存在するということになります。

#include <iostream>

using namespace std;

int open[] = {4, 9, 2, 7, 11, 4};

long long ret;

void solve(int p, int mask) {
    if (p == 6) {
        if (mask+1 == (1<<3))
            ++ret;
        return;
    }

    // mute this string
    solve(p+1, mask);

    for (int i = 0; i <= 24; i++) {
        int q = (open[p] + i) % 12;
        if (q == 0)
            solve(p+1, mask|1<<0);
        else if (q == 4)
            solve(p+1, mask|1<<1);
        else if (q == 7)
            solve(p+1, mask|1<<2);
    }
}

int main() {
    solve(0, 0);
    cout << ret << endl;
    return 0;
}

包除原理を使った解

ド、ミ、ソをA, B, Cと表すことにします。また、ミュートはDで表すことにします。S(x,y, ..)を{x, y, ...}内の要素のみを用いたときに6弦を弾くパターンの数とします。

 包除原理により求めたい数は、
 S(A, B, C, D)
- S(B, C, D) - S(A, C, D) - S(A, B, D) - S(A, B, C)
+ S(A, B) + S(A, C) + S(A, D) + S(B, C) + S(B, D) + S(C, D)
-S(A) - S(B) - S(C) - S(D)　　（式１）
となります。と思っていましたが、実は違います。

上の値で求めたのは、少なくとも一つの弦がミュートされていて、かつ、コードCがなっているフォームの数です。ミュートしないフォームについては数えられていません。

ミュートしないフォームのパターン数は、
S(A, B, C) - S(A, B) - S(A, C) - S(B, C) + S(A) + S(B) + S(C)　　（式２）
です。

（式１）+ （式２）より求めたい数は、
 S(A, B, C, D)
- S(B, C, D) - S(A, C, D) - S(A, B, D)
+ S(A, D) + S(B, D) + S(C, D)
- S(D) 　　（式３）

となります。がんばれば手で計算できそうですが、プログラミングで解きました。（式３）から{A, B, C}の冪集合は打ち消しあって答えに影響しないので、Dを含まない要素はスキップしています。

#include <iostream>

using namespace std;

int cnt[][4] = {    // cnt[i] is (# of C, # of E, # of G, # of mute) in the i-th string
    {2, 3, 2, 1},
    {2, 2, 2, 1},
    {2, 2, 2, 1},
    {2, 2, 3, 1},
    {2, 2, 2, 1},
    {2, 3, 2, 1}
};

int main() {
    long long ret = 0;

    for (int k = 0; k < 1<<4; k++) {
        if (!(k & 1 << 3)) continue;     // cancelled term
        long long ptr = 1;
        for (int i = 0; i < 6; i++) {
            int sum = 0;
            for (int j = 0; j < 4; j++) {
                if (k & 1 << j)
                    sum += cnt[i][j];
            }
            ptr *= sum;
        }

        if (__builtin_popcount(k) % 2)
            ret -= ptr;
        else
            ret += ptr;
    }

    cout << ret << endl;
    return 0;
}

正多面体群の回転軸、位数、彩色パターン数について考察

　前回まではコンピュータの力を使って正多面体をt色以下で塗るパターン数を数えていました。今回は手計算で計算してみようと思います。そのためにはまず回転を表す対称群をうまく列挙するために回転軸について考える必要があります。群の元を回転軸ごとに分けて解析することで、計算がしやすくなります。

回転軸について

回転軸について考えてみたことを書きます。
回転軸を点(a, b)を結ぶ直線で表します（※a, bは正多面体の表面上の点とします）。このときa, bはどのような点か考えると、軸を簡単に探すことができます。(a, b)は

• 頂点である
• 辺上の点である
• 面上の点である

のいずれかになります。もし辺上の点だったとして、辺をx : y (x != y)に分割するような点になることがありえるでしょうか？この辺は360度回転しないと元に戻らないのでこのような点が回転軸上の点となるような群の元は存在しないと考えることができます。

辺の中点であれば、180度回転で元の形と一致するのでこれは軸上の点になりえます。

次に面上の点について考えます。辺のときと同じように考えると、点は面の重心になければいけないと分かります。この面は、360 / n度回転すると元の形に戻ります。（nは正多面体の一面にある頂点数）

最後にaまたはbが頂点のときを考えてみます。頂点は点なので軸上にあれば回転によって形がくずれません。頂点に隣接する面の数をxとすると、360 / x度回転すると、頂点のまわりの面の形が元に戻ります。

さらに、(a, b)を結ぶ直線は多面体の重心を通る必要があります。

よって回転軸は

• (頂点, 多面体の重心)を結ぶ直線
• (辺の中点, 多面体の重心)を結ぶ直線
• (面の重心, 多面体の重心)を結ぶ直線

のいずれかになります。

また、頂点, 辺の中点, 面の重心のいずれかから重心に線を引き延ばしていくと、逆側の頂点, 辺の中点, 面の重心にたどりつきます。

以上を踏まえると回転軸の数は、

(V + E + F) / 2 = (V + (V + F - 2) + F) / 2 = V + F - 1.　　（式1）

と考えられます。

正四面体の対称群について

正四面体を例に考えてみます。対称群の位数、各元のサイクル数（閉じた置換の数。軌道と呼ぶ？？）を数えて、正四面体をt色で塗るパターン数を計算してみます。

Tetrahedron, transparent, slowly turning It's a variation of Cyp's animated polyhedrons, turning by approx. 8 frames/sec instead of 20. Variated by -- Peter Steinberg, March 12 2005

1. 回転軸 = (頂点, 面の中点)の場合
この軸の選び方のパターンは頂点の数 = 4個あります。
回転角度のパターンは、120度、240度、360度の3種類です。
サイクル数は2です。

2. 回転軸 = (辺の中点, 辺の中点)の場合
この軸の選び方のパターンは辺数 / 2 = 3個あります。
回転角度のパターンは、180度、360度の2種類です。
サイクル数は2です。

回転軸の数は7個となり、(式1)と合致しています。

以上より、
Tetrahedral Coloring = (t⁴ + 8t² + 3t²) / 12 = (t⁴ + 11t²) / 12.

ちょっと詳しめに説明すると、t⁴は回転しない場合(identity element)に対する不動点の数です。8t²の部分は上記1.の場合で数えたもので、(1.の軸の選び方) * (回転角のパターン) = 4 * (3-1) = 8が係数部分です。3から1引いているのは360度回転だとidentity elementと同じだからです。t²の2はサイクル数 = 自由に選べる色の数です。
同様に、3t²は上記2.の場合で数えたものです。最後に12で割っていますが、これは群の位数です。

正六面体の対称群について

正六面体についても同様のことをやってみます。

Hexahedron, transparent and slowly turning It's a variation of Cyp's animated polyhedrons, turning by approx. 8 frames/sec instead of 20. Variated by -- Peter Steinberg, March 12 2005

1. 回転軸 = (頂点, 頂点)の場合
この軸の選び方のパターンは頂点の数/2 = 4個あります。
回転角度のパターンは、120度、240度、360度の3種類です。
サイクル数は2です。

2. 回転軸 = (辺の中点, 辺の中点)の場合

この軸の選び方のパターンは辺の数/2 = 6個あります。
回転角度のパターンは、180度、360度の2種類です。
サイクル数は3です。

3. 回転軸 = (面の重心, 面の重心)の場合

この軸の選び方のパターンは面の数/2 = 3個あります。
回転角度のパターンは、90度、180度、270度、360度の4種類です。
サイクル数は180度回転の場合は4、それ以外の場合は3です。

回転軸の数は13個となり、(式1)と合致しています。

t色以下で面を塗った場合のパターン数は、
hexahedral coloring = (t⁶ + 8t² + 6t³ + 3(t⁴ + 2t³)) / 24 = (t⁶ + 3t⁴ + 12t³ + 8t²) / 24
です。

正八面体の対称群について

次、正八面体。

Spinning octahedron, slowly turning by Cyp

1. 回転軸 = (頂点, 頂点)の場合
この軸の選び方のパターンは頂点の数/2 = 3個あります。
回転角度のパターンは、90度、180度、270度、360度の4種類です。
サイクル数は180度回転の場合は4、それ以外の場合は2です。

2. 回転軸 = (辺の中点, 辺の中点)の場合

この軸の選び方のパターンは辺の数/2 = 6個あります。
回転角度のパターンは、180度、360度の2種類です。
サイクル数は4です。

3. 回転軸 = (面の重心, 面の重心)の場合

この軸の選び方のパターンは面の数/2 = 4個あります。
回転角度のパターンは、120度、240度、360度の4種類です。
サイクル数は4です。

回転軸の数は13個となり、(式1)と合致しています。

t色以下で面を塗った場合のパターン数は、
octahedral coloring = (t⁸ + 3(t⁴ + 2t²) + 6t⁴ + 8t⁴) / 24 = (t⁸ + 17t⁴ + 6t²) / 24
です。

正十二面体の対称群について

なんか難しくなった。紙で実物をつくるのが難しいので透視図を見ながら考えてみる。

A rotating Dodecahedron. Animated GIF image. Created by en:User:Cyp and copied from the English Wikipedia.

1. 回転軸 = (頂点, 頂点)の場合
この軸の選び方のパターンは頂点の数/2 = 10個あります。
回転角度のパターンは、120度、240度、360度の3種類です。
サイクル数は12/3 = 4です。

2. 回転軸 = (辺の中点, 辺の中点)の場合

この軸の選び方のパターンは辺の数/2 = 15個あります。
回転角度のパターンは、180度、360度の2種類です。
サイクル数は12/2 = 6です。

3. 回転軸 = (面の重心, 面の重心)の場合

この軸の選び方のパターンは面の数/2 = 6個あります。
回転角度のパターンは、72度、144度、216度、288度、360度の5種類です。
サイクル数は4です。

回転軸の数は31個となり、(式1)と合致しています。

t色以下で面を塗った場合のパターン数は、
dodecahedral coloring = (t¹² + 20t⁴ + 15t⁶ + 24t⁴) / 60 = (t¹² + 15t⁶ + 44t⁴) / 60
です。

正二十面体の対称群について

最後、正二十面体。 紙で作ってクルクルしなくても解けるようになったけど、透視図か回転している映像がないとイメージ沸かないなー。平面グラフに直して二次元で解く方法とかありそうだけどなー。無いのかなー。

A rotating Icosahedron. Animated GIF image. Created by en:User:Cyp and copied from the English Wikipedia. Text from en: Spinning icosahedron, made by me using POV-Ray, see en:image:poly.pov for source.

1. 回転軸 = (頂点, 頂点)の場合
この軸の選び方のパターンは頂点の数/2 = 6個あります。
回転角度のパターンは、72度、144度、216度、288度、360度の5種類です。
サイクル数は4です。

2. 回転軸 = (辺の中点, 辺の中点)の場合

この軸の選び方のパターンは辺の数/2 = 15個あります。
回転角度のパターンは、180度、360度の2種類です。
サイクル数は20/2 = 10です。

3. 回転軸 = (面の重心, 面の重心)の場合

この軸の選び方のパターンは面の数/2 = 10個あります。
回転角度のパターンは、120度、240度、360度の3種類です。
サイクル数は2 + 18/3 = 8です。

回転軸の数は31個となり、(式1)と合致しています。

t色以下で面を塗った場合のパターン数は、
icosahedron coloring = (t²⁰ + 24t⁴ + 15t¹⁰ + 20t⁸) / 60 = (t²⁰ + 15t¹⁰ + 20t⁸ + 24t⁴) / 60
です。

よーし、全部解けた。答え[1]もあってた。

あと(式1)で予想した回転軸の数もすべての正多面体でマッチした。

参考URL

[1] Polyhedron Coloring -- from Wolfram MathWorld

正八面体をn色で塗るパターン数

　正八面体（regular octahedron）をn色で塗るパターン数をポリアの数え上げで数えてみました。もはや正多面体に愛情すら感じるようになってきました。まさにプラトニックラブです。

前回の隣接する面が写像前後で変わらないという性質を使って、対称群（symmetric group）から正多面体の回転に対応する元のみ抽出するという技を応用して正八面体に適用してみました。

対称群すべて列挙して・・というやり方ではなくて、DFSで枝刈りしながら探索としてみました。かなり大きな枝刈りが出来ていると思うので、この方法を使えば正二十面体も攻略できそうです。

検算の結果[1, 2]、正しく計算できていました。

#include <cassert>
#include <iostream>
#include <cstring>
#include <vector>

using namespace std;

// faces that shares an edges (clockwise)
int adjFaces[][3] = {
    {3, 1, 7},
    {0, 2, 4},
    {1, 3, 5},
    {0, 6, 2},
    {1, 5, 7},
    {2, 6, 4},
    {3, 7, 5},
    {0, 4, 6}
};

vector<vector<int> > group;

bool isRotate(int a[]) {
    for (int i = 0; i < 8; i++) {
        if (a[i] == -1)          // this face is not aligned yet
            continue;

        int p = -1;
        for (int j = 0; j < 3; j++) {
            if (a[adjFaces[i][j]] != -1) {
                p = j;
                break;
            }
        }

        if (p == -1)              // no adjacent faces are not aligned yet
            continue;

        int q = -1;
        for (int j = 0; j < 3; j++) {
            if (a[adjFaces[i][p]] == adjFaces[a[i]][j]) {
                q = j;
                break;
            }
        }
        
        if (q == -1)              // wrong alignment
            return false;

        for (int j = 0; j < 3; j++) {
            if (a[adjFaces[i][(p+j)%3]] == -1)   // not aligned yet
                continue;
            if (a[adjFaces[i][(p+j)%3]] != adjFaces[a[i]][(q+j)%3])
                return false;
        }
    }

    return true;
}

void dfs(int p[], int k, int mask) {
    if (!isRotate(p))
        return;

    if (__builtin_popcountll(mask) == 8)
        group.push_back(vector<int>(p, p+8));

    for (int i = 0; i < 8; i++) {
        if (mask >> i & 1)
            continue;
        p[k] = i;
        dfs(p, k+1, mask|1<<i);
        p[k] = -1;
    }
}

long long pow(long long x, long long p) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x;
        x = x * x;
        p >>= 1;
    }
    return ret;
}

int cycleNum(const vector<int> &p) {
    bool vis[8] = {};
    
    int ret = 0;
    for (int i = 0; i < 8; i++) {
        if (vis[i])
            continue;
        ++ret;
        for (int j = i; !vis[j]; j = p[j])
            vis[j] = true;
    }

    return ret;
}

long long solve(int color) {
    int p[8];
    memset(p, -1, sizeof(p));

    dfs(p, 0, 0);
    assert(group.size() == 24);
    
    long long sum = 0;
    for (size_t i = 0; i < group.size(); i++) {
        sum += pow(color, cycleNum(group[i]));
    }

    return sum / group.size();
}

long long getAnswer(int c) {
    return (6*pow(c, 2) + 17*pow(c, 4) + pow(c, 8)) / 24;
} 

int main(int argc, char **argv) {
    int c;
    cin >> c;

    long long ret = solve(c);
    assert(ret == getAnswer(c));
    cout << ret << endl;

    return 0;
}

参考URL

[1] Polyhedron Coloring -- from Wolfram MathWorld
[2] Octahedral symmetry - Wikipedia, the free encyclopedia

正多面体の頂点・辺・面の数を数える

　前回正六面体をn色で塗るパターンを計算機を使って数えましたが、どうやら手計算で出来るっぽいです。
すぐには分からないので、とりあえず正多面体（Platonic Solid）について考察することから始めてみました。

• 正多面体は何種類あるのか？
• 各正多面体の頂点・辺・面の数はいくつなのか？

を考えてみました。

正多面体はいくつあるか？

一つの頂点に接している面がa個あるとします。このとき
a >= 3　　（式1）
となります。a=1, a=2だと立体にならないからです。
正多面体の面が正n角形であるとします。このとき正n角形の内角は180 - 360/n 度になることから、
a (180 - 360/n) < 360　　（式2）
という式が満たされなければならないことが分かります。当然nは、
n >= 3　　（式3）
です。
式1-3を満たす整数(n, a)のペアは、(3, 3)、(3, 4)、(3, 5)、(4, 3)、(5, 3)となります。

正多面体の頂点・辺・面の数は？

これは、オイラーの多面体定理

V - E + F = 2

を使えば解けます。ここで、V, E, Fはそれぞれ、頂点、辺、面の数です。

面が正n角形で、頂点に接する面の数がaであるような正多面体について考えてみます。
面の数をxとおくと、頂点の数はn*x/a（※1）、辺の数はn*x/2（※2）となります。

（※1）1つの面に頂点はn個。面はx個あるのでn*xとなるが、これだと1つの頂点がa回重複して数えられるのでaで割る。
（※2）1つの面に辺はn個。面はx個あるのでn*xとなるが、これだと1つの辺が2回重複して数えられるので2で割る。

これをオイラーの多面体定理に代入すると、
n*x/a - n*x/2 + x = 2
となります。

これを解くと、
V = 4n / (2n - an + 2a)
E = 2an / (2n - an + 2a)
F = 4a / (2n - an + 2a)

となります。

正六面体をn色で塗るパターン数

ポリアの数え上げを使って、正六面体をn色で塗るパターン（回転して同じになるものは同じとみなす）を数えてみました。

回転群の要素を列挙するクールな方法が思いつかなかったのですが、面の数が６個なので以下のようにしました。

(0,1,2,3,4,5)の順列をすべて挙げ、それが正六面体の回転群の元になっているかどうか判定する。
判定条件は以下のとおり。
『すべての面fについて、fを手前に向けてfに接する面を時計回りに挙げていったときに得られる列が写像前と写像後で（シフトすれば）一致していれば、回転群の要素である。』

上の判定で正しく回転群を抽出できていることは以下のように確認しました。

回転群の集合をA, 上記の判定条件を満たす集合をBとする。
要素xがAに含まれるなら、xはBに含まれる。が言えるため、AはBの部分集合になる。
Aの要素数は、上面にどれを持ってくるか * 側面の回転 = 6 * 4 = 24となる。
つまり、Bの要素数が24となれば、BとAは同じ集合であると言える。実際にBの要素数は24となったためA = B。

以上をふまえて書いたプログラムが以下です。
実は今回計算した数は多項式で計算できるようなので、それを使って検算しました。
 正しく計算できてました。

#include <algorithm>
#include <cassert>
#include <iostream>

using namespace std;

int adjacentFaces[][4] = {
    {1, 4, 3, 5}, 
    {2, 4, 0, 5}, 
    {3, 4, 1, 5}, 
    {0, 4, 2, 5}, 
    {1, 2, 3, 0}, 
    {3, 2, 1, 0}
};

bool isRotate(int a[]) {
    for (int i = 0; i < 6; i++) {
        int p = -1;
        for (int j = 0; j < 4; j++) {
            if (a[adjacentFaces[i][j]] == adjacentFaces[a[i]][0]) {
                p = j;
                break;
            }
        }
        if (p == -1)
            return false;

        for (int j = 0; j < 4; j++)
            if (adjacentFaces[a[i]][j] != a[adjacentFaces[i][(p+j)%4]])
                return false;
    }
    return true;
}

long long pow(long long x, long long p) {
    long long ret = 1;
    while (p > 0) {
        if (p & 1)
            ret *= x;
        x = x * x;
        p >>= 1;
    }
    return ret;
}

int cycleNum(int p[]) {
    bool vis[6] = {};
    
    int ret = 0;
    for (int i = 0; i < 6; i++) {
        if (vis[i])
            continue;
        ++ret;
        for (int j = i; !vis[j]; j = p[j])
            vis[j] = true;
    }

    return ret;
}

long long solve(int color) {
    int p[] = {0, 1, 2, 3, 4, 5};

    int cnt = 0;
    long long sum = 0;
    do {
        if (!isRotate(p))
            continue;
        ++cnt;
        sum += pow(color, cycleNum(p));

    } while (next_permutation(p, p+6));

    return sum / cnt;
}

int main(int argc, char **argv) {
    int c;
    cin >> c;
    
    long long ret = solve(c);
    assert(ret == (pow(c, 6) + 3*pow(c, 4) + 12*pow(c, 3) + 8*pow(c, 2)) / 24);
    cout << ret << endl;
    
    return 0;
}

GDBでC++デバッグするときはpretty-printersが便利

pretty-printersというツールを使うと、GDBデバッグ時にSTLコンテナの中身を”人間にも読める形”で表示させることが出来ます。

Before --pretty-printersを使わない場合

以下のソースをサンプルとして使います。

int main(int argc, char **argv) {
    int buf[] = {1,2,3,4,5};
    vector<int> vec(buf, buf+5);

    cout << accumulate(buf, buf+5, 0) << endl;
    cout << accumulate(vec.begin(), vec.end(), 0) << endl;

    return 0;
}

デバッグ中に、bufおよびvecの中身を表示すると以下のようになります。

(gdb) p buf
$1 = {1, 2, 3, 4, 5}
(gdb) p vec
$2 = {<std::_Vector_base<int, std::allocator<int> >> = {_M_impl = {<std::allocator<int>> = {<__gnu_cxx::new_allocator<int>> = {<No data fields>}, <No data fields>}, _M_start = 0x804c008, _M_finish = 0x804c01c, _M_end_of_storage = 0x804c01c}}, <No data fields>}

配列は中身が表示されますが、vectorはコンテナに格納した値が表示されません。
p *(vec._M_impl._M_start)@5とすれば中身が見れますが、毎回タイプするのは面倒です。

After --pretty-printersを使った場合

pretty-printersを使うと、先ほどの例がこうなります。

(gdb) p buf
$1 = {1, 2, 3, 4, 5}
(gdb) p vec
$2 = std::vector of length 5, capacity 5 = {1, 2, 3, 4, 5}

きちんと欲しい情報が表示されました。
他にもいろいろ試してみましたが、vector<vector<int> >やset<pair<int, int> >などの複雑な型にも対応しているようです。

インストール方法

1. svnからチェックアウト

svn co svn://gcc.gnu.org/svn/gcc/trunk/libstdc++-v3/python

2. ~/.gdbinitに以下を追記

python
import sys
sys.path.insert(0, '/home/maude/gdb_printers/python')
from libstdcxx.v6.printers import register_libstdcxx_printers
register_libstdcxx_printers (None)
end

※赤字のところは、チェックアウトしたディレクトリに置き換えてください。

参考URL

[1] debugging - How do I print the elements of a C++ vector in GDB? - Stack Overflow
[2] STLSupport - GDB Wiki