BZOJ3930-莫比乌斯反演+杜教筛

题目的意思很简单，求给定区间内的gcd=k的个数，这应该是传统的莫比乌斯反演了。
有两种思路，一种是直接将里面变成gcd=1，然后里面看作元函数用莫比乌斯函数和恒等函数展开，然后改变求和顺序。
还有一种是构造两个函数，一个是f(x)表示x|gcd的数对个数，一个是g(x)表示x=gcd的数对个数。则f(x)等于g(d)求和，其中x|d，然后再用莫比乌斯反演得到g(x)的表达式，为了缩小求值范围我们可以也将gcd变成1，然后求g(1)。

通过这两种方法都不难得到最后的表达式，即mu(x)*f(x)求和，问题就在于如何求mu(x)上，因为这里x特别的大，我们显然是不能线性求得的，这里就要用到杜教筛。然后筛得的结果用map储存。

详见代码（因为时间不太多，然后写数学公式太复杂所以以后有时间再进行更新）

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<ctime>
#include<climits>
#include<queue>
#include<vector>
#include<set>
#include<map>
using namespace std;typedef long long ll;
const int INF=0x3f3f3f3f;
const int MAXN=1e7+5;
const int mod=1000000007;
int prime[MAXN],mobius[MAXN],sum[MAXN];
bool check[MAXN]; int tot;
map<int,int> Sum;ll quick_pow(ll a,ll b,ll p)
{ll ret=1; a%=p;while(b){if(b&1) ret=ret*a%p; a=a*a%p; b>>=1;}return ret;
}void pre()
{tot=0; mobius[1]=1; sum[1]=1; ll x;for(int i=2;i<MAXN;i++){if(!check[i]){prime[tot++]=i; mobius[i]=-1;}for(int j=0;j<tot && (ll)prime[j]*i<(ll)MAXN;j++){x=prime[j]*i; check[x]=true;if(i%prime[j]) mobius[x]=-mobius[i];else{ mobius[x]=0;break;}}sum[i]=sum[i-1]+mobius[i];}}ll getSum(int x)
{if(x<MAXN) return sum[x];if(Sum[x]) return Sum[x];	//如果已经计算过的直接返回ll ret=1;for(int l=2,r;l<=x;l=r+1){r=x/(x/l); ret-=(r-l+1)*getSum(x/l);//杜教筛，也用分块处理}return Sum[x]=ret;
}int N,K,L,H;int main()
{pre();while(~scanf("%d%d%d%d",&N,&K,&L,&H)){L=(L-1)/K; H=H/K;	//这里要求的gcd为L-H范围内的，通过这样可以得到H/k向下取整和L/K向上取整-1的效果。实际计算的式子是H/k-L/K+1ll ans=0;for(int l=1,r;l<=H;l=r+1){r=H/(H/l);//除法分块if(l<=L) r=min(r,L/(L/l));//当有两个取整函数的时候要注意需要一个一个判断一下，因为有可能把除数变成0，这种错误很隐蔽，需要特别注意ans+=(getSum(r)-getSum(l-1))*quick_pow(H/l-L/l,N,mod); ans%=mod;}printf("%lld\n",(ans+mod)%mod);}return 0;
}