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\)
处理完整问题
设\(S(n,i)=\sum\limits_{k=1}^n[LPF_k>p_i]f(k)\),则我们的目标是求解\(S(n,0)\)
考虑如何转移:
- \(p_i*p_i > n\)
- \(p_i*p_i\leq n\)
可以直接按照\(\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\)
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}\)
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\)
#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;
}
标签:prime,right,limits,int,min25,笔记,学习,sum,left
From: https://www.cnblogs.com/clapp/p/17603088.html