POJ 1091 容斥原理

链接：

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

题意：

给你两个正整数n，m，让你求长度为n+1的满足条件的一个等式:a[1]*x1+a[2]*x2+a[3]*x3+...+a[n]*xn+a[n+1]*x(n+1)=1 (0<=a[i]<=m&&a[n+1]=m)

让你求一共有多少种情况满足这个条件。

要使得 a[1]*x1+a[2]*x2+a[3]*x3+...+a[n]*xn+a[n+1]*m=1 (0<=a[i]<=m)，那么a[1],a[2],a[3]....a[n+1]的最大公约数为1.

题解：

要解决此题，你需要知道的知识有扩展欧几里得，鸽巢原理，以及递归求所有的排列组合。

许多博客都举了这么一个例子：

例如:n=2,m=360  
360=3^2*2^3*5  所有不满足条件的数列，最大公约数是360质因子的乘积，只要将这些组合去掉，就是要求的答案(不懂的慢慢揣摩)

那么就要先求出m的所有质因子，然后求出总的排列组合的个数，即题目中说的M^N，最后根据鸽巢原理求得最后答案。

公式为：ans=M^N-(有奇数个公因数的n元组)+(有偶数个公因数的n元组)。拿上面的例子来说就是

ans=m^n-( 有公因数2的n元组)- (有公因数3的n元组)- (有公因数5的n元组)+ (有公因数2，3的n元组) +(有公因数2，5的n元组)+ (有公因数3，5的n元组)- (有公因数2，3，5的n元组).

有公因数d的n元组，每个位置上有 (m/d)个选择（1 ~ m里面有m/d个d的倍数），根据乘法原理，可以得出有公因数d的n元组有 (m/d)^n 个.

代码：

31 ll factor[100], sz;
32
33 void prime_factor(ll n) {
34     for (ll i = 2; i*i <= n; i++) if (n%i == 0) {
35         factor[sz++] = i;
36         while (n%i == 0) n /= i;
37     }
38     if (n > 1) factor[sz++] = n;
39 }
40
41 ll mod_pow(ll x, ll n) {
42     ll res = 1;
43     while(n) {
44         if (n & 1) res *= x;
45         x *= x;
46         n >>= 1;
47     }
48     return res;
49 }
50
51 ll n, m;
52 ll p[100], sum;
53
54 void dfs(ll id, ll step, ll num) {
55     if (step == num) {
56         ll x = m;
57         rep(i, 0, num) x /= p[i];
58         sum += mod_pow(x, n);
59         return;
60     }
61     rep(i, id, sz) {
62         p[step] = factor[i];
63         dfs(i + 1, step + 1, num);
64     }
65 }
66
67 int main() {
68     cin >> n >> m;
69     prime_factor(m);
70     ll ans = mod_pow(m, n);
71     rep(i, 1, sz + 1) {
72         sum = 0;
73         dfs(0, 0, i);
74         if (i & 1) ans -= sum;
75         else ans += sum;
76     }
77     cout << ans << endl;
78 }