之前做题要用到min25筛，就断断续续地学了几下，用后即忘，简直就是浪费时间。不如现在好好记下来，巩固一下记忆。
在找博客学习过程中发现了一个写得非常好的博客：Min-25筛学习笔记 | LNRBHAW，配合Min_25 筛 - OI Wiki (oi-wiki.org)食用，效果很好。
作用和适用范围
min25可以以低于线性的复杂度计算出一类积性函数(f)的前缀和（不一定是局限于积性函数，见题目Hasse Diagram）。
要求对于素数(p)，(f(p))是一个关于(p)的低次多项式，并且(f(p^c))可以快速计算出来。例如(f(p)=p^2-p)，(f(p^k)=p^k(p^k-1))这样的函数。
注意，这个min25筛只能得到单点的值，没办法预处理出(sqrt{n})的值。
实现
定义：
-
({ m lpf}(n))为(n)的最小质因数。例如({ m lpf}(6)=2)，特别地，({ m lpf}(1)=1)。
-
(p_k)代表第(k)大的质数。
-
(S(n,k)=sumlimits_{i=2}^{n}{[p_k<{ m lpf}(i)]f(i)})
-
(F_{ m prime}(n)=sumlimits_{p_ile n}{f(p_i)})
这样，(S(n,1)+f(1))就是我们要求的前缀和。
(S(n,k))中可以分为两个部分分别计算，一个是素数的部分，一个合数部分。
-
素数：即为(S(n,k)=F_{ m prime}(n)-F_{ m prime}(p_{k-1}))
-
合数：通过枚举每个(i)的最小质因子及其次数就可以求得递推式（枚举最小质因子保证了不重不漏）。
素数部分
因为(f(p)=sum{a_kp^k})，有
令(g(p)=p^k)，可以发现(g(p))是完全积性函数（注意(g(p))系数是1，否则就不是完全积性函数了）。所以现在要求的是(F_{ m prime}(n)=sumlimits_{p_ile n}{g(p_i)})
设函数(G(n,i))为用前(i)个质数做埃式筛剩下的(g)值之和，这样(F_{ m prime}(n)=G(n,lfloorsqrt{n} floor))
它由有两个部分组成：
-
(p_i^2>n)，筛完了，(G(n,i)=G(n,i-1))
-
(p_i^2le n)，(p_i)的在区间([p_i,lfloor frac{n}{p_i} floor])且没被前(i-1)个素数筛掉的倍数会被筛掉，要减去。而这些倍数位于(G(lfloor frac{n}{p_i} floor,i-1))中，除此之外，(G(lfloor frac{n}{p_i} floor,i-1))还包含了(p_1,p_2,...,p_{i-1})这些值没被筛去，要从中减去(F_{ m prime}(p_{i-1}))
综上有
注意(g(p)=p^k,p{ m 是质数})，是完全积性函数，故解析式为(g(x)=x^k,xin Z^+)，可以O(1)求出其(g)前缀和(G(n,0))，(F_{ m prime}(p_{i-1}))也可以在O((sqrt{n}))预处理出来。
可以发现都是形如(lfloor frac{n}{p_i} floor)，可以使用数论分块+递推预处理出(2sqrt{n})项(F_{ m prime})值，分别为(F_{ m prime}(i),i=1,...,sqrt{n})和(F_{ m prime}(frac{n}{i}),i=1,...,sqrt{n})。这(2sqrt{n})个就是下面要用的值。
合数部分
枚举最小质数可得递推式如下：
由于(f(x))是积性函数，(f(p_i^c)S(frac{n}{p_i^c},i+1))相当于最小质因数为(p_i^c)的合数对应的值的和，最后还要补上一个(f(p_i^{c+1}))。
合并
最终结果即为：
直接递归计算即可。
细节
-
实现中将(F_{ m prime}(i),i=1,...,sqrt{n})这小于(sqrt{n})的映射到( m id1[i])，(F_{ m prime}(frac{n}{i}),i=1,...,sqrt{n})这些大于(n)的映射到( m id2[i])，相当于hashmap。
-
求素数部分用的是非递归写法，直接递推上去即可。
例题
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10;
const int M = 1e9 + 7;
const int INF = 0x3f3f3f3f;
typedef long long ll;
ll num[N], ansp1[N], ansp2[N], id1[N], id2[N];
int pri[N], cnt;
bool isnp[N];
void init() {
isnp[1] = 1;
pri[cnt++] = 1;
for(int i = 2; i < N; i++) {
if(!isnp[i]) {
pri[cnt++] = i;
for(int j = 2 * i; j < N; j += i)
isnp[j] = 1;
}
}
}
ll r2 = 500000004;
ll r6 = 166666668;
ll g1(ll x) { return x % M; }
ll sum1(ll x) { // 求和，用于求初始值S(n,0)
x %= M;
return x * (x + 1) % M * r2 % M;
}
ll g2(ll x) { return x % M * x % M; }
ll sum2(ll x) {
x %= M;
return x * (x + 1) % M * (2 * x + 1) % M * r6 % M;
}
ll f(ll x) { // 目标函数，x为p的幂次
x %= M;
return x * (x - 1) % M;
}
namespace min25 {
ll n, sqn;
int ID(ll x) { // hash映射
if(x <= sqn) return id1[x];
return id2[n / x];
}
void solve(ll (*g)(ll), ll (*sum)(ll), ll ansp[]) {
static ll sump[N]; // 预处理素数和
int cur = 0, mx = 0;
for(sqn = 0; sqn * sqn <= n; sqn++) {
if(pri[mx] <= sqn) {
mx++;
sump[mx] = (sump[mx - 1] + g(pri[mx])) % M;
}
}
ll i = 1;
while(i <= n) { // 数论分快求出映射关系和初始值
num[cur] = n / i;
ansp[cur] = (sum(num[cur]) - g(1)) % M;
if(num[cur] <= sqn) id1[num[cur]] = cur;
else id2[n / num[cur]] = cur;
cur++;
i = n / (n / i) + 1;
}
for(int i = 1; i < mx; i++) { // 直接递推
for(int j = 0; j < cur && 1ll * pri[i] * pri[i] <= num[j]; j++) {
ansp[j] -= g(pri[i]) * (ansp[ID(num[j] / pri[i])] - sump[i - 1]) % M;
ansp[j] = (ansp[j] % M + M) % M;
}
}
}
ll F(ll n, int k, ll (*f)(ll)) {
ll res = 0;
for(int i = k; 1ll * pri[i] * pri[i] <= n; i++) {
ll p = pri[i];
while(p * pri[i] <= n) {
res += (f(p) * F(n / p, i + 1, f) % M + f(p * pri[i])) % M;
res %= M;
p = p * pri[i];
}
}
// 合并素数部分
res += (ansp2[ID(n)] - ansp1[ID(n)]);
if(k - 1) res -= ansp2[ID(pri[k - 1])] - ansp1[ID(pri[k - 1])];
return (res % M + M) % M;
}
}
int main() {
init();
while(cin >> min25::n) {
min25::solve(g1, sum1, ansp1);
min25::solve(g2, sum2, ansp2);
cout << (min25::F(min25::n, 1, f) + 1) % M << endl;
}
}