Page List

Search on the blog

2010年11月27日土曜日

空間認知能力

今日練習した問題で面白いのがあったので紹介。

問題はこちら。

簡単に訳するとこんな感じ。
立方体をn*n*n個の小さな立方体に分ける。そのあと、適当に小さな立方体を抜き取る。
そして、この立方体の写真を三方向(x軸、y軸、z軸方向)から取る。
3枚の写真を見て、抜き取られた可能性のある小さな立方体の最小個数を求めよ。

最近買ったルーブックキューブと睨めっこ・・・。
結局解けなかった。。この問題は、テクニックというより、生まれ持った才能による部分が大きい気がする・・・・。そういえば、高校のとき「東大に行きたければ、数学は空間系の問題に強くならないとダメだ。」って先生が言ってたのを思い出した。

まー、話を戻して、友達の回答を見て自分なりに書いたコードはこんな感じ。

int removeCubes(vector<string> A, vector<string> B, vector<string> C) {
// A x-y
// B x-z
// C y-z

int sz = (int)A.size();

memset(cube, 1, sizeof(cube));
REP(i, sz) {
REP(j, sz) {
if (A[i][j] == 'N') REP(k, sz) cube[i][j][k] = 0;
if (B[i][j] == 'N') REP(k, sz) cube[i][k][j] = 0;
if (C[i][j] == 'N') REP(k, sz) cube[k][i][j] = 0;
}
}

// valid-check
REP(i, sz) {
REP(j, sz) {
if (A[i][j] == 'Y') {
int t = 0;
REP(k, sz)
t += cube[i][j][k];
if (!t) return -1;
}
if (B[i][j] == 'Y') {
int t = 0;
REP(k, sz)
t += cube[i][k][j];
if (!t) return -1;
}
if (C[i][j] == 'Y') {
int t = 0;
REP(k, sz)
t += cube[k][i][j];
if (!t) return -1;
}
}
}

int ret = 0;
REP(i, sz)REP(j,sz)REP(k,sz)
if (!cube[i][j][k])
++ret;

return ret;
}

なるほど!そういうことか!
ちょっと気になったのが、valid-checkのところ。なんか短くできそうな・・。
Editorialに赤い人のコードが載っていたので、それを参考に書いてみた。
こんな感じ。


int removeCubes(vector<string> A, vector<string> B, vector<string> C) {
int sz = (int)A.size();
bool a[sz][sz], b[sz][sz], c[sz][sz];

memset(a, 0, sizeof(a));
memset(b, 0, sizeof(b));
memset(c, 0, sizeof(c));

int ret = sz*sz*sz;
REP(i,sz)REP(j,sz)REP(k,sz)
if (A[i][j] == 'Y' && B[i][k] == 'Y' && C[j][k] == 'Y') {
--ret;
a[i][j] = 1;
b[i][k] = 1;
c[j][k] = 1;
}


REP(i,sz)REP(j,sz) {
if (A[i][j] == 'Y' && !a[i][j]) return -1;
if (B[i][j] == 'Y' && !b[i][j]) return -1;
if (C[i][j] == 'Y' && !c[i][j]) return -1;
}

return ret;

}
ほう、逆の発想をすることで、判定のところがシンプルに書けるわけか。。結構このコード理解するのに時間かかりました(笑)

ちょっと、この才能の差を埋めるにはかなり時間がかかりそうだな~~~。。。


2010年11月24日水曜日

誤差にまつわるエトセトラ

今日は、浮動小数点小数の数値計算誤差について書こうと思う。


まず、以下のコードを実行してみましょう。



#define EPS 1e-7

int main() {
// case 1
double x = 1e16;
printf("%lf\n", x + 1);

// case 2
double y1 = 123456123456.1234588623046875;
double y2= 123456123456.1234741210937500;

printf("%d\n", ABS(y1 - y2) < EPS);
printf("%d\n", ABS(y1 - y2) < EPS * y1);

// case 3
double z1 = 1e-1072;
double z2 = -1e-1072;
printf("%d\n", ABS(z1 - z2) < EPS * z1);

return 0;
}



結果は以下のようになります。


>10000000000000000.000000
0
>1
>0
>1.000000


まず、case1から見ていきましょう。
10000000000000000.000000に1を足しても
10000000000000001.000000にはなりません。
これは、doubleの有効桁数が2^52 ~10^15だからです。最初の15桁以下は丸められてしまいます。逆に言うと、10^15以下の整数ならdoubleで正確に表すことができます。
この辺の話は昔投稿しているので、こちらをご覧ください。


