[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
\(1\leqslant N,K\leqslant10^9,1\leqslant L\leqslant H\leqslant10^9,H-L\leqslant10^5\)

---

题目要求（令\(H\)变为\(R\)）

\[\sum\limits_{a_1=L}^R\sum\limits_{a_2=L}^R...\sum\limits_{a_n=L}^R[\gcd\limits_{i=1}^n(a_i)=K] \]

同样的，我们令所有数整除一个\(K\)，注意到\(L\leqslant a_i\leqslant R\)，那么再次枚举的\(a_i'\)要满足\(L\leqslant a_i'×K\leqslant R\)，所以有

\[\sum\limits_{a_1=\lfloor\frac{L-1}{d}\rfloor+1}^{\lfloor\frac{R}{d}\rfloor}...\sum\limits_{a_n=\lfloor\frac{L-1}{d}\rfloor+1}^{\lfloor\frac{R}{d}\rfloor}[\gcd\limits_{i=1}^n(a_i)=1] \]

对最后一部分反演，并且将枚举\(x\)提前，得到

\[\sum\limits_{x=1}^{\lfloor\frac{R}{d}\rfloor}\mu(x)\sum\limits_{a_1=\lfloor\frac{L-1}{dx}\rfloor+1}^{\lfloor\frac{R}{dx}\rfloor}...\sum\limits_{a_n=\lfloor\frac{L-1}{dx}\rfloor+1}^{\lfloor\frac{R}{dx}\rfloor}1 \]

后面那一坨长的都一个样子，而且可以给他们容斥一下，得到

\[\sum\limits_{x=1}^{\lfloor\frac{R}{d}\rfloor}\mu(x)(\lfloor\dfrac{R}{dx}\rfloor-\lfloor\dfrac{L-1}{dx}\rfloor)^n \]

如果我们能够求出\(\sum\limits_{i=1}^n\mu(i)\)，那么我们就可以进行整除分块了

线筛？死了……\(n\leqslant 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\n",(Ans+p)%p);
    return 0;
}