【CQOI2015】选数
题面
Description
我们知道，从区间[L,H]（L和H为整数）中选取N个整数，总共有(H-L+1)^N种方案。小z很好奇这样选出的数的最大公约数的规律，他决定对每种方案选出的N个整数都求一次最大公约数，以便进一步研究。然而他很快发现工作量太大了，于是向你寻求帮助。
你的任务很简单，小z会告诉你一个整数K，你需要回答他最大公约数刚好为K的选取方案有多少个。由于方案数较大，你只需要输出其除以1000000007的余数即可。
Input
输入一行，包含4个空格分开的正整数，依次为N，K，L和H。
Output
输出一个整数，为所求方案数。
Sample Input
2 2 2 4
Sample Output
3
Hint
【样例解释】
所有可能的选择方案：(2, 2), (2, 3), (2, 4), (3, 2), (3, 3), (3, 4), (4, 2), (4, 3), (4, 4)
其中最大公约数等于2的只有3组：(2, 2), (2, 4), (4, 2)
【数据范围】
对于30%的数据，N≤5，H-L≤5
对于100%的数据，1≤N,K≤109，1≤L≤H≤109，H-L≤10^5
题目分析
设\(r=\lfloor\frac HK\rfloor,l=\lfloor\frac {L-1}K\rfloor\)
根据套路：\(\displaystyle ans=\sum_{d=1}^r\mu(d)(\lfloor\frac rd\rfloor-\lfloor\frac ld\rfloor)^N\)
由于\(r\)可能很大，需要用杜教筛处理\(\mu\)的前缀和。
杜教筛：
\[\begin{split}
(g*f)(i)&=\sum_{d|i}g(d)f(\frac id)\\
\Rightarrow g(1)S(n)&=\sum_{i=1}^n(g*f)(i)-\sum_{i=2}^ng(i)S(\frac ni)
\end{split}
\]
其中，\(S(x)\)为\(f()\)的前缀和。
这次，我们的\(f\)为\(\mu\)，根据杜教筛的套路，取\(g(x)=1\)。
\[\begin{split}
S(n)=1-\sum_{i=2}^nS(\frac ni)
\end{split}
\]
可以用线性筛预处理一部分\(\mu\)的前缀和，剩下的用杜教筛记忆化搜索即可。
代码实现
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#include<map>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=1e7+5,M=1e7,mod=1000000007;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int mu[N],prime[N];
bool vis[N];
map<int,int>smu;
int Smu(int x){
if(x<=M)return mu[x];
if(smu[x])return smu[x];
int ret=1;
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
ret-=(r-l+1)*Smu(x/l);
}
return smu[x]=ret;
}
LL ksm(LL x,LL k){
LL ret=1;
while(k){
if(k&1)ret=ret*x%mod;
x=x*x%mod,k>>=1;
}
return ret;
}
int main(){
mu[1]=1;
for(int i=2;i<=M;i++){
if(!vis[i])prime[++prime[0]]=i,mu[i]=-1;
for(int j=1;j<=prime[0]&&i*prime[j]<=M;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0)break;
mu[i*prime[j]]=-mu[i];
}
mu[i]+=mu[i-1];
}
int n=Getint(),K=Getint(),L=(Getint()-1)/K,R=Getint()/K;
int ans=0;
for(int l=1,r;l<=R;l=r+1){
r=R/(R/l);
if(l<=L)r=min(r,L/(L/l));
ans=(ans+1ll*(Smu(r)-Smu(l-1))*ksm(R/l-L/l,n)%mod)%mod;
}
cout<<(ans+mod)%mod;
return 0;
}
