牛牛的最大公约数

题目：

https://ac.nowcoder.com/acm/problem/21546
牛牛有一个区间$[L,R]$,需要选择$N$个数,这$N$个数都在这个区间范围内,那么我们知道一共有$(R-L+1{)}^N$种选法,假如我们想要这$N$个数的最大公约数恰好是$K$.请问一共有多少种选法，输出答案对$10^9+7$取模.
$N,K,L,R(1\le N,K,L\le 10^9,L\le R\le 10^9)0\le R-L\le 10^5$

思路：
$$\begin{array}{rl}& \sum _{i_1=l}^r...\sum _{i_n=l}^r[(i_1....i_n)=k]\\ =& \sum _{i_1=\lceil \frac{l}{k}\rceil }^{\lfloor \frac{r}{k}\rfloor }...\sum _{i_n=\lceil \frac{l}{k}\rceil }^{\lfloor \frac{r}{k}\rfloor }[(i_1....i_n)=1]令l=\lceil \frac{l}{k}\rceil ,r=\lfloor \frac{r}{k}\rfloor \\ =& \sum _{i_1=l}^r...\sum _{i_n=l}^r[(i_1....i_n)=1]\\ =& \sum _{i_1=l}^r...\sum _{i_n=l}^r\sum _{d|(i_1...i_2)}\mu (d)\\ =& \sum _{d=1}^r\sum _{i_1=l}^r...\sum _{i_n=l}^r[d|(i_1...i_n)]\mu (d)\\ =& \sum _{d=1}^r\mu (d)\sum _{i_1=l}^r[d|i_1]...\sum _{i_n=l}^r[d|i_n]只要另i_j是d的倍数就满足d|(i_1...i_n)\\ =& \sum _{d=1}^r\mu (d)(\lfloor \frac{r}{d}\rfloor -\lfloor \frac{l-1}{d}\rfloor {)}^n分块+杜教筛\end{array}$$
注意
分块时考虑$d>l-1$的情况，会除$0$

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,k,l,r;
const ll maxn=1e7;//预处理量
const ll mod=1e9+7;
map<ll,ll> mp_mu;
ll vis[maxn+2],p[maxn],mu[maxn+2],tot=0,mu_sum[maxn+2];
void init(){//预处理
	mu[1]=1;
	for(int i=2;i<=maxn;i++){
		if(!vis[i]) p[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*p[j]<=maxn;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){
				mu[i*p[j]]=0;
				break;
			}
			mu[i*p[j]]=-mu[i];
		}
	}
}
ll qpow(ll a,ll n){
	ll ans=1;
	while(n){
		if(n&1) ans=(ans*a)%mod;
		a=(a*a)%mod;
		n>>=1;
	}
	return ans;
}
ll getMu(ll x){
	if(x<=maxn) return mu_sum[x];
	if(mp_mu.find(x)!=mp_mu.end()) return mp_mu[x];
	ll res=1ll;
	for(ll i=2,j;i<=x;i=j+1){
		j=x/(x/i);
		res-=getMu(x/i)*(j-i+1);
	}
	return mp_mu[x]=res;
}
ll du(ll L,ll R){//杜教筛
	ll res=0;
	for(int i=1,j;i<=R;i=j+1){
		if(L-1<i) j=R/(R/i);//防止除零意外
		else j=min(R/(R/i),(L-1)/((L-1)/i));//取小的那个：毕竟是分块嘛
		ll x=getMu(j)-getMu(i-1);
		ll y=qpow(R/i-(L-1)/i,n);
		res=(res+(x*y)%mod)%mod;
	}
	res=(res+mod)%mod;//要是正才行
	return res;
} 
int main(){
	int L,R;
	init();
	scanf("%lld%lld%d%d",&n,&k,&L,&R);
	l=(L+k-1)/k;r=R/k;
	for(int i=1;i<=min(r,maxn);i++) mu_sum[i]=mu_sum[i-1]+mu[i];
	ll ans=du(l,r);
	printf("%lld\n",ans);
}