Codeforces Round #174 Cows and Primitive Roots

Div2 Aだけど面白い問題。

問題概要

素数pが与えられる。法pにおける原始根の数を数えよ。(p < 2000)

解法

各整数1 <= x < pについて、原始根になるかどうか愚直に調べればいい。計算量はO(p²)。
これで十分間に合うがもっと効率的な解法がある。

別解

ここにあるとおり。

まず、pは素数なので必ず少なくとも１つの原始根が存在する。それをgとおく。
gは原始根なので、{g, g², g³,...,g^p-1}と{1,2,...p-1}はmod pにおいて同じ集合を表す。
よって、{1,2,...,p-1}の中で原始根になる数の個数と{g, g², g³,...,g^p-1}の中で原始根になる数の個数は同じ。

{g, g², g³,...,g^p-1}の中で原始根になるのはどれか考える。
gcd(i, p-1) > 1であるとき、g^{i * (p-1)/gcd(i, p-1)} = g^{(p-1) * i/gcd(i, p-1)} = 1^{i/gcd(i, p-1)} = 1
なので、gⁱは原始根ではない。

gcd(i, p-1) = 1のとき、(p-1) | (i * x) となるような最小のxは(p-1)なので、gⁱは(p-1)乗したときにはじめて1になる。よってこれは原始根。

以上のことから、φ(p-1)が答えとなる。ここでφはオイラーのtotient function。
計算量はsqrt(p)。

Codeforces Round #173 Yet Another Number Game

問題概要

n個の非負整数が与えられる（1 <= n <= 3）。この非負整数を使って２人がゲームをする。それぞれのプレイヤーは自分の手番にて以下のいずれかの処理を行う。

1. n個の非負整数から任意の正数を選ぶ（これをxとおく）。その正数から1以上x以下の正数を引く。
2. n個の非負正数から最小のものを選ぶ（これをxとおく）。すべての非負正数から1以上x以下の正数を引く。

これを交互に繰り返し、いずれの操作も行えなくなった方が負け。お互い最善手を指したときどちらが勝つか答えよ。

解法

最初、動的計画法でといたけど間に合わなかった。以下のように場合分けするとうまくいく。

■n=1のとき
これは自明。非負整数が0であれば先手の負け。そうでなければ先手の勝ち。

■n = 2のとき
これは動的計画法で解ける。

■n = 3のとき
Nimと同様になる。

もうちょっと掘り下げて考えてみる。n=2のときなぜNimと同じにならないか？これはNimにおける負けポジションから負けポジションへ遷移することができるから。例えば、(2, 2)という負けポジションで回ってきても処理2を使って(1, 1)や(0, 0)という負けポジションで相手に返すことができる。

じゃあn=3のときはなぜNimでうまくいくか？
まずNimにおける勝ちポジションからは、処理1を使って負けポジションへ遷移することができる。
よって、Nimにおける負けポジションからは、処理2を使っても負けポジションに遷移させることはできないことを示せばいい。

処理2で引く数をxとおく。このときxを二進数で表したとき1番左にある1のビットを考える。このビットをbビットと呼ぶことにする。
今ゲームで使用している3つの非負整数x, y, zのbビット目のパターンは、負けポジションなので、(0, 0, 0)か(0, 1, 1)である。

■(0, 0, 0)のとき
(0, 0, 0) → (1, 1, 1)となるので、x ^ y ^ z = 0とはならない。

■(0, 1, 1)のとき
(0, 1, 1) → (1, 0, 0)となるので、x ^ y ^ z = 0とはならない。

よって、n = 3のときは処理2を使ってもNimの負けポジションから負けポジションに遷移することはない。したがってn = 3のときはNimと同様に考えることができる。

ソースコード

#define REP(i,n) for(int i=0; i<(int)(n); i++)
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)

bool dp[300][300];

void init() {
    REP (i, 300) REP (j, 300) {
        dp[i][j] = false;
        REP (k, i) if (!dp[k][j]) dp[i][j] = true;
        REP (k, j) if (!dp[i][k]) dp[i][j] = true;
        int x = min(i, j);
        FOR (k, 1, x+1) if (!dp[i-k][j-k]) dp[i][j] = true;
    }
}

int main() {
    int n;
    cin >> n;
    vector<int> a(n);
    REP (i, n)
        cin >> a[i];
    init();
    if (n == 1) {
        if (a[0])
            puts("BitLGM");
        else
            puts("BitAryo");
    }
    else if (n == 2) {
        if (dp[a[0]][a[1]])
            puts("BitLGM");
        else
            puts("BitAryo");
    }
    else {
        if (a[0] ^ a[1] ^ a[2])
            puts("BitLGM");
        else
            puts("BitAryo");
    }
    return 0;
}

2013年の目標

もう3月ですが、2013年の目標。

英語

• 週に4話英語の漫画を読む。
• 「めぞん一刻」英語版アニメを全話見る。
• ESL PODのビハインド分（1年半分くらい）をきちんとやる。

競技プログラミング

• 土日でその週に行われたCodeforcesの問題をすべて(div2のみ)解く。

英語はいい感じに出来ているが、競技プログラミングはダメダメだ。今週からは、土日は競技プログラミングに専念しようと思う。

Kaggle Digit Recognizer

Kaggleというパターン認識のcompetitionに参加してみました。

入門用の「Digit Recognizer」という問題にチャレンジしました。名前のとおり手書きの数字を判定するという問題です。

