当前位置：首页 > news >正文

【题解】P3172 [CQOI2015] 选数（倍数莫反做法）

news 2025/8/6 7:22:16

题面：从 $[L,H]$ 选 $N$ 个数（可能重复），一共有 $(H-L+1)^N$ 种取法。问这些取法中满足 $N$ 个数最大公约数为 $K$ 的有多少种。

第一眼就知道要把 $L=\left \lfloor \frac{L}{K} \right \rfloor+[L\,mod\,K\neq 0]$ ， $H=\left \lfloor \frac{H}{K} \right \rfloor$ 。

然后问题就变成求满足 $N$ 个数最大公约数为 $1$ 有多少种。

直接来解，复习下莫反公式：

倍数形式： $F(n)=\sum_{n|d}^{}f(d) \Rightarrow f(n) = \sum_{n|d}^{}\mu (\frac{d}{n})F(d)$ ，不会点这里。

考虑构造 $F(n)=\sum_{n|d}^{}f(d)$ ，结合题目则有两个函数的定义：

$F(n)$ ：从 $[L,H]$ 中选 $N$ 个数的最大公约数是 $n$ 或者 $n$ 的倍数的情况数。

$f(n)$ ：从 $[L,H]$ 中选 $N$ 个数的最大公约数是 $n$ 的情况数。

那么就有： $f(n)=\sum_{n|d}^{}\mu (\frac{d}{n})F(d)$ 。

我们只需要最大公约数为 $1$ 的情况，所以：

$f(1)=\sum_{d=1}^{H}\mu (d)F(d)$

（其实这里的 $d$ 严格意义上是无限大的，但本题最大值只能取到 $H$ ）

现在我们来看看怎么求 $F(d)$ ：

回顾定义，发现应该是： $(\left \lfloor \frac{H}{d} \right \rfloor-(\left \lfloor \frac{L}{d} \right \rfloor+[L\,mod\,d\neq 0])+1)^N$

很好，差一点就能分块加速了。。

其实我们可以让 $L-1$ ，这样直接 $(\left \lfloor \frac{H}{d} \right \rfloor-\left \lfloor \frac{L-1}{d} \right \rfloor)^N$

（当 $n=1$ 的情况也符合，大家可以自己试下）

那最后就变成了： $f(1)=\sum_{d=1}^{H}\mu (d)(\left \lfloor \frac{H}{d} \right \rfloor-\left \lfloor \frac{L-1}{d} \right \rfloor)^N$

代码：

#include<bits/stdc++.h>
using namespace std;
template<typename T> void qread(T &x){x=0; int f=1; char c=getchar();for(; !isdigit(c); c=getchar()) if(c=='-') f=-1;for(; isdigit(c); c=getchar()) x=x*10+(c-'0');x*=f;
}
typedef long long LL;
const int N=5e6+10;
const LL P=1e9+7;
int pr, prime[N];
LL mu[N];
bool v[N]; 
void init(){pr=0; memset(v, 0, sizeof(v));mu[1]=1; mu[0]=0;for(int i=2; i<=N-10; i++){if(!v[i]){pr++; prime[pr]=i;mu[i]=-1;} for(int j=1; (j<=pr) && (i*prime[j]<=N-10); j++){int p=prime[j];v[i*p]=1;if(i%p==0){mu[i*p]=0;break;}else mu[i*p]=-mu[i]; }}for(int i=1; i<=N-10; i++) mu[i]+=mu[i-1];
}
LL q_pow(LL a, LL b){LL res=1;while(b){if(b&1) res=res*a%P;a=a*a%P; b/=2;}return res;
}
LL n, K, L, H;
map<LL, LL> hsm;
LL calc(LL x){if(x<=N-10) return mu[x];if(hsm.find(x)!=hsm.end()) return hsm[x];LL res=1;for(int i=2, j; i<=x; i=j+1){j=x/(x/i);res=(res-calc(x/i)*(j-i+1)%P+P)%P;}return hsm[x]=res;
}
LL solve(){LL res=0; L--;for(int i=1, j; i<=H; i=j+1){if(!(L/i)) j=H/(H/i); //到了一定位置大家的 L/i都变成 0，没有参考价值了 else j=min(L/(L/i), H/(H/i));res=(res+q_pow(H/i-L/i, n)*((calc(j)-calc(i-1))%P+P)%P)%P;}return res;
}
int main(){init();qread(n); qread(K); qread(L); qread(H);H=H/K; L=((L%K==0)? L=L/K: L=L/K+1);printf("%lld\n", solve());return 0;
}

查看全文

http://www.lryc.cn/news/610557.html