min_25筛法-亚线性时间内求解积性函数的前缀和
min_25筛
参考链接：
https://www.luogu.com.cn/blog/wucstdio/solution-p5325
https://www.cnblogs.com/zkyJuruo/p/13445596.html#处理整个函数的前缀和
wucstdio同学讲的很清楚了，但是中间递推g的式子错了，应该是\(p_j^k(g(\lfloor\frac{n}{p_j}\rfloor, j-1)-g(p_{j-1}, j-1))\).
min_25筛可以用来在亚线性时间内求解满足条件的积性函数的前缀和。
条件为： 当\(p\)为质数时，\(f(p)\)是一个关于\(p\)的项数较少的多项式，或可以快速求值；\(f(p^c)\)可以快速求值
设\(S(n,j)=\sum_{i=2}^{n}f(i)\times[i的最小质因子比第j个质数大]\)，那么我们要求的就是\(f(1)+S(n,0)\)；其中\(f\)是满足条件的函数。
\(S\)的求解，被分为了\(i\)为质数和\(i\)为合数俩部分。
1.质数部分的前缀和
假设完全积性函数\(f'\)在质数处取值与\(f\)相同。
定义\(g(n,j)\)为\([1,n]\)中用前\(j\)个质数执行埃拉斯托尼筛法后，还未被筛去的数对答案的贡献。
即
那么质数部分对答案的贡献
其中\(P\)表示\([1,n]\)中全体质数的集合。
\(g(n,0)=\sum_{i=2}^{n}f'(i)\)
\(g(n,|P|)\)和\(f'\)在合数处的取值无关，故在中间递推时可以用\(f'\)代替\(f\)。
设\(p_j\)为第\(j\)个质数。
当\(p_j^2>n\)时，显然有\(g(n,j)=g(n,j-1)\) (没有以\(p_j\)作为最小质因数的合数能够被筛去)
当\(p_j^2<=n\)时，在从\(g(n,j-1)\)向\(g(n,j)\)的转移中新被筛去的数，是恰好以\(p_j\)作为最小质因数的合数。
利用完全积性函数的性质把它提取出来，得到
提取了一个\(p_j\)后，剩下的数中合数的最小质因子都大于等于\(p_j\)。
但\(g(\lfloor\frac{n}{p_j}\rfloor,j-1)\)还包括了前\(j-1\)个质数对答案的贡献，应该减去。
最终得到了上述转移方程，能够用于计算\(g(n,|P|)\)
我们还需要处理出所有的\(\lfloor\frac{i}{p_j}\rfloor\)，这需要利用到以下性质：
这意味着我们不需要求出转移过程中所有的\(n\)，只需要根据整除分块的那套理论对\(n\)求出所有\(\lfloor\frac{n}{x}\rfloor\)即可。
根据整除分块的那套理论，对于\(i\)，满足\(\lfloor\frac{n}{i}\rfloor=k\)的所有的整数\(i\)所在区间的左端点记为\(l\)的话，右端点就是\(\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor=\lfloor\frac{n}{k}\rfloor\)，记为\(r\)
证明：
先证当\(i\)取\(\lfloor\frac{n}{k}\rfloor\)时满足\(\lfloor\frac{n}{i}\rfloor=k\)，
因为\(\lfloor\frac{n}{k}\rfloor\le\frac{n}{k}\)
故\(\lfloor\frac{n}{\lfloor\frac{n}{k}\rfloor}\rfloor\ge k\)
再证\(r\ge i\):
再证\(r=\lfloor\frac{n}{k}\rfloor\)
又因为\(i\)是正整数，故\(r=\lfloor\frac{n}{k}\rfloor\)
转移过程中的数都是\(\lfloor\frac{n}{x}\rfloor\)的形式，数量大概是\(O(\sqrt n)\)级别。故我们可以手动离散化一下，用\(idx1\)记录\(\lfloor\frac{n}{x}\rfloor\le \sqrt{n}\)时的\(\lfloor\frac{n}{x}\rfloor\)的下标，\(idx2\)记录另一种情况时的下标.
2.合数部分的前缀和
枚举每一个数的最小质因数及其次数，由于这里满足互质的条件，可以直接把这个质因数的幂提取出来。
那么合数部分的前缀和为
\([e>1]\)表示\(p_k^e\)自身对答案的贡献(\(S\)没有计算f(1)的贡献，并且质数部分已经计算过了)
那么
那么最终答案是\(S(n,0)+f(1)\)
即
实操：洛谷P5325 【模板】Min_25筛
符合条件的函数为\(f(x)\)，且\(f(p^k)=p^k(p^k-1)\)
把多项式\(f(p^k)\)表示为若干个单项式的和，得
依次取\(f'(x)=x^2,f'(x)=x\).
在求解\(g\)时，我们要计算\(f'\)在小于等于\(\sqrt n\)的质数处的前缀和\(g(p_{j-1},j-1)\)，分别将一次幂和二次幂的前缀和记为sum1和sum2，这个可以在线性筛的时候顺便完成。
然后用整除分块求出所有\(\lfloor\frac{n}{x}\rfloor = k\)，并离散化存一下.
对于所有的\(k\)， 我们需要求出\([1,k]\)的前缀和。
把\(f'\)拆分成一次项和二次项，分别求前缀和即可。但这个前缀和不能递推，需要用等差数列求和公式以及二次幂求和公式快速计算出\(k\)处的前缀和。
求解\(g\)的时候要注意，需要保证\(f'\)是完全积性函数，也就是说，必须分别求出一次幂作为\(f'\)和二次幂作为\(f'\)时的\(g\)，不妨分别设为\(g_1,g_2\)
注意\(l，r\)得为\(long long\)
首先整除分块并求出\(g(i, 0)\)
for (ll l = 1, r; l <= n; l = r + 1) {
ll k = n / l;
r = n / k;
w[++totw] = k;
//求g(i,0)
//f'(i)=i^2-i
//g(i,0) = \sum_{j=2}^{n}f'(j) = 一次幂的前缀和减去二次幂的前缀和减去f'(1)
//这题f'(1)是0 不用减了。。。
ll tk = k % md;
g2[totw] = tk * (tk+1) % md * inv2 % md * (2*tk+1) % md * inv3 % md;
g2[totw]--;
g1[totw] = (tk + 1ll) * tk % md * inv2 % md;
g1[totw]--;
if (k <= sqrtn) {
idx1[k] = totw;
} else idx2[r] = totw;
}
然后\(dp\)求出\(g(i, |P|)\)
for (int i = 1; i <= tot; i++) {
for (int j = 1; j <= totw; j++) {
if (pri[i] * pri[i] > w[j]) break;
ll k = w[j] / pri[i];
ll pre = k <= sqrtn ? idx1[k] : idx2[n / k];
g2[j] = (g2[j] - pri[i] * pri[i] % md * (g2[pre] - sum2[i-1] + md) % md + md) % md;
g1[j] = (g1[j] - pri[i] * (g1[pre]-sum1[i-1] + md) % md + md) % md;
}
}
接下来求\(S\)
ll S(ll i, ll j) {
//printf("i == %lld j == %lld\n", i, j);
//S(n,j)=\sum_{i=2}^{n}f(i)[i的最小质因子比第j个质数大]
if (pri[j] >= i)
return 0;
//质数部分的前缀和，pri[j]在sqrtn以内，一定是由idx1表示。
int id = i > sqrtn ? idx2[n/i] : idx1[i];
ll sp = ((g2[id] - g1[id] + sum1[j] - sum2[j]) % md + md) % md;
for (ll k = j + 1; k <= tot && pri[k] * pri[k] <= i; k++) { //超过sqrti的对答案没贡献
for (ll p = pri[k]; p <= i; p *= pri[k]) {
ll tmp = S(i/p, k) + (p > pri[k]);
ll tp = p % md;
sp = (sp + tp * (tp - 1 + md) % md * tmp) % md;
}
}
return sp;
}
最终目标是求出\(S(n, 0) + f(1)\)，由积性函数性质知\(f(1)=1\)，故答案为\((S(n, 0) + 1ll) \mod p\)
洛谷ac代码如下
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll maxn = 1e6 + 7, md = 1e9 + 7;
bool isp[maxn];
ll n, pri[maxn], idx1[maxn], idx2[maxn], sum1[maxn], sum2[maxn], w[maxn], sqrtn, g1[maxn], g2[maxn], inv2, inv3;
int tot, totw;
void sieve () {
memset(isp, true, sizeof(isp));
isp[1] = 0;
for (int i = 2; i <= 1000000; i++) {
if (isp[i]) {
pri[++tot] = i;
sum1[tot] = (sum1[tot-1] + pri[tot]) % md;
sum2[tot] = (sum2[tot-1] + pri[tot] * pri[tot] % md) % md;
}
for (int j = 1; j <= tot && i * pri[j] <= 1000000; j++) {
isp[i*pri[j]] = 0;
if (i % pri[j] == 0)
break;
}
}
}
ll ksm(ll a, int b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % md;
a = a * a % md;
b >>= 1;
}
return res % md;
}
ll S(ll i, ll j) {
//printf("i == %lld j == %lld\n", i, j);
//S(n,j)=\sum_{i=2}^{n}f(i)[i的最小质因子比第j个质数大]
if (pri[j] >= i)
return 0;
//质数部分的前缀和，pri[j]在sqrtn以内，一定是由idx1表示。
int id = i > sqrtn ? idx2[n/i] : idx1[i];
ll sp = ((g2[id] - g1[id] + sum1[j] - sum2[j]) % md + md) % md;
for (ll k = j + 1; k <= tot && pri[k] * pri[k] <= i; k++) { //超过sqrti的对答案没贡献
for (ll p = pri[k]; p <= i; p *= pri[k]) {
ll tmp = S(i/p, k) + (p > pri[k]);
ll tp = p % md;
sp = (sp + tp * (tp - 1 + md) % md * tmp) % md;
}
}
return sp;
}
int main() {
cin >> n;
sqrtn = sqrt(n);
inv2 = ksm(2, md - 2);
inv3 = ksm(3, md - 2);
sieve();
for (ll l = 1, r; l <= n; l = r + 1) {
ll k = n / l;
//printf("k == %lld\n", k);
r = n / k;
w[++totw] = k;
//求g(i,0)
//f'(i)=i^2-i
//g(i,0) = \sum_{j=2}^{n}f'(j) = 一次幂的前缀和减去二次幂的前缀和减去f'(1)
//这题f'(1)是0 不用减了。。。
ll tk = k % md;
g2[totw] = tk * (tk+1) % md * inv2 % md * (2*tk+1) % md * inv3 % md;
g2[totw]--;
g1[totw] = (tk + 1ll) * tk % md * inv2 % md;
g1[totw]--;
if (k <= sqrtn) {
idx1[k] = totw;
} else idx2[r] = totw;
}
//printf("tot == %d totw == %d\n", tot, totw);
//g数组滚动掉了第二维，因此要先枚举第二维
for (int i = 1; i <= tot; i++) {
for (int j = 1; j <= totw; j++) {
if (pri[i] * pri[i] > w[j]) break;
ll k = w[j] / pri[i];
ll pre = k <= sqrtn ? idx1[k] : idx2[n / k];
g2[j] = (g2[j] - pri[i] * pri[i] % md * (g2[pre] - sum2[i-1] + md) % md + md) % md;
g1[j] = (g1[j] - pri[i] * (g1[pre]-sum1[i-1] + md) % md + md) % md;
}
}
cout << (S(n, 0) + 1ll) % md;
return 0;
}
