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;
}