Min-25筛
Min-25筛学习笔记
现在写了。
拜谢oi-wiki。
波特好闪。
liuhangshin是我们的红太阳。
Chery写完OJ When？
启动！
问题
其中 \(f\) 是积性函数，\(f(p)\) 是低阶多项式，\(f(p^c)\) 能快速求值。条件其实非常宽松。
约定和记号
- \(p_k\) ：第 \(k\) 小的质数。
- \(mp(n)\) ：\(n\) 的最小质因数。如果\(n = 1\)，定义为 \(1\) 。
- \(F_p(n)\) ：\(\sum_{p \leq n}f(p)\)
- \(F_k(n)\) ：\(\sum_{i = 2} ^ {n} [p_k \leq mp(i)]f(i)\) ，也就是满足最小质因数不小于 \(p_k\) 的数的 \(f\) 值和。
答案
最后加上的那个1是1的贡献，注意这个算法随时都要注意1的问题。
\(part 1\)
根据定义，大力递推。
分别解释一下：
第一步转化是枚举每个数的最小质因子和它出现了几次，那么这部分的贡献根据积性函数的性质是可以分开算的。注意我们只需要枚举 \(\leq\sqrt{n}\) 的质因子。如果最小质因子都大于 \(\sqrt{n}\) 那么这个数肯定是个质数，所以右边我们需要质数处的前缀和。
左边枚举的是每个质数的幂和不是质数幂的情况。一定注意 \(F_k\) 的求和下界是 \(2\) 。
第二步是定义。
第三步很有意思，它直接排除掉了最小质因子满足 \(p_i^c \leq n < p_i^{c + 1}\) 的情况。这是为什么呢？因为如果 \(i\) 乘上它的最小质因子都不能在 \(n\) 以内的话，乘以其它数就更加没有办法保持在 \(n\) 以内了。
这个式子的边界条件是 \(F_k(n)(p_k > n) = 0\) 。
这个看起来非常暴力啊！恭喜你答对了，这个做法的复杂度是 \(\Theta(n^{1-\epsilon})\) ，意思是在 \(n\) 增大的过程中这个玩意的效率越来越接近 \(\Theta(n)\) 。
但是有什么关系呢？反正这玩意跑的飞快。
\(\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\ell\)
手动分割（被打）
\(part 2\)
那么我们要求出 \(F_p(n)\) 。
这个怎么做呢？当然是Meissel-Lehmer 不可能的。
注意到我们求值的所有地方全都是 \(\frac{n}{x}\) 和质数。这个质数的值一定比较小，因为如果质数超过 \(\sqrt{n}\)，在递推进入下一层搜索的 \(\frac{n}{p_i}\) 就小于 \(\sqrt{n}\) ，就会直接返回。所以这个我们直接暴力做前缀和就可以了。
那么我们又知道，所有 \(\frac{n}{x}\) 的取值只有 \(\sqrt{n}\) 个，我们先通过一次整出分块找到所有这样的值并编号，然后把 \(<\sqrt{n}\) 的那部分直接用一个数组保存它的编号，把 \(\geq\sqrt{n}\) 的那部分用另一个数组在 \(\frac{n}{\frac{n}{x}}\) 处保存它的编号。这部分的预处理就做完了。
根据我们的 \(f(p)\) 是低阶多项式，我们拆每一项计算它的贡献。那么我们只需要算 \(\sum_{i = 2}^{n}[isprime(i)]i^x\)，记 \(g(i) = i^x\) 。注意 \(g_i\) 是完全积性函数。
剩下的部分我们 模拟埃筛，记 \(G_k(n) = \sum_{i = 1}^n[isprime(i) \lor mp(i) > p_k]g(i)\) ，也就是埃筛用第 \(k\) 个质数筛完之后剩下函数值的和。
我们写出递推式：
怎么理解？这个式子其实是模拟埃筛用最小的质因子筛掉一个合数的过程。所以首先 \(p_k^2 > n\) 的部分一定不是最小质因子，可以排除。
那么我们先拿所有含有质因子 \(p_k\) 的数出来，那么这些数是 \(p_k,2p_k,3p_k\cdots\) ，在这些数 \(x\) 中，\(\frac{x}{p_k}\) 没有被前 \(k - 1\) 个质数筛掉的部分才具有保留的意义。所以减去的部分是 \(G_{k - 1}(n / p_k)\) 。
但是如果 \(\frac{x}{p_k}\) 就是前 \(k - 1\) 个质数之一呢？这里也会被减掉。所以我们需要加上 \(G_{k - 1}(p_{k} - 1)\) ，注意这也就是前 \(k - 1\) 个质数处的 \(g\) 值。因为 \(p_k^2 \leq n\) ，所以这个也是可以预处理的。
发现这里访问到的 \(G\) 的取值也只会是 \(\frac{n}{x}\)直接按照这个dp，复杂度 \(\Theta(\frac{n^\frac{3}{4}}{\log(n)})\) 。显然可以滚动数组。
复杂度证明先咕了，反正 oi-wiki 上也有。
这下应该是真没了。
上代码。
code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int P = 1e9+7, inv2 = (P + 1) >> 1;
const int N = 1e6+5;
template<typename T>inline void read(T &a)
{
a = 0;
int f = 1;
char c = getchar();
while(!isdigit(c))
{
if(c == '-') f = -1;
c = getchar();
}
while(isdigit(c))
a = a * 10 + c - 48,c = getchar();
a *= f;
}
template<typename T,typename ...L> inline void read(T &a,L &...l)
{
read(a),read(l...);
}
inline int Mul(ll a, ll b) {
return (a % P) * (b % P) % P;
}
inline int Add(int a, int b) {
return a + b >= P ? a + b - P : a + b;
}
inline int Sub(int a, int b) {
return a - b < 0 ? a - b + P : a - b;
}
inline void Inc(int &a, int b) {
a = Add(a, b);
}
const int inv6 = Mul((P + 1) >> 1, (P + 1) / 3);
vector<int> primes;
bool isprime[N];
ll n, val[N];
int le[N], ge[N], ptr, g1[N], g2[N];
int h1[N], h2[N];//质数的值的前缀和
ll fs(ll x) {
return Mul(x % P, x % P - 1);
}
void xxs() {
for(int i = 2; i < N; i++) {
h1[i] = h1[i - 1], h2[i] = h2[i - 1];
if(isprime[i]) {
primes.emplace_back(i);
Inc(h1[i], i);
Inc(h2[i], Mul(i, i));
}
for(auto it:primes) {
if(i * it > N) break;
isprime[i * it] = 0;
}
}
}
int F(int k, ll n) {
if(primes[k] > n) return 0;
int res = 0;
for(int i = k; i < int(primes.size()); i++) {
if(primes[i] > n / primes[i]) break;
ll d = primes[i];
while(d <= n / primes[i]) {
Inc(res, Add(Mul(fs(d), F(i + 1, n / d)), fs(d * primes[i]))), d *= primes[i];
}
}
int id = 0;
if(n > N) {
id = ge[::n / n];
}
else {
id = le[n];
}
return Add(res, Sub(Sub(g2[id], g1[id]), Sub(h2[primes[k] - 1], h1[primes[k] - 1])));
}
int main() {
read(n);
memset(isprime, 1, sizeof isprime);
isprime[1] = 0;
xxs();
for(ll l = 1; l <= n; l = (n / (n / l)) + 1) {
ll d = (n / l);
val[++ptr] = d;
if(d < N) {
le[d] = ptr;
}
else {
ge[n / d] = ptr;
}
g1[ptr] = Mul(d - 1, Mul(d + 2, inv2));
g2[ptr] = Sub(Mul(Mul(Mul(d, d + 1), 2 * d + 1), inv6), 1);//减掉1的值
}
for(int i = 1; i < N; i++) {
h1[i] = h1[i - 1], h2[i] = h2[i - 1];
if(isprime[i]) {
Inc(h1[i], i);
Inc(h2[i], Mul(i, i));
}
}
for(auto it:primes) {
for(int i = 1; i <= ptr; i++) {
if(it * ll(it) > val[i]) break;
ll c = val[i] / it;
int pos = 0;
if(c >= N) {
pos = ge[n / c];
}
else {
pos = le[c];
}
g1[i] = Sub(g1[i], Mul(it, g1[pos]));
g1[i] = Add(g1[i], Mul(it, h1[it - 1]));
g2[i] = Sub(g2[i], Mul(Mul(it, it), g2[pos]));
g2[i] = Add(g2[i], Mul(Mul(it, it), h2[it - 1]));
}
}
printf("%d\n", Add(F(0, n), 1));
}
/*
start coding at:2023/12/18 22:09
finish debugging at:2023/12/18 19:36
stubid mistakes:d没有累乘，算G时fs的参数是it，元素在表中是倒序的，预处理g和h的时候使用定义式，少开ll
*/
/*
吾日三省吾身：
你排序了吗？
你取模了吗？
你用%lld输出long long 了吗？
1LL<<x写对了吗？
判断=0用了abs()吗？
算组合数判了a<b吗？
线段树build()了吗？
.size()加了(signed)吗？
树链剖分的DFS序是在第二次更新的吗？
修改在询问前面吗？
线段树合并到叶子结点return了吗？
__builtin_ctz后面需要加ll吗？
*/
