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
(1leqslant N,Kleqslant10^9,1leqslant Lleqslant Hleqslant10^9,H-Lleqslant10^5)
题目要求（令(H)变为(R)）
[sumlimits_{a_1=L}^Rsumlimits_{a_2=L}^R...sumlimits_{a_n=L}^R[gcdlimits_{i=1}^n(a_i)=K]
]
同样的，我们令所有数整除一个(K)，注意到(Lleqslant a_ileqslant R)，那么再次枚举的(a_i')要满足(Lleqslant a_i'×Kleqslant R)，所以有
[sumlimits_{a_1=lfloorfrac{L-1}{d}
floor+1}^{lfloorfrac{R}{d}
floor}...sumlimits_{a_n=lfloorfrac{L-1}{d}
floor+1}^{lfloorfrac{R}{d}
floor}[gcdlimits_{i=1}^n(a_i)=1]
]
对最后一部分反演，并且将枚举(x)提前，得到
[sumlimits_{x=1}^{lfloorfrac{R}{d}
floor}mu(x)sumlimits_{a_1=lfloorfrac{L-1}{dx}
floor+1}^{lfloorfrac{R}{dx}
floor}...sumlimits_{a_n=lfloorfrac{L-1}{dx}
floor+1}^{lfloorfrac{R}{dx}
floor}1
]
后面那一坨长的都一个样子，而且可以给他们容斥一下，得到
[sumlimits_{x=1}^{lfloorfrac{R}{d}
floor}mu(x)(lfloordfrac{R}{dx}
floor-lfloordfrac{L-1}{dx}
floor)^n
]
如果我们能够求出(sumlimits_{i=1}^nmu(i))，那么我们就可以进行整除分块了
线筛？死了……(nleqslant 10^9)；于是我们杜教筛一波，可以参考浅谈算法——杜教筛
然后这题就做完了
/*program from Wolfycz*/
#include<map>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Fi first
#define Se second
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
static char buf[1000000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
inline int frd(){
int x=0,f=1; char ch=gc();
for (;ch<'0'||ch>'9';ch=gc()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<3)+(x<<1)+ch-'0';
return x*f;
}
inline int read(){
int x=0,f=1; char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+ch-'0';
return x*f;
}
inline void print(int x){
if (x<0) putchar('-'),x=-x;
if (x>9) print(x/10);
putchar(x%10+'0');
}
const int N=1e6,p=1e9+7;
map<int,int>Sum;
int prime[N+10],miu[N+10],sum[N+10];
bool inprime[N+10];
void prepare(){
miu[1]=sum[1]=1; int tot=0;
for (int i=2;i<=N;i++){
if (!inprime[i]) miu[prime[++tot]=i]=-1;
for (int j=1;j<=tot&&i*prime[j]<=N;j++){
inprime[i*prime[j]]=1;
if (i%prime[j]==0){
miu[i*prime[j]]=0;
break;
}miu[i*prime[j]]=-miu[i];
}
sum[i]=sum[i-1]+miu[i];
}
}
int S(int n){
if (n<=N) return sum[n];
map<int,int>::iterator it=Sum.find(n);
if (it!=Sum.end()) return it->Se;
int res=1;
for (int i=2,pos;i<=n;i=pos+1){
pos=n/(n/i);
res=(res-1ll*(pos-i+1)*S(n/i)%p)%p;
}
res=(res+p)%p;
Sum.insert(map<int,int>::value_type(n,res));
return res;
}
int mlt(int a,int b){
int res=1;
for (;b;b>>=1,a=1ll*a*a%p) if (b&1) res=1ll*res*a%p;
return res;
}
int main(){
prepare();
int n=read(),k=read(),L=read(),R=read(),Ans=0;
--L/=k,R/=k;
for (int i=1,pos;i<=R;i=pos+1){
pos=min(L/i?L/(L/i):inf,R/(R/i));
Ans=(Ans+1ll*(S(pos)-S(i-1))*mlt(R/i-L/i,n)%p)%p;
}
printf("%d
",(Ans+p)%p);
return 0;
}