Description：

求从 $[L,R]$ 中选 $n$ 个可以相同的数字，他们的最大公约数为 $k$ 的方案数。

$1\leqslant n,k,L,R\leqslant 10^9,R-L\leqslant 10^5$

Solution：

经典的莫比乌斯反演： $$ \begin{align} &\sum_{i_1=L}^R\sum_{i_2=L}^R\cdots\sum_{i_n=L}^R[\gcd(i_1,i_2,\dots,i_n)=k]\\ =&\sum_{i_1=l}^r\sum_{i_2=l}^r\cdots\sum_{i_n=l}^r\sum_{d|i_1,i_2,\dots,i_n}\mu(d)(l=(L-1)/k+1,r=R/k)\\ =&\sum_{d=1}^r\mu(d)(r'-l')^n(r'=r/d,l'=(l-1)/d) \end{align} $$ 右边可以除法分块，求 $\mu$ 的前缀和可以用杜教筛，注意 $[1,l]$ 和 $[l+1,r]$ 应该分开求，因为 $L/l$ 为 $0$ 时会出错。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll n,k,L,R;

#define MAXN 1000010

bool isprime[MAXN];

int prime[MAXN],tot = 0;

ll mu[MAXN];

ll sum[MAXN];

map<int,ll> mu_;

#define MOD 1000000007

ll calc(int k)

{

	if(k < MAXN)return sum[k];

	if(mu_.find(k) != mu_.end())return mu_[k];

	if(k == 1)return 1;

	if(k == 0)return 0;

	ll res = 1;

	for(int l = 2,r;l <= k;l = r + 1)

	{

		r = k / (k / l);

		res = (res - (r - l + 1) * calc(k / l) % MOD + MOD) % MOD;

	}

	mu_[k] = res;

	return res;

}

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % MOD;

		a = a * a % MOD;

		b = b >> 1;

	}

	return res;

}

int main()

{

	scanf("%lld%lld%lld%lld",&n,&k,&L,&R);

	L = (L - 1) / k + 1;R = R / k;

	for(int i = 2;i < MAXN;++i)isprime[i] = true;

	mu[1] = sum[1] = 1;

	for(int i = 2;i < MAXN;++i)

	{

		if(isprime[i])mu[i] = -1,prime[++tot] = i;

		for(int j = 1;j <= tot && i * prime[j] < MAXN;++j)

		{

			isprime[i * prime[j]] = false;

			if(i % prime[j] == 0){mu[i * prime[j]] = 0;break;}

			else mu[i * prime[j]] = -1 * mu[i];

		}

		sum[i] = sum[i - 1] + mu[i];

	}

	ll ans = 0;

	--L;

	for(int l = 1,r;l <= L;l = r + 1)

	{

		r = min(R / (R / l),L / (L / l));

		ans = (ans + (calc(r) - calc(l - 1)) * power(R / l - L / l,n) % MOD + MOD) % MOD;

	}

	for(int l = L + 1,r;l <= R;l = r + 1)

	{

		r = R / (R / l);

		ans = (ans + (calc(r) - calc(l - 1)) * power(R / l,n) % MOD + MOD) % MOD;

	}

	cout << ans << endl;

	return 0;

}