次に、case2です。
これは、おもしろいです。y1とy2をビット列で表すとそれぞれ以下のようになります。


01000010 00111100 10111110 10001110 11110010 01000000 00011111 10011011
01000010 00111100 10111110 10001110 11110010 01000000 00011111 10011100


実は、この2つの数は、doubleの世界では連続する数値です。それにも関らず、(y1-y2)は1e-5程度です。
これは、大変です。連続する値なので、ほんの少しのエラーで、ある値がもう片方の値になりえます。
数値計算を実施する際は、この2つは、同じ値とみなしてよいでしょう。
しかし、ABS(y1-y2) < eps (eps = 1e-7)
なんてやっちゃうと、この2つの値は、同値とは判定できません。

何がまずかったのでしょう・・。
doubleは浮動小数点なので、小数点の位置は固定ではありません。つまり、絶対的な値を誤差の許容範囲とするには無理があります。
ここでは、
ABS(y1-y2) < eps * y1
とするとよいでしょう。

最後にcase3。これは、やっかいです。case2で
ABS(y1-y2) < eps * y1
を判定条件に用いればよいことを示しましたが、case3では、この判定条件は適切ではありません。
同値とみなすべきz1とz2が同値とみなされません。
符号が異なる2つの微量な数値を比較する場合は、case2の場合は適さないようです。
では、どうすればいいか・・・。
case2とcase3を2つ使って判定しましょう。結構面倒ですが・・。。

競技系プログラミングでは、小数が現れないように分母を他辺に掛けるテクニックが有効のようです。また、比較する値の絶対値が比較的小さいものであれば、絶対評価(case2)だけ十分でしょう。


出典)
http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=integersReals
読んで理解したことを書こうと思ったら、ほぼ和訳みたいになりました(笑)
英語が出来る人は、上の原文を読んだ方が分かりやすいはずです。

2010年11月20日土曜日

C++で関数型プログラミング:accumulate編

マルチパラダイム言語として知られるC++。

競技系プログラミングばかりやっていると、手続き型でばかり書いてしまう。
職場では、ぱっと見オブジェクト指向のようなスパゲッティソースを保守している。。

最近、C++で関数型プログラミングっぽいこともできることを発見。
今日は「畳み込み関数」をC++でやってみます。
畳み込み関数は、Pythonでいうreduce()です。。

では、早速コードと実行結果を。

#include<numeric>

#define INF 999999999

int x[] = {1,2,3,4,5,6,7,8,9,10};
int y[] = {15, 12, 99, 27};
int z[] = {10, 20, 150, 100};
int w[] = {1,3,10,100, -12, 2, 4};

int multiply(int x, int y) {
return x*y;
}

int gcd(int a, int b) {
if (!b)
return a;
return gcd(b, a%b);
}

int lcd(int a, int b) {
return a*b/gcd(a,b);
}

int myMin(int a, int b) {
return (a < b) ? a : b;
}

int myMax(int a, int b) {
return (a > b) ? a : b;
}

int main() {
cout << accumulate(x, x+SIZE(x), 0) << endl;
cout << accumulate(x, x+SIZE(x), 1, multiply) << endl;
cout << accumulate(y, y+SIZE(y), *y, gcd) << endl;
cout << accumulate(z, z+SIZE(z), *z, lcd) << endl;
cout << accumulate(w, w+SIZE(w), INF, myMin) << endl;
cout << accumulate(w, w+SIZE(w), -INF, myMax) << endl;

return 0;
}

[実行結果]
55
3628800
3
300
-12
100

accumulate()のシグネチャーはこんな感じ。
T accumulate(InputIterator first, InputIterator last, T init, BinaryOperation binary_op );
  1. input は畳み込み対象の開始位置
  2. lastは畳み込み対象の終了位置
  3. initは初期値
  4. binary_opは引数を2つもつ関数のポインタ
です。第4引数を省略すると、operator+()が畳み込み関数として使用されます。
注意ポイントは、accumulate()を使用する場合は、numericをincludeしないといけないという点です。algorithmではないので間違えないように!!

やばいぞ、C++。C++好き度が上がりました。。

2010年11月15日月曜日

トポロジカルソート

今日は、SRMの過去問をトポロジカルソートで解いてみた。
有向グラフのトポロジーに注目したソートですが、トポロジーって一体・・・、
なんか他にも似たような言葉がたくさんあるような・・・。

