min_25筛学习笔记
作用
对积性函数 \(F(n)\) 求前缀和。
需要满足的条件：
- 对于 \(p\) 是素数，\(F(p)\) 是多项式。
- 对于 \(p\) 是素数，\(F(p^k)\) 可以快速计算（由k从小到大递推计算也可以）。
复杂度据说是 \(O(\frac{n^{3/4}}{logn})\) ，在 \(n \le 10^{12}\) 基本上都能快速出解。空间复杂度是 \(O(\sqrt n)\)
原理
分为两个阶段：预处理、计算答案。思路是先求所有素数的函数和，再求合数的，最后考虑 \(F(1)\) 。
预处理
由于 \(F(p)\) 是多项式，所以可以拆分成若干个 \(p^k\) 的形式。记 \(f(p)=p^k\) （这里的 p 可以不是素数），下面考虑求 \(\sum_{i=1}^n f(i)[\text{i is prime}]\)
\(g(n,j)\) 表示小于等于 n 中所有质数的 \(f\) 和最小素因子大于 \(P_j\) 的合数的 \(f\) 之和。易得
对于下面那个式子，减号后面利用了 \(f\) 是完全积性函数的性质。即把 \(P_j^k\) 提出来，剩下的数字的最小质因数也要大于 \(P_j\) ，于是中括号内是最小质因子大于 \(P_j\) 的所有数的贡献。实际上也可以写成
由于 \(\lfloor\frac{n}{i}\rfloor\) 只会取根号级别个，所以只需要对这些数计算 \(g\) 即可。而第二维，由于 \(P_j^2>n\) 的所有 j 没有贡献，所以只会取到 \(O(\sqrt n / \ln n)\) 。于是以 j 为阶段就可以进行 dp。
最终要求的 \(\sum_{i=1}^n f(i)[\text{i is prime}]\) 就是 \(g(n,|P|)\) 。
初始化：\(g(n,0) = \sum_{i=2}^nf(i)\)
计算答案
与预处理中类似符号的定义不同，这里定义 \(S(n,j)\) 是小于等于 n 中最小素因子大于等于 \(P_j\) 的数的 \(F\) 之和。则
其中 \(G(n,|P|)\) 为 \(\sum_{i=1}^n F(i)[\text{i is prime}]\) 。这其实是合并了预处理中若干个 \(f(p)=p^k\) 的结果。
这个式子即：大于等于 \(P_j\) 的素数的 \(F\) 之和，加上最小质因子大于等于 \(P_j\) 的合数的 \(F\) 之和。
于是最终答案为 \(S(n,1) + F(1)\) 。
应用
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef long double LD;
typedef pair<int,int> pii;
typedef pair<LL,int> pli;
const int SZ = 1e7 + 10;
const int INF = 1e9 + 10;
const int mod = 1e9 + 7;
const LD eps = 1e-8;
LL read() {
LL n = 0;
char a = getchar();
bool flag = 0;
while(a > '9' || a < '0') { if(a == '-') flag = 1; a = getchar(); }
while(a <= '9' && a >= '0') { n = n * 10 + a - '0',a = getchar(); }
if(flag) n = -n;
return n;
}
// 此题：f(p^k)=p^k(p^k-1)
struct Min_25 {
LL n;
int gen; // sqrt(n)
LL w[SZ]; // 所有的 n/i 下取整
int id1[SZ],id2[SZ]; // w中对应数字的标号
int sump[SZ],sump2[SZ]; // 前 i 个素数的和、平方和
int g1[SZ],g2[SZ]; // 计算完毕后为小于等于 w[i] 的所有质数的对应函数和
bool vis[SZ];
int pri[SZ],tot;
void shai(int n) {
tot = 0;
for(int i = 2;i <= n;i ++) {
if(!vis[i]) pri[++ tot] = i;
for(int j = 1,m;j <= tot && (m = i * pri[j]) <= n;j ++) {
vis[m] = 1;
if(i%pri[j] == 0) break;
}
}
for(int i = 1;i <= tot;i ++) {
sump[i] = (sump[i-1] + pri[i]) % mod; /// modify
sump2[i] = (sump2[i-1] + 1ll*pri[i]*pri[i]%mod) % mod; /// modify
}
}
int getid(LL x) {
if(x <= gen) return id1[x];
return id2[n/x];
}
int S(LL x,int y) {
if(x <= 1 || pri[y] > x) return 0;
int k = getid(x);
int ans = (g2[k]-g1[k]-(sump2[y-1]-sump[y-1])) % mod; /// modify
for(int i = y;i <= tot && 1ll * pri[i] * pri[i] <= x;i ++) {
LL p = pri[i];
LL t1 = p,t2 = p*p;
for(int e = 1;t2 <= x;e ++) {
(ans += (S(x/t1,i+1) * (t1 % mod)%mod * ((t1-1) % mod) % mod /// modify
+ (t2%mod) * ((t2-1)%mod) % mod)%mod) %=mod;
t1 *= p; t2 *= p;
}
}
return ans;
}
int f1(LL x) {
x %= mod;
return x * (x+1) / 2 % mod;
}
int f2(LL x) {
x %= mod;
int ni6 = 166666668;
return x * (x+1) % mod * (2*x+1) % mod * ni6 % mod;
}
int work(LL nn) {
n = nn;
gen = sqrt(n);
shai(gen);
int m = 0;
for(LL i = 1,r;i <= n;i = r + 1) {
r = n/(n/i); w[++m] = n/i;
g1[m] = f1(w[m])-1; /// modify
g2[m] = f2(w[m])-1; /// modify
if(w[m] <= gen) id1[w[m]]=m;
else id2[r]=m;
}
for(int j = 1;j <= tot;j ++) {
LL p = pri[j];
for(int i = 1;i<=m && p*p<=w[i];i ++) {
int k = getid(w[i]/p);
(g1[i] -= p * (g1[k] - sump[j-1]) % mod) %= mod; /// modify
(g2[i] -= p * p % mod * (g2[k] - sump2[j-1]) % mod) %= mod; /// modify
}
}
int ans = S(n,1) + 1;
ans += mod; ans %= mod;
return ans;
}
}min25;
int main() {
LL n = read();
printf("%d\n",min25.work(n));
}
