包除原理とその実装

日曜のSRM、no-contestとは。。悲しすぎる。。

本番は解けなかったHardの問題。包除原理を使うことは分かっていたけど、実装が間に合わなかった。

今日は、簡単な包除原理の復習。英語では、Inclusion-exclusion principle.

そのまんま。

この原理を用いると、集合A1, A2, A3, ... Anがあるときに、すべての集合の和集合の要素数を求めることができます。

一番簡単な例で、A1とA2の2つしか集合が無い場合を考えます。これは簡単。

|A1| + |A2| - |A1 and A2|

ですね。ここで、|x|は集合xの要素数です。

じゃ、A1, A2, A3と集合が3つになったら？？

これはベン図を書いて、少し考えると分かります。

|A1| + |A2| + |A3| - |A1 and A2| - |A1 and A3| - |A2 and A3| + |A1 and A2 and A3|

です。

じゃあ、集合がn個あるときは？を説明したのがInclusion-exclusion principleです。

いろいろ解説サイトがありますが、wikiの説明が1番分かりやすいです。

http://en.wikipedia.org/wiki/Inclusion%E2%80%93exclusion_principle

証明も載ってます。

それでは、実装してみましょう。POJの問題を解いてみます。

ぱっと見簡単そうですが、入力値が10^9のときも時間内で解くためには、工夫が必要です。包除原理を使って解いてみます。

http://poj.org/problem?id=2407

set<int> primes;

void fact(int n) {
 FOR (i, 2, n) {
     if (n % i == 0) {
         primes.insert(i);
         n /= i;
         i = 1;
     }
 }
 if (n > 1) primes.insert(n);
}

int solve(int n, int bit, VI &vec) {
 int on = 0, mul = 1;

 REP(i, 32)
     if (bit >> i & 1)
         ++on, mul*= vec[i];

 return (on % 2 ? 1 : -1) * n / mul;
}

int main() {
 int n;

 while (cin >> n, n) {
     primes.clear();
     fact(n);

     int ret = 0;

     vector<int>vec(ALL(primes));
     REP(bit, 1<<vec.size()) {
         if (!bit) continue;
         ret += solve(n-1, bit, vec);
     }
     cout << n-1-ret << endl;
 }

 return 0;
}

部分集合をビット演算を使って実装するのがポイントでしょうか～。

2^n - 1パターンの部分集合を考えないといけないので、nが大きい場合は何らかの工夫が必要なようです。数学的アルゴリズムの世界は奥が深い。。