包除原理とその実装

日曜の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

view plainprint?

1. set<int> primes;  
2.   
3. void fact(int n) {  
4.  FOR (i, 2, n) {  
5.      if (n % i == 0) {  
6.          primes.insert(i);  
7.          n /= i;  
8.          i = 1;  
9.      }  
10.  }  
11.  if (n > 1) primes.insert(n);  
12. }  
13.   
14. int solve(int n, int bit, VI &vec) {  
15.  int on = 0, mul = 1;  
16.   
17.  REP(i, 32)  
18.      if (bit >> i & 1)  
19.          ++on, mul*= vec[i];  
20.   
21.  return (on % 2 ? 1 : -1) * n / mul;  
22. }  
23.   
24. int main() {  
25.  int n;  
26.   
27.  while (cin >> n, n) {  
28.      primes.clear();  
29.      fact(n);  
30.   
31.      int ret = 0;  
32.   
33.      vector<int>vec(ALL(primes));  
34.      REP(bit, 1<<vec.size()) {  
35.          if (!bit) continue;  
36.          ret += solve(n-1, bit, vec);  
37.      }  
38.      cout << n-1-ret << endl;  
39.  }  
40.   
41.  return 0;  
42. }  

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が大きい場合は何らかの工夫が必要なようです。数学的アルゴリズムの世界は奥が深い。。