容斥经典：求1-m中与n互素的数的个数（hdu4153）

hdu4153是一个很经典的容斥原理的题目，先贴一下链接：http://acm.hdu.edu.cn/showproblem.php?pid=4135

题意：求区间a-b内与n互素的数的个数

思路：考虑一下逆问题：区间a-b内不与n互素的数的个数，不互素那么一定至少有一个相同的素因子，这样我们就可以枚举n的素因子，来探讨区间a-b中有的n的某个素因子个数，但是并不好计算，所以我们考虑1-(a-1)和1-b两个区间不与n互素的数的个数，那么最后两个数相减就是a-b中不与n互素的个数了。所以就转化成了一个容斥原理里非常经典的问题：求1-m中与n互素的数的个数。

如上所说，如果1-m中某个数与n不互素，那么一定可以被n的某个素因子整除，我们枚举n的素因子。举个例子，比如n=12，m=8。12的素因子是2和3，1-8中有几个数是2的倍数呢，很容易就可以发现是8/2=4，同理1-8中为3倍数的个数为8/3=2，那么1-8中与12不互素的数的个数就是4+2=6个吗？结果显然不是，1-8中与12互素的数是2、3、4、6、8，是五个，这是因为同为2倍数和3倍数的6被计算了两次，所以我们还要再减去一次6的倍数，即8/6=1，所以最后的结果应该是8/2+8/3-8/(2*3)，如果多举几个例子就可以发现存在一个奇偶规律了，如果是奇数个数的组合，那么符号是加，如果是偶数个数的组合，那么符号是减，这便是容斥原理了。

代码采用了二进制的方法，先选出n的素因子，数据范围是1e9，据说素因子个数不会超过10个……最关键的便是如何实现多个素因子的组合，先贴一下模板

vector<LL> p;
LL fun(LL n,LL m)
{
    p.clear();
    for(LL i=2;i*i<=n;i++)
    {
        if(n%i==0)
        {
            p.push_back(i);
            while(n%i==0)
                n/=i;
        }
    }
    if(n>1) p.push_back(n);///求n的素因子
    LL S=0;
    for(LL i=1;i<(1<<p.size());i++)///从1一直枚举到1<<素因子个数
    {
        LL cnt=0;
        LL mul=1;
        for(LL j=0;j<p.size();j++)///枚举每个素因子
        {
            if(i&(1<<j))///若对应位置为1那么计入
            {
                cnt++;
                mul*=p[j];
            }
        }
        if(cnt&1)///容斥原理，选取的个数为奇数个，加
            S+=m/mul;
        else S-=m/mul;///选取的个数为偶数个，减
    }
    return m-S;///S求得的是1-m中与n不互素的数的个数，返回的是1-m中与n互素的数的个数
}

比如枚举的i为1110，那么也就是选取第二个、第三个、第四个素因子的组合，这样就能保证所有的素因子的组合都被计入了，然后再用容斥原理就可以求出1-m中与n互素的数的个数了。

完整代码：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
typedef long long LL;
const int N=1e6+5;
vector<LL> p;
LL fun(LL n,LL m)
{
    p.clear();
    for(LL i=2;i*i<=n;i++)
    {
        if(n%i==0)
        {
            p.push_back(i);
            while(n%i==0)
                n/=i;
        }
    }
    if(n>1) p.push_back(n);
    LL S=0;
    for(LL i=1;i<(1<<p.size());i++)
    {
        LL cnt=0;
        LL mul=1;
        for(LL j=0;j<p.size();j++)
        {
            if(i&(1<<j))
            {
                cnt++;
                mul*=p[j];
            }
        }
        if(cnt&1)
            S+=m/mul;
        else S-=m/mul;
    }
    return m-S;
}
int main()
{
    int t;
    scanf("%d",&t);
    int num=0;
    while(t--)
    {
        num++;
        LL a,b,n;
        scanf("%lld%lld%lld",&a,&b,&n);
        LL ans=fun(n,b)-fun(n,a-1);
        printf("Case #%d: %lld\n",num,ans);
    }
    return 0;
}