Page List

Search on the blog

2011年5月27日金曜日

ビット演算魔術、部分集合を高速に出すの巻

 以下のビット演算が何を意味するのか知ってる人は読まなくていいです。
--y &= x

 次の問題を解いているときに出会いました。
[SRM474 SquaresCovering]
2次元平面上に分散するN個の点を正方形で多いたい。正方形はK種類あり、それぞれ大きさと使用するためのコストが異なる。すべての点を覆うには最低いくらコストが必要か?ただし、正方形は無限個あり、同じ正方形を複数回使用しても構わない。

N<=16, K<=50
0 < xi < 10^9, i=1,..,n
0 < yi < 10^9, i=1,..,n

 ぱっと見普通のビットDPで解けそう。が、2秒縛りがキツイ。straightforwardに解くとO(2^{2N})なので、なんとかしたい。
何とかできなかった。。
以下rngさんの回答。やっぱ数学オリンピック金メダリストは違う。。と思ったけど、アリーナからコピペ出来なかったので、私のヘタレソースで我慢してください。(横にはみ出していて、見にくい場合はview plainで見てください。)
const long long int INF = 999999999999999LL;
long long int dp[1<<16];

class SquaresCovering {
public:
int minCost(vector<int> x, vector<int> y, vector<int> cost, vector<int> sides) {
int N = x.size();
int K = sides.size();

REP(mask, 1<<N) {
LL xmin = INF, xmax = -1;
LL ymin = INF, ymax = -1;

REP(i, N) {
if (mask >> i & 1) {
xmin = MIN(x[i], xmin);
xmax = MAX(x[i], xmax);
ymin = MIN(y[i], ymin);
ymax = MAX(y[i], ymax);
}
}

LL val = INF;
int diff = max(xmax-xmin, ymax-ymin);
REP(i, K)
if (diff <= sides[i])
val = MIN(val, cost[i]);
dp[mask] = val;
}

REP(x, 1<<N) {
for (int y = x; y > 0; --y &= x)
dp[x] = min(dp[x], dp[y]+dp[x^y]);
}

return (int)dp[(1<<N)-1];
}
};

 --y & xは、集合xの部分集合をすべて列挙することができます。これはかなりカッコいいです。ピンと来ない人は次のソースを実行してみましょう。
#include <bitset>
int main() {
int x = 202;

for (int y = x; y > 0; --y &= x)
cout << bitset<8>(y) << endl;
}
上の実行結果を見れば何がしたいか分かるはずです。

 最後に、比較してみましょう。SRMの問題を愚直な2重ループでやると、4*10^9くらいかかります。
この方法を使用すれば、2^16 * 2^8 = 2^24 ~16*10^6くらいには収まるかなという予想です。
int main() {
LL x = 0;
REP(mask1, 1<<16)
REP(mask2, 1<<16)
x++;

LL y = 0;
REP(mask1, 1<<16)
for (int mask2 = mask1; mask2 > 0; --mask2 &= mask1)
y++;

cout << x << endl;
cout << y << endl;

return 0;
}
上のプログラムを実行した場合、出力される数値は、
4294967296
42981185
です。まあこんなもんでしょうか。100倍くらい速くなります。



0 件のコメント:

コメントを投稿