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)
对于100%的数据，1≤N,K≤10^9，1≤L≤H≤10^9，H-L≤10^5
正解：莫比乌斯反演+杜教筛。
看到$PoPoQQQ$设了两个函数，于是照着模仿，自己推了一下。。
设$f(k)$表示$[l,r]$中选数，$gcd=d$的方案数。
设$g(k)$表示$[l,r]$中选数，$d|gcd$的方案数，易知$g(k)=(left lfloor frac{r}{k} ight floor-left lfloor frac{l-1}{k} ight floor)^{n}$。
因为$g(k)=sum_{k|d}f(d)$，由莫比乌斯反演，$f(k)=sum_{k|d} mu(frac{d}{k})g(i)$。
令$d=kQ$，$f(k)=sum_{Q=1}^{left lfloor frac{r}{k} ight floor}mu(Q)(left lfloor frac{r}{kQ} ight floor-left lfloor frac{l-1}{kQ} ight floor)^{n}$。
此时我们可以把$l-1,r$同除以$k$，用杜教筛求出莫比乌斯函数的前缀和，然后数论分块就能计算出答案了。
1 //It is made by wfj_2048~ 2 #include <algorithm> 3 #include <iostream> 4 #include <cstring> 5 #include <cstdlib> 6 #include <cstdio> 7 #include <vector> 8 #include <cmath> 9 #include <queue> 10 #include <stack> 11 #include <map> 12 #include <set> 13 #define rhl (1000000007) 14 #define N (3000010) 15 #define inf (1<<30) 16 #define il inline 17 #define RG register 18 #define ll long long 19 #define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout) 20 21 using namespace std; 22 23 int mu[N],vis[N],prime[N],n,k,l,r,cnt,maxn; 24 ll ans; 25 26 map <int,int> f,vi; 27 28 il int gi(){ 29 RG int x=0,q=1; RG char ch=getchar(); 30 while ((ch<'0' || ch>'9') && ch!='-') ch=getchar(); 31 if (ch=='-') q=-1,ch=getchar(); 32 while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar(); 33 return q*x; 34 } 35 36 il ll qpow(RG ll a,RG ll b){ 37 RG ll ans=1; 38 while (b){ 39 if (b&1) ans=ans*a%rhl; 40 a=a*a%rhl,b>>=1; 41 } 42 return ans; 43 } 44 45 il void sieve(){ 46 mu[1]=1; 47 for (RG int i=2;i<=maxn;++i){ 48 if (!vis[i]) prime[++cnt]=i,mu[i]=rhl-1; 49 for (RG int j=1,k;j<=cnt;++j){ 50 k=i*prime[j]; if (k>maxn) break; vis[k]=1; 51 if (i%prime[j]) mu[k]=rhl-mu[i]; else break; 52 } 53 } 54 for (RG int i=2;i<=maxn;++i){ 55 mu[i]+=mu[i-1]; if (mu[i]>=rhl) mu[i]-=rhl; 56 } 57 return; 58 } 59 60 il ll du(RG int n){ 61 if (n<=maxn) return mu[n]; if (vi[n]) return f[n]; 62 RG ll ans=1; RG int pos; vi[n]=1; 63 for (RG int i=2;i<=n;i=pos+1){ 64 pos=n/(n/i),ans-=(ll)(pos-i+1)*du(n/i)%rhl; 65 if (ans<0) ans+=rhl; 66 } 67 return f[n]=ans; 68 } 69 70 il void work(){ 71 n=gi(),k=gi(),l=(gi()-1)/k,r=gi()/k,maxn=min(3000000,r),sieve(); 72 for (RG int q=1,p,pos;q<=r;q=pos+1){ 73 p=l/q; if (!p) pos=r/(r/q); else pos=min(l/(l/q),r/(r/q)); 74 (ans+=(du(pos)-du(q-1)+rhl)*qpow(r/q-l/q,n))%=rhl; 75 } 76 printf("%lld ",ans); return; 77 } 78 79 int main(){ 80 File("number"); 81 work(); 82 return 0; 83 }