ちょっとまとめてみます。
  • トポロジー   幾何学的な相対位置、ネットワークの接続状況など
  • エントロピー 物事の煩雑さ
  • オントロジー 物事の概念を意味・関係などで定義したもの
で、今回はトポロジカルソートです。
簡単に言うと、どの接点も自分の子孫より先に来るようにソートすることです。
Makeでコンパイルする順序を決めたり、PERTのスケジュール管理で使われるみたいです。

ではさっそく解いてみます。
問題は、こちら。

ソースは、これ。
定型的なアルゴリズムがあるようですが、とりあえず、自己流で出力辺の無い接点を取り出して、その接点はグラフ集合からのぞいて・・・を繰り返しています。



bool dn[64];
class CorporationSalary {
public:
long long totalSalary(vector<string> relations) {
memset(dn, 0, sizeof(dn));

VI idx;
int sz = relations.size();
while (accumulate(dn, dn+sz, 0) < sz) {
REP(i, relations.size()) {
if (dn[i]) continue;
bool lst = true;
REP(j, relations[i].size()) {
if (dn[j]) continue;
if (relations[i][j] == 'Y') {
lst = false;
break;
}
}
if (lst) {
idx.PB(i);
dn[i] = 1;
}
}
}

long long ret = 0LL;
long long sl[sz];
fill(sl, sl+sz, 0LL);

REP(i, sz) {
long long m = 0LL;
REP(j, sz) {
if (relations[idx[i]][j] == 'Y')
m += sl[j];
}
sl[idx[i]] = m ? m : 1;
ret += sl[idx[i]];
}
return ret;
}
};

2010年11月14日日曜日

知ってると便利なSTL(5) fill, fill_n

今日は、便利なSTLについて。
  • fill
  • fill_n
を紹介します。その名の通り、配列やコンテナをある値で”埋める”ことができる関数です。初期化のときに重宝します。

同じような標準関数にmemset()がありますが、こちらは1byte単位で値を初期化するため、charの初期化には有効ですが、intの初期化にはあまり都合がよくありません。(0で初期化するのであれば問題ありませんが、1や-1で初期化するのは面倒。詳しくは、memset()とmemcpy()を参照。)

これに対して、fill()、fill_n()は、その型に対応した単位で初期化を実行できるためとても便利です。

それでは、実際に動かしてみます。
まずは、fill()から。
void fill ( ForwardIterator first, ForwardIterator last, const T& value );
という形で使用し、[first, last)区間の要素をvalueにセットします。



int main() {
int x[10];
vector<string> names(7);

fill(x, x+10, 100);
fill(names.begin(), names.end(), "ken-ken");

REP(i, 10)
printf("%d\n", x[i]);

REP(i, 7)
printf("%s\n", names[i].c_str());

return 0;
}



次に、fill_n。
void fill_n ( OutputIterator first, Size n, const T& value );
という形で使用し、[first, first+n)までの要素をvalueにセットします。


int main() {
int x[10];
vector<int> y(10);
char z[10];

fill_n(x, SIZE(x), 10);
fill_n(y.begin(), 10, -1);
fill_n(z, 10, 'o');

REP(i, 10)
printf("%d\n", x[i]);

REP(i, 10)
printf("%d\n", y[i]);

REP(i, 10)
printf("%c\n", z[i]);

return 0;
}


要素数が可変なコンテナ(vectorなど)に対してはfill()を、要素数が固定長の配列の場合は、fill_n()を使うといいかと思います。

2010年11月10日水曜日

初めての・・・・Hard AC!

人生初の、TopCoder Hard(Division 2)自力AC。
直観でこの問題は簡単と思って、1時間もあれば・・って思ったけど、落とし穴があって結局3時間かかってしまった(泣)

でも、初めて何にも頼らず自力で解けたというのは大きな進歩だと思う。
これが、赤色への記念すべき歴史的第一歩となることを願うばかりである。。

ヘタレすぎるコードですが、記念ということで貼っておきます。
問題は、こちら。基本的な動的計画なのですが、・・・、同じ高さのところの処理をどうするかでかなり悩みました。解説や赤色の人のコードを参考にスマートな解法を探ってみます。

Single Round Match 404 Round 1 - Division II, Level Three:


