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 个. 

代码：

ll factor[100], sz;

void prime_factor(ll n) {
    for (ll i = 2; i*i <= n; i++) if (n%i == 0) {
        factor[sz++] = i;
        while (n%i == 0) n /= i;
    }
    if (n > 1) factor[sz++] = n;
}

ll mod_pow(ll x, ll n) {
    ll res = 1;
    while(n) {
        if (n & 1) res *= x;
        x *= x;
        n >>= 1;
    }
    return res;
}

ll n, m;
ll p[100], sum;

void dfs(ll id, ll step, ll num) {
    if (step == num) {
        ll x = m;
        rep(i, 0, num) x /= p[i];
        sum += mod_pow(x, n);
        return;
    }
    rep(i, id, sz) {
        p[step] = factor[i];
        dfs(i + 1, step + 1, num);
    }
}

int main() {
    cin >> n >> m;
    prime_factor(m);
    ll ans = mod_pow(m, n);
    rep(i, 1, sz + 1) {
        sum = 0;
        dfs(0, 0, i);
        if (i & 1) ans -= sum;
        else ans += sum;
    }
    cout << ans << endl;
}