min25筛学习笔记
min25筛
min25筛用于求一类数论函数的前缀和，适用于函数在素数处的取值可以用一个关于此素数的多项式来表示的数论函数。
处理质数部分
这部分我们需要解决\(\sum\limits_{p \subseteq prime}f(p)\)，这里简单起见，假设\(f(p)=p^t\)
用\(s_i\)表示前\(i\)个质数之和，用\(LPF_i\)表示\(i\)的最小质因子，\(p_i\)表示第\(i\)个质数，设$$g(n,i)=\sum\limits_{k=1}^n\left[ LPF_k>p_i \vee k \subseteq prime \right]k^t$$
因此有:$$g(n,0)=\sum\limits_{k=1}nkt$$
假设\(\sqrt{n}\)内有\(cnt\)个质数，那么
\[g(n,cnt)=\sum\limits_{k=1}^n[k \subseteq prime]k^t
\]
因此我们的目标就是求解\(g(n,cnt)\)，接下来考虑转移，即\(g(n,i)\)和\(g(n,i-1)\)的关系。注意到：
\[g(n,i-1)=\sum\limits_{k=1}^n\left[ LPF_k>p_i \vee k \subseteq prime \right]k^t+\sum\limits_{k=1}^n\left[ LPF_k=p_i \wedge k \not\subseteq prime \right]k^t
\]
因此有转移：
- \(p_i*p_i>n\)\[g(n,i)=g(n,i-1) \]
- \(p_i*p_i\leq n\)
\[\begin{aligned}
g(n,i)&=g(n,i-1)-\sum\limits_{k=1}^n\left[ LPF_k=p_i \vee k \not\subseteq prime \right]k^t\\
&=g(n,i-1)-p_i^t\cdot \sum\limits_{k=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}[LPF_k\ge p_i]k^t\\
&=g(n,i-1)-p_i^t\cdot \left[ g\left(\left\lfloor \frac{n}{p_i} \right\rfloor,i-1\right)- s_{i-1}\right]
\end{aligned}
\]
处理完整问题
设\(S(n,i)=\sum\limits_{k=1}^n[LPF_k>p_i]f(k)\)，则我们的目标是求解\(S(n,0)\)
考虑如何转移：
- \(p_i*p_i > n\)
\[S(n,i-1)=g(n,i-1)-s_{i-1}
\]
- \(p_i*p_i\leq n\)
\[\begin{aligned}
S(n,i-1)&= \sum\limits_{p_j^e\leq n\wedge j\ge i}f(p_j^e)\cdot\left([e>1]+\sum\limits_{k=1}^{\left\lfloor \frac{n}{p_j^e}\right\rfloor}[LPF_k>p_j]f(k) \right)+g(n,cnt)-s_{i-1}\\
&=g(n,cnt)-s_{i-1}+\sum\limits_{p_j^e\leq n\wedge j\ge i}f(p_j^e)\cdot \left( S\left(\left\lfloor \frac{n}{p_j^e} \right\rfloor,j\right)+[e>1] \right){\kern 5pt}\left( 1 \right) \\
&=S(n,i)+\sum\limits_{p_i^e\leq n}f(p_i^e)\cdot \left( S\left(\left\lfloor \frac{n}{p_i^e} \right\rfloor,i\right)+1 \right){\kern 86pt}\left( 2 \right)
\end{aligned}
\]
可以直接按照\(\left( 1 \right)\)式暴力求解，也可以按照\(\left( 2 \right)\)式仿照质数部分DP递推求解。
时间复杂度分析
对于质数部分的求解，我们只关注\(t \subseteq \left\{ \left\lfloor \frac{n}{k} \right\rfloor|k \subseteq \mathbb{N}^*\right\}\)处\(g\left( t,i \right)\)的取值，而只有当\(p_i*p_i \leq t\)时才会产生转移。因此复杂度可以估计为：
\[\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{i}}{\ln{\sqrt{i}}}+\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{\frac{n}{i}}}{\ln{\sqrt{\frac{n}{i}}}}=O\left( \int_1^{\sqrt{n}}\frac{\sqrt{\frac{n}{x}}}{\log_2\sqrt{\frac{n}{x}}}dx \right)=O\left( \frac{n^{\frac{3}{4}}}{\log_2{n}} \right)
\]
对于第二部分的复杂度分析类似，值得一提的是，按照\(\left( 1 \right)\)式暴力求解的时间复杂度也是对的，而且这种写法更加简洁。
一些例子
- 求\(\sum\limits_{i=1}^n\sum\limits_{i=1}^nf^k\left( \operatorname{gcd}(i,j) \right)\)其中\(f(x)\)表示\(x\)的次大质因子，重复的质因子计算多次，例如\(f(4)=f(6)=2\)规定\(f(1)=0,f(p)=1\)其中\(p\)为质数。\(n,k\leq 2\times 10^9\)
\[\begin{aligned}
&{\kern 10pt}\sum\limits_{i=1}^n\sum\limits_{i=1}^nf^k\left( \operatorname{gcd}(i,j) \right)\\
&=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\sum\limits_{j=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\left[ \operatorname{gcd}(i,j)=1 \right] f^k(d)\\
&=\sum\limits_{d=1}^nf^k(d)\sum\limits_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\mu(i)\cdot \left\lfloor \frac{n}{d\cdot i} \right\rfloor^2\\
&令g(n)=\sum\limits_{i=1}^n\mu(i)\cdot \left\lfloor \frac{n}{i} \right\rfloor^2，则原式=\sum\limits_{d=1}^nf^k(d)\cdot g\left(\left\lfloor\frac{n}{d} \right\rfloor\right)\\
&于是接下来要做的事情是求g的前缀和以及求f^k的前缀和，g的前缀和可以用数论分块+杜教筛做。\\
&然后考虑求f的前缀和，设S(n,i)=\sum\limits_{j=1}^n\left[ LPF_j>p_i \vee j \subseteq prime \right]f^k(j)\\
&若p_i\cdot p_i>n,S(n,i-1)=\sum\limits_{j=1}^n\left[ j \subseteq prime \right]\\
&若p_i\cdot p_i\leq n,S(n,i-1)=S(n,i)+S\left(\left\lfloor \frac{n}{p_i} \right\rfloor,i-1\right)-\sum\limits_{j=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}\left[ j \subseteq prime \right]+p_i\cdot \left( \left( \sum\limits_{j=1}^{\left\lfloor \frac{n}{p_i} \right\rfloor}[j \subseteq prime] \right)-i+1 \right)\\
&对于n以内质数个数可以直接套用前面min25筛的质数部分
\end{aligned}
\]
code：
#include <cstdio>
#include <iostream>
#define int unsigned int
using namespace std;
const int M = 3e5 + 5;
int prime[M], cnt = 0, a[M], g[M << 1], f[M << 1], primeMi[M];
int n, k, v[M << 1], ss = 0, mu[M], mus[M];
inline int getIdx(int o) { return (o > M) ? (M + n / o) : o; }
int mi(int o, int p) {
int yu = 1;
while (p) {
if (p & 1) yu *= o;
o *= o;
p >>= 1;
}
return yu;
}
int cal(int o) {
if (o < M) return mu[o];
if (mus[n / o]) return mus[n / o];
int yu = 1;
for (int l = 2, r; l <= o; l = r + 1) {
r = o / (o / l);
yu -= (r - l + 1) * cal(o / l);
}
return mus[n / o] = yu;
}
int find(int o) {
int l = 1, r = ss, res = 0;
while (l <= r) {
int mid = (l + r) >> 1;
if (o <= v[mid] / o) {
res = mid;
l = mid + 1;
} else
r = mid - 1;
}
return res;
}
signed main() {
scanf("%u%u", &n, &k);
mu[1] = 1;
for (int i = 2; i < M; i++) {
if (!a[i]) {
prime[++cnt] = i;
primeMi[cnt] = mi(prime[cnt], k);
mu[i] = -1;
}
for (int j = 1; j <= cnt && i * prime[j] < M; j++) {
a[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
} else {
mu[i * prime[j]] = -mu[i];
}
}
}
mu[0] = 0;
for (int i = 1; i < M; i++) mu[i] += mu[i - 1];
for (int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
v[++ss] = n / l;
g[getIdx(n / l)] = n / l - 1;
}
for (int i = 1; i <= cnt; i++) {
for (int j = 1; j <= ss; j++) {
int t = v[j];
if (prime[i] > t / prime[i]) break;
g[getIdx(t)] -= g[getIdx(t / prime[i])] - i + 1;
}
}
for (int j = 1; j <= ss; j++) {
f[getIdx(v[j])] = g[getIdx(v[j])];
}
for (int i = cnt; i >= 1; i--) {
int idx = find(prime[i]);
for (int j = idx; j >= 1; j--) {
int t = v[j];
if (prime[i] > t / prime[i]) continue;
int p1 = getIdx(t);
int p2 = getIdx(t / prime[i]);
f[p1] += f[p2] - g[p2] + primeMi[i] * (g[p2] - i + 1);
}
}
int ans = 0;
f[0] = 0;
for (int l = 1, r; l <= n; l = r + 1) {
int t = n / l, yu = 0;
r = n / t;
for (int l1 = 1, r1; l1 <= t; l1 = r1 + 1) {
r1 = t / (t / l1);
yu += (cal(r1) - cal(l1 - 1)) * (t / l1) * (t / l1);
}
ans += (f[getIdx(r)] - f[getIdx(l - 1)]) * yu;
}
cout << ans << endl;
return 0;
}
- 对于正整数\(n\)，定义如下操作每次选择\(p_i,p_j\)满足\(p_i|n \wedge p_j|n\wedge p_i\neq p_j\)将\(n\)变为\(\frac{n}{p_i\cdot p_j}\)。求有多少n以内的正整数满足可以通过若干次操作变为\(1\)。\(n\leq 10^{11}\)
\[\begin{aligned}
&首先可以通过若干次操作变为1的正整数一定满足质因数的次数之和为偶数且质因数的最大幂次不超过质因数总幂次的一半。\\
&设f(n,j,m,s)=\sum\limits_{i=1}^n\left[ LPF_i>p_j\wedge 当前最大质因子次数为m，质因子此数之和为s 且加入i后(m,s)能到达合法状态\right]\\
&则我们要求f(n,0,0,0)\\
&则有转移f(n,j,m,s)=[在插入一个新质因子后，(m,s)能到达合法状态]\cdot \left( \left( \sum\limits_{i=1}^n[i \subseteq prime] \right) -j\right)+\\
&\sum\limits_{i=j+1}\sum\limits_{a=1}\left( f\left( \left\lfloor \frac{n}{p_i^a} \right\rfloor,i,\max(m,a),s+a \right) +[a>1\wedge 插入p_i^a后(m,s)能到达合法状态] \right)
\end{aligned}
\]
code:
#pragma GCC optimize(2)
#include <chrono>
#include <cstdio>
#include <ctime>
#include <iostream>
using namespace std;
const int M = 317005;
int prime[M], cnt = 0, a[M], ss = 0;
long long g[M << 1], n, v[M << 1];
int getIdx(long long o) { return (o > M) ? (n / o + M) : o; }
int check(int sum, int mx) { return ((mx << 1) <= sum) && ((sum & 1) == 0); }
long long S(long long o, int p, int sum, int mx) {
if (o < prime[p + 1]) return 0;
long long yu = (check(sum + 1, max(mx, 1))) ? (g[getIdx(o)] - p) : 0;
for (int i = p + 1; i <= cnt; i++) {
if (prime[i] > o / prime[i]) break;
int t = 1;
for (long long j = prime[i]; j <= o; j *= prime[i], t++) {
yu += S(o / j, i, sum + t, max(mx, t));
yu += check(sum + t, max(mx, t)) * (t > 1);
}
}
return yu;
}
signed main() {
for (int i = 2; i < M; i++) {
if (!a[i]) prime[++cnt] = i;
for (int j = 1; j <= cnt && i * prime[j] < M; j++) {
a[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
}
}
scanf("%lld", &n);
for (long long l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
g[getIdx(n / l)] = n / l - 1;
v[++ss] = n / l;
}
for (int i = 1; i <= cnt; i++) {
for (int j = 1; j <= ss; j++) {
long long t = v[j];
if (prime[i] > t / prime[i]) break;
g[getIdx(t)] -= g[getIdx(t / prime[i])] - i + 1;
}
}
cout << S(n, 0, 0, 0) << endl;
return 0;
}
- 求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^mf(i\cdot j)\),令积性函数\(f(n)\)满足\(f(p^c)=p^{\operatorname{gcd}(c,k)}\)，其中\(p\)为给定常数，\(p\)为素数，\(c\)为正整数。\(n\leq 10^5,m\leq 10^{10},k\leq 10^9\)
\[\begin{aligned}
&令F(n,x)=\sum\limits_{i=1}^nf(i\cdot x),则ans=\sum\limits_{i=1}^nF(m,i)\\
&我们仿照min25筛的思想，每次提取出一个质数的贡献，假设提取质因子p的贡献，p在x中次数为c\\
&F(n,x)=\sum\limits_{p^e\leq n}f(p^e)\cdot \left( F\left( \left\lfloor \frac{n}{p^e} \right\rfloor,\frac{x}{p^c} \right) -F\left( \left\lfloor \frac{n}{p^{e+1}} \right\rfloor,\frac{x}{p^{c-1}} \right) \right)\\
&这样看起来状态数有O(n\cdot \sqrt{m}),但如果每次都取x的最大质因子，状态数和转移数都在可接受的范围内\\
&极限数据下状态数仅有10^6级别，转移数在10^7次级别。\\
&接下来需要解决的问题是F(n,1)\\
&F(n,1)=\sum\limits_{i=1}^nf(i),注意到f和id在质数处的取值都一样，f的前缀和可以很方便的用Powerful Number解决。\\
&设f=id*h,则S_f(n)=\sum\limits_{d \subseteq PN}h(d)\cdot \frac{\left\lfloor \frac{n}{d} \right\rfloor\cdot (\left\lfloor \frac{n}{d}+1 \right\rfloor)}{2},h(p^c)=p^{\operatorname{gcd}(c,k)}-\sum\limits_{t=1}^cp^t\cdot h(p^{c-t})
\end{aligned}
\]
#include <assert.h>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <vector>
#define int long long
using namespace std;
const int mo = 1e9 + 7;
const int inv2 = (mo + 1) >> 1;
const int M = 1e5 + 5;
const int N = 2e6 + 5;
int n, m, prime[200005], cnt = 0, vis[N + 5], k, tmp = 0, pns[1000005],
cs[N + 5];
vector<int> h[10005];
vector<pair<int, int>> pn;
vector<int> pp[100005];
map<int, int> mp[100005];
inline int gcd(int o, int p) { return (!p) ? o : gcd(p, o % p); }
inline int getIdx(int o) { return (o <= M) ? o : (M + m / o); }
inline int Mod(int o) { return (o >= mo) ? (o - mo) : o; }
int mi(int o, int p) {
int yu = 1;
while (p) {
if (p & 1) yu = yu * o % mo;
o = o * o % mo;
p >>= 1;
}
return yu;
}
vector<array<int, 3>> factor(int o) {
vector<array<int, 3>> ans;
for (int i = 1; i <= cnt; i++) {
if (o % prime[i] == 0) {
int yu = 0;
int tt = 1;
while (o % prime[i] == 0) {
yu++;
tt *= prime[i];
o /= prime[i];
}
ans.push_back(array<int, 3>{yu, prime[i], tt});
}
if (prime[i] * prime[i] >= o) break;
}
if (o > 1) {
ans.push_back(array<int, 3>{1, o, o});
}
return ans;
}
void init() {
memset(pns, -1, sizeof(pns));
memset(cs, -1, sizeof(cs));
for (int i = 2; i <= N; i++) {
if (!vis[i]) prime[++cnt] = i;
for (int j = 1; j <= cnt && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
break;
}
}
}
for (int i = 1; i <= cnt; i++) {
if (prime[i] > 100000) break;
h[i].push_back(1);
for (int j = prime[i], t = 1; j <= m; j *= prime[i], t++) {
int sum = 0;
for (int l = 0; l < t; l++) {
sum = Mod(sum + mi(prime[i], t - l) * h[i][l] % mo);
}
h[i].push_back(Mod(mi(prime[i], gcd(t, k)) - sum + mo));
}
}
}
void pre(int o, int p, int q) {
cs[o] = p % mo;
for (int j = q; j <= cnt; j++) {
if (o * prime[j] > N) break;
for (int t = prime[j] * o, s = 1; t <= N; t *= prime[j], s++) {
pre(t, p * mi(prime[j], gcd(s, k)) % mo, j + 1);
}
}
}
void getPN(int o, int p, int q) {
pn.push_back(make_pair(o, p));
for (int j = q; j <= cnt; j++) {
if (prime[j] > 100000) break;
if (prime[j] * prime[j] > m / o) break;
for (int t = prime[j] * prime[j] * o, s = 2; t <= m;
t *= prime[j], s++) {
getPN(t, p * h[j][s] % mo, j + 1);
}
}
}
inline int s2(int o) { return o * (o + 1) % mo * inv2 % mo; }
int calH(int o) {
int t = getIdx(o);
if (pns[t] != -1) return pns[t];
int ans = 0;
for (auto [l, r] : pn) {
if (l > o) break;
ans = Mod(ans + r * s2(o / l % mo) % mo);
}
return pns[t] = ans;
}
int cal(int o, vector<array<int, 3>> p, int q) {
if (o == 0) return 0;
if (o == 1) return cs[q];
if (p.empty()) return calH(o);
if (o * q <= N) return pp[q][o];
if (mp[q].find(o) != mp[q].end()) return mp[q][o];
tmp++;
int ans = 0;
int q1 = q / p.back()[2];
int q2 = q1 * p.back()[1];
for (int i = 1, t = 0; i <= o; i *= p.back()[1], t++) {
vector<array<int, 3>> p1 = p;
vector<array<int, 3>> p2 = p;
p1.pop_back();
p2.pop_back();
p2.push_back(array<int, 3>{1, p.back()[1], p.back()[1]});
ans =
(ans + mi(p.back()[1], gcd(k, t + p.back()[0])) *
(cal(o / i, p1, q1) - cal(o / i / p.back()[1], p2, q2)) %
mo) %
mo;
ans = Mod(ans + mo);
}
return mp[q][o] = ans;
}
signed main() {
scanf("%lld%lld%lld", &n, &m, &k);
init();
pre(1, 1, 1);
getPN(1, 1, 1);
sort(pn.begin(), pn.end());
for (int i = 1; i <= n; i++) {
pp[i].push_back(0);
for (int j = 1; i * j <= N; j++) {
int yu = Mod(pp[i][pp[i].size() - 1] + cs[i * j]);
pp[i].push_back(yu);
}
}
int ans = 0;
for (int i = 1; i <= n; i++) {
ans = (ans + cal(m, factor(i), i)) % mo;
}
cout << ans << endl;
return 0;
}
