[CQOI2015] 选数

前言

这东西好歹得写点题解吧。

题目

[CQOI2015]选数

题目大意：

从 \([L,R]\) 中选 \(N\) 个数(可重复，有序)使得这 \(N\) 个数的最大公因数为 \(K\)，求方案数，对 \(10^9+7\) 取模。

讲解

我们可以先转化题意为从 \((\lfloor\frac{L-1}{K}\rfloor,\lfloor\frac{R}{K}\rfloor]\) 中选 \(N\) 个数使得它们的最大公因数为 \(1\)。为了方便，下文的 \(L,R\) 为转换后的题意。

令 \(f(n)\) 表示从 \((L,R]\) 中选 \(N\) 个数，这 \(N\) 个数的最大公因数为 \(n\) 的方案数。

\(F(n)\) 表示从 \((L,R]\) 中选 \(N\) 个数，这 \(N\) 个数的最大公因数为 \(n\) 的倍数的方案数，那么显然有：

\(F(n)=\sum_{n|d}f(d)=(\lfloor\frac{R}{n}\rfloor-\lfloor\frac{L}{n}\rfloor)^N\)

于是我们直接反演：

\[f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)\\ f(1)=\sum_{i=1}^{R}\mu(i)F(i)\\ f(1)=\sum_{i=1}^{R}\mu(i)(\lfloor\frac{R}{i}\rfloor-\lfloor\frac{L}{i}\rfloor)^N\\ \]

前面 \(\mu\) 的前缀和用杜教筛随便求，后面整除分块即可。

时间复杂度 \(O(R^{\frac{2}{3}})\)，速度一般。

代码

//12252024832524
#include <bits/stdc++.h>
#define TT template<typename T>
using namespace std; 

typedef long long LL;
const int MAXN = 32005;
const int MOD = 1e9 + 7;
int n,k,L,R,ans;

LL Read()
{
	LL x = 0,f = 1;char c = getchar();
	while(c > '9' || c < '0'){if(c == '-')f = -1;c = getchar();}
	while(c >= '0' && c <= '9'){x = (x*10) + (c^48);c = getchar();}
	return x * f;
}
TT void Put1(T x)
{
	if(x > 9) Put1(x/10);
	putchar(x%10^48);
}
TT void Put(T x,char c = -1)
{
	if(x < 0) putchar('-'),x = -x;
	Put1(x); if(c >= 0) putchar(c);
}
TT T Max(T x,T y){return x > y ? x : y;}
TT T Min(T x,T y){return x < y ? x : y;}
TT T Abs(T x){return x < 0 ? -x : x;}

int qpow(int x,int y)
{
	int ret = 1;
	while(y){if(y & 1) ret = 1ll * ret * x % MOD;x = 1ll * x * x % MOD;y >>= 1;}
	return ret;
}
int prime[MAXN],pn,mu[MAXN],pmu[MAXN];
bool vis[MAXN];
void sieve(int x)
{
	pmu[1] = mu[1] = 1;
	for(int i = 2;i <= x;++ i)
	{
		if(!vis[i]) prime[++pn] = i,mu[i] = -1;
		for(int j = 1;j <= pn && i * prime[j] <= x;++ j)
		{
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) break;
			mu[i*prime[j]] = -mu[i];
		}
		pmu[i] = pmu[i-1] + mu[i];
	}
}
unordered_map<int,int> pm;
int premu(int x)
{
	if(x <= Min(31622,R)) return pmu[x];
	if(pm.count(x)) return pm[x];
	int ret = 1;
	for(int l = 2,r;l <= x;l = r+1)
	{
		r = x / (x/l);
		ret -= (r-l+1) * premu(x/l);
	}
	return pm[x] = ret;
}

int main()
{
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	n = Read(); k = Read(); L = (Read()-1) / k; R = Read() / k;
	sieve(Min(31622,R));
	for(int l = 1,r;l <= R;l = r+1)
	{
		if(!(L/l)) r = R/(R/l);
		else r = Min(L/(L/l),R/(R/l));
		ans = (ans + 1ll * (premu(r) - premu(l-1)) * qpow(R/l-L/l,n)) % MOD;
	}
	Put((ans+MOD)%MOD,'\n');
	return 0;
}