HDU 5768 Lucky7（中国剩余定理+容斥原理）

Description
给出区间[l,r]，问[l,r]中能整除7且模pi不等于ai的数的个数
Input
第一行一整数T表示用例组数，每组用例首先输入三个整数n,l,r,之后n行每行两个整数pi,ai
(T<=20,0<=n<=15,0< l< r< 10^18,pi互不相同且pi!=7,p1*…*pn<=10^18,0<=ai< pi<=10^5)
Output
对于每组用例，输出[l,r]中能够整除7且模pi不等于ai的数的个数
Sample Input
2
2 1 100
3 2
5 3
0 1 100
Sample Output
Case #1: 7
Case #2: 14
Solution
令S为[1,N]中能够整除7的数的集合，A为[1,N]中能够整除7且模pi不等于ai的数的集合，Ai为[1,N]中能够整数7且模pi等于ai的数的集合，由容斥定理有

后面若干项Ai的并每一项都可以用一次中国剩余定理解决，注意中间有一处乘法可能会爆longlong所以要用快速乘
Code

#include<cstdio>#include<iostream>using namespace std;typedef long long ll;typedef unsigned long long ull;#define maxn 25ll a[maxn],p[maxn],c[maxn],m[maxn];ll extend_gcd(ll a,ll b,ll &x,ll &y)//扩展欧几里得,求ax+by=gcd(a,b) {     ll d=a;    if(b!=0)    {        d=extend_gcd(b,a%b,y,x);        y-=(a/b)*x;    }    else    {        x=1;        y=0;    }    return d;}ll mul(ll a,ll b,ll p)//快速乘,求a*b(mod p) {    a=(a%p+p)%p,b=(b%p+p)%p;    ll ans=0;    if(a<b)swap(a,b);    while(b)    {        if(b&1)ans=(ans+a)%p;        a=(a+a)%p;        b>>=1;    }    return ans;}ll CRT(int n,ll N)//a数组表示余数,m数组表示模数,n表示方程个数,N表示解的上限 {    ll M=1,Mi,ans=0;    ll x,y,d;    for(int i=0;i<n;i++) M*=m[i];    for(int i=0;i<n;i++)    {        Mi=M/m[i];        d=extend_gcd(Mi,m[i],x,y);        ans=(ans+mul(mul(Mi,x,M),a[i],M))%M;    }    if(ans<0) ans=(ans%M+M)%M;    if(ans>N)return 0;    return (N-ans)/M+1;}ll Solve(int n,ll N){    if(N==0)return 0;    ll ans=N/7;    for(int i=1;i<(1<<n);i++)    {        int res=0;        for(int j=0;j<n;j++)            if(i&(1<<j))                a[res]=c[j],m[res++]=p[j];        a[res]=0,m[res++]=7;            ll cnt=CRT(res,N);        if(res%2)ans+=cnt;        else ans-=cnt;    }    return ans;}int main(){    int T,n,Case=1;    ll x,y;    scanf("%d",&T);    while(T--)    {        scanf("%d%I64d%I64d",&n,&x,&y);        for(int i=0;i<n;i++)            scanf("%I64d%I64d",&p[i],&c[i]);            ll ans=Solve(n,y)-Solve(n,x-1);        printf("Case #%d: %I64d\n",Case++,ans);    }    return 0;}