class Prb {
public:
int sweet;
int x;
int y;
int len;

Prb(int s, int _x, int _y, int l) {
sweet = s;
x = _x;
y = _y;
len = l;
}
bool operator<(const Prb &prb) const {
if (this->y != prb.y)
return this->y < prb.y;
if (this->x != prb.x)
return this->x < prb.x;
return this->sweet < prb.sweet;
}
};

int dp[64];
bool dn[64];
vector<Prb> pp;

class GetToTheTop {
public:
int collectSweets(int K, vector<int> sweets, vector<int> x, vector<int> y, vector<int> stairLength) {
pp.clear();
memset(dp, 0, sizeof(dp));

pp.PB(Prb(0, 1, 0, 11000));
REP(i, sweets.size())
pp.PB(Prb(sweets[i], x[i], y[i], stairLength[i]));
SORT(pp);
// same y
FOR (i, 1, pp.size()) {
if (pp[i].y != pp[i-1].y)
continue;
int l1 = pp[i].x, r1 = l1+pp[i].len;
int l2 = pp[i-1].x, r2 = l2+pp[i-1].len;

double r;
if (l1 >r2)
r = dist(l1, 0, r2, 0);
else
r = dist(l2, 0, r1, 0);
if (r <= K) {
pp[i].sweet += pp[i-1].sweet;
pp[i-1].sweet = -1;
}
}

FOR (i, 1, pp.size())
if (pp[pp.size()-1-i].sweet == -1)
pp[pp.size()-1-i].sweet = pp[pp.size()-i].sweet;

int posy = 0;
REP(i, pp.size()) {
if (posy != pp[i].y) {
posy = pp[i].y;
meanSameY(i, posy, pp.size(), K);
}
if (i && !dp[i])
continue;
int l1 = pp[i].x;
int r1 = l1 + pp[i].len;

REP(j, pp.size()) {
if (i == j) continue;
if (pp[i].y >= pp[j].y) continue;

int l2 = pp[j].x;
int r2 = l2 + pp[j].len;
double r;
if (r1 < l2) r = dist(r1,pp[i].y, l2, pp[j].y);
else if (l1 > r2) r = dist(l1,pp[i].y, r2, pp[j].y);
else r = dist(0, pp[i].y, 0, pp[j].y);

if (r <= K)
dp[j] = max(dp[j], dp[i] + pp[j].sweet);
}
}

return *max_element(dp, dp+pp.size());
}

double dist(int x1, int y1, int x2, int y2) {
return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
}

void meanSameY(int pos, int h, int sz, int K) {
VVI idx;
VI tmp;
FOR (i, pos, sz) {
if (i >= sz pp[i].y != h) {
idx.PB(tmp);
tmp.clear();
break;
}
if (pp[i].x - (pp[i-1].x + pp[i-1].len) > K) {
idx.PB(tmp);
tmp.clear();
}
tmp.PB(i);
}
REP(i, idx.size()) {
int ret = 0;
tmp = idx[i];
REP(j, tmp.size())
ret = max(ret, dp[tmp[j]]);
REP(j, tmp.size())
dp[tmp[j]] = ret;
}
}
};

2010年11月9日火曜日

ヒープを作ってみよう!

今日は、ヒープを自分で実装してみます。
昔、「二分木は配列で実装できます。」みたいなことを読んだことがありました。
そのときは、??だったが、今ならできる!!

では、早速ソースを。



class myHeap {
public:
    myHeap() {
        this->tail = 0;
    }

    void add(int n) {
        hp[tail++] = n;
        int pos = tail-1;
        while (pos) {
            if (hp[(pos-1)/2] < hp[pos])
                swap(hp[(pos-1)/2], hp[pos]);
            pos = (pos-1)/2;
        }
    }

    int pop() {
        if (isEmpty()) {
            cerr << "Heap is empty!" << endl;
            return -1;
        }

        int ret = hp[0];
        hp[0] = hp[--tail];
        int pos = 0;

        while (2 * pos + 2 <= tail) {
            if (hp[pos] > max(hp[2*pos+1], hp[2*pos+2]))
                break;
            if (hp[2*pos+1] > hp[2*pos+2]) {
                swap(hp[pos], hp[2*pos+1]);
                pos = 2*pos+1;
            }
            else {
                swap(hp[pos], hp[2*pos+2]);
                pos = 2*pos+2;
            }
        }
        return ret;
    }

    bool isEmpty() {
        return tail <= 0;
    }

private:
    int hp[1<<10];
    int tail;
};


で、ちゃんと動くか確認してみましょう。