https://www.kaggle.com/c/digit-recognizer

有名なデータセットを使っているらしく、どのアルゴリズムを使うとどれくらいの識別率がえられるかという情報が下のページに載っています。

http://yann.lecun.com/exdb/mnist/index.html

お試しのsubmitということで、K-Nearest Neighborsを使って解いてみました。

結果は識別率96.8%。
ピクセルの明暗の値を特徴ベクトルに使っているだけなのに、こんなに高い識別率が出せたのには驚きました。

しかし、Kaggleなかなかいいですねー。賞金300万ドルの問題とかあるし・・。
次の機械学習勉強会のネタにしようかな。

Octaveで非線形SVMを解く

前回の線形SVMに引き続き、今回はカーネルトリックを使って非線形SVMを解いてみます。
カーネルを使うために主問題を双対問題に変換します。下のプログラムを見ると分かるように、主問題よりも双対問題の方がQPソルバーに乗せるのは簡単です。

#--------------------------------------------------------------------
# procedure:
#   solve a dual svm
# input: 
#   x feature vectors of training data
#   y labels of training data
#   C upper bound for alpha
#   kernel kernel function to be used
# output:
#   ret [w; b] where w x + b = 0 is the obtained hyperplane
# -------------------------------------------------------------------
function ret = dsvm(x, y, C, kernel)
    l = rows(x);           # number of training data

    x0 = zeros(l, 1);
    H = zeros(l);
    for i = 1:l
        for j = 1:l
            H(i, j) = y(i) * y(j) * kernel(x(i,:), x(j,:));
        endfor
    endfor
    q = -ones(l, 1);
    A = transpose(y);
    b = [0];
    lb = zeros(l, 1);
    ub = C * ones(l, 1);
    Al = [-inf];
    Ai = ones(1, l);
    Au = [inf];
    
    [sol, obj, info, lambda] = qp (x0, H, q, A, b, lb, ub, Al, Ai, Au);
    ret = sol;
    
endfunction;

線形分離不可能な学習データを学習させてみました。
まずは、円を使うと分けられるパターン。カーネル関数にはK(u, v) = (uv + 1)^2を使いました。

次に円でも分離できないパターン。手裏剣のような形を使うとうまく分離できそうです。
ガウシアンカーネルK(u, v) = exp(- |u-v| / 20)を使って学習してみました。

うお、そう分けるか。。すごい。

Octaveで線形SVMを解く

OctaveでハードマージンSVMを解くプログラムを書いてみました。主問題を解きます。

#--------------------------------------------------------------------
# procedure:
#   solve a hard-margin prime svm
# input: 
#   x feature vectors of training data
#   y labels of training data
# output:
#   ret [w; b] where w x + b = 0 is the obtained hyperplane
# -------------------------------------------------------------------
function ret = svm(x, y)
    l = rows(x);           # number of training data
    n = columns(x);        # demension of feature vectors

    x0 = zeros(n+1, 1);
    H = eye(n+1);
    H(n+1, n+1) = 0;
    q = zeros(n+1, 1);
    A = [];
    b = [];
    lb = [];
    ub = [];
    Ai = [x, ones(l, 1)];
    Al = -inf * ones(l, 1);
    Au = inf * ones(l, 1);    
    for i = 1:l
        if y(i) == 1
            Al(i) = 1;
        else
            Au(i) = -1; 
        endif  
    endfor

    [sol, obj, info, lambda] = qp (x0, H, q, A, b, lb, ub, Al, Ai, Au);
    ret = sol;
    
endfunction;

動作確認用のユーティリティ関数を書いてみました。
まず学習データをプロットする関数。

#--------------------------------------------------------------------
# procedure:
#   plot training data
# input: 
#   x feature vectors of training data
#   y label of training data
# output:
#   n/a
# -------------------------------------------------------------------
function plotData(x, y)
    l = rows(x);  # number of training data
    hold on;
    for i = 1:l
        if y(i) == 1
            plot (x(i, 1), x(i, 2), 'ro', 'markersize', 10);
        else
            plot (x(i, 1), x(i, 2), 'go', 'markersize', 10);
        endif
    endfor
endfunction;

そして、直線wx + b = 0をプロットする関数。

#--------------------------------------------------------------------
# procedure:
#   plot a line w x + b = 0
# input: 
#   w weight vector
#   b bias
# output:
#   n/a
# -------------------------------------------------------------------
function plotLine(w, b)
    ezplot (@(x, y) w(1) * x + w(2) * y + b);
    title ("");
endfunction;

使ってみます。

octave> x = [
> 0 0 
> 1 0
> 0 1
> ];
octave> y = [
> 1
> -1
> -1
> ];
octave> plotData(x, y);
octave> ret = svm(x, y);
octave> plotLine(ret(1:2), ret(3));

とやると、以下のようなグラフが表示されます。

ふむ。ちゃんと出来てるっぽいですね。
他の学習データパターンでもやってみました。

Octaveメモ（4）

f(x, y) = 0をプロットする方法。

ezplotを使う。例えば2*x + y - 10 = 0を範囲[-10, 10]でプロットしたい場合は、

octave> ezplot (@(x, y) 2*x + y - 10, [-10, 10]);

とする。

 関数fは線形じゃなくてもいい。

octave> ezplot (@(x, y) x.^2 + y.^2 - 1);
octave> axis square