int main() {
    myHeap hp = myHeap();

    REP(i, 100)
        hp.add(rand() % 1000);

    while (!hp.isEmpty())
        printf("%d\n", hp.pop());

    return 0;
}


はい、動きました。
実際に自分で作ってみると、ヒープのオーダーについての理解も深まります。
  • 最大値の参照は、O(1)
  • push()およびpop()は、O(log n)
なのは明らかですね。

2010年11月7日日曜日

部分和を列挙する

今日は、昨日練習したTopCoderの問題について書きます。

SRM 403 devision2 Hard:

簡単に要約すると、
4と7のみから構成される数字をlucky numberと呼ぶ。
整数nが与えられたときlucky numberのみを足し合わせてnをつくることができる場合は、そのlucky numberを列挙せよ。ただし、そのような組み合わせが複数個ある場合は、個数が最小のもの、さらにそれが複数個ある場合は、辞書順で最小になるものを求めよ。

詳しくは、上のTopCoderのサイトを参照ください。

ここでポイントは、
  1. lucky numberを探索してキャッシュすること
  2. nがlucky numberのみで構成できるかどうかDPを使って判定すること
  3. 構成出来る場合は、その解を問題の定義どおりに列挙すること
2.までは、出来たんだけど・・・

赤い人のコードを参考に復習しました。
まず、1.ですが再帰で書くとかなりスマートに書けます。
2.は部分和問題と同等の方法で書けます。dp[i]には整数iを構成するのに必要なlucky numberの最小値をキャッシュします。(私はdp[i]にiを構成するために必要なluckyu numberをvectorにぶち込んでキャッシュしてました。これだと、メモリ制限や時間制限でoutになりました。。)
3.は2.で作った配列を元に経路復元みたいなことをやります。その際、なるべく小さい数をたくさん取るように(辞書順で小さくなるように)greedyで復元します。

以下、(赤い人を参考に)私が書いたコードを張っておきます。


int dp[1000001];

class TheSumOfLuckyNumbers {
public:
    void makeLucky(vector<int>&A, int x, int n) {
        if (x)
            A.PB(x);
        if (x * 10 + 4 <= n)
            makeLucky(A, x * 10 + 4, n);
        if (x * 10 + 7 <= n)
            makeLucky(A, x * 10 + 7, n);
    }

    vector<int> sum(int n) {
        VI A;
        makeLucky(A, 0, n);
        SORT(A);

        // DP
        EACH(itr, A) {
            dp[*itr] = 1;
            REP(i, n+1) {
                if (!dp[i])
                    continue;
                if (i + *itr > n)
                    continue;
                if (!dp[i + *itr] || dp[i + *itr] >= dp[i] + 1)
                    dp[i + *itr] = dp[i] + 1;
            }
        }

        VI ret;
        int pos = n;
        // Greedy
        EACH(itr, A) {
            while (dp[pos] - dp[pos - *itr] == 1 && pos) {
                if (!dp[pos - *itr] && pos - *itr)
                    break;
                ret.PB(*itr);
                pos -= *itr;
                if (pos - *itr < 0) break;
            }
            if (!pos)
                break;
        }

        return ret;
    }
};



あと、気付きですが、TopCoderの過去問を解くのはPOJを解くよりも良いかもしれません。
その理由は、
  1. 問題のレベル分けがきちんとなされている
  2. 赤い人のコードや解説(editorial)が読める
  3. すべての問題が種類別にカテゴライズされている。
3.のカテゴライズは今日気付きました。
Problem Archiveのところにcategoryってのがあるので、そこで「DP」とか「String manipulation」とか「recursive」とか選択できるみたいです。

2010年11月4日木曜日

知ってると便利なSTL(4) swap

STLを使いこなすことは、アルゴリズムコンテストを勝ち抜くうえでとても重要である。
まだまだ、しらないSTLがたくさんあるようだ。。

今日のSTLはSTL Algorithmsのswap()。
その名のとおり、変数の値をswapできます。

まずは、簡単な例から。

int main() {
int a = 10;
int b = 3;

swap(a, b);
printf("%d %d\n", a, b);

return 0;
}

次に、swapを用いて選択ソートをしてみます。

int main() {
int x[] = {3,4,5,1,2,10,9,7,8,6};

REP(i, 10)
FOR (j, i, 10)
if (x[i] > x[j])
swap(x[i], x[j]);
REP(i, 10)
printf("%d ", x[i]);

return 0;
}
これは、かなり便利そう。
これから、重宝しそうです。