亚线性筛
所谓的亚线性筛，是指一类用低于线性复杂度求出积性函数前缀和的方法。
杜教筛
杜教筛可以在\(O(n^{\frac{2}{3}})\)的时间复杂度求出\(\sqrt{n}\) 个点值，原理和实现都比较简单。
原理与实现
对于数论函数\(f\)，要求计算\(S(n)=\sum\limits_{i=1}^{n}{f(i)}\)。对于任意的数论函数\(g\)，有
如果可以快速求出\(\sum\limits_{i=1}^n{(g*f)(i)}\)以及\(g\)的\(\sqrt{n}\)个前缀和，再用数论分块递归求解即可得到\(S(n)\)。因此关键是找到合适的\(g\)，使得\(g*f\)可以快速求和。
假设求和计算复杂度为\(O(1)\)，直接计算时间复杂度为\(O(\int_0^{n^{\frac{1}{2}}}{\sqrt{\frac{n}{x}}=O(n^{\frac{3}{4}})})\)，如果使用线性筛预处理前\(n^{\frac{2}{3}}\)项的结果，时间复杂度就变为\(O(\int_0^{n^{\frac{1}{3}}}{\sqrt{\frac{n}{x}}=O(n^{\frac{2}{3}})})\)，可以处理\(10^{11}\)级别的数据。
常用的乘性函数的求和：
- 莫比乌斯函数：\(\mu * 1=\epsilon\)，\(\epsilon(n) = [n=1]\)
- 欧拉函数：\(\varphi * 1 = \operatorname{id}\)，\(\operatorname{id}(n)=n\)
- 约数个数：\(\operatorname{d} * \mu=1\)(需要\(\mu\)前缀和)
Powerful Number筛
Powerful Number(下面简称PN)筛可以说是杜教筛的扩展，可以结合杜教筛求一些积性函数的前缀和。
要求:
假设要求\(F(n)=\sum\limits_{i=1}^{n}{f(i)}\)，如果存在函数\(g\)满足：
- \(g\)是积性函数
- \(g\)前缀和\(G\)容易求得
- 对于质数\(p\)，\(g(p)=f(p)\)
那么就可以较快的求出\(F(n)\)。
原理与实现
\(f\)可以写成\(f=h*g\)，那么有(\(h\)和\(g\)都是积性函数)
如果可以使得\(p\)为质数时有\(h(p)=0\)，即\(f(p)=h(p)+g(p)=g(p)\)时，就有当\(d\)中含有次数为1的质因子时上面式子中\(h(d)=0\)，相当于把含有次数为1质因子的\(d\)全部都“筛掉”了。将所有不含1次质因子的数称为PN，就有
这是很强力的筛选，因为\(n\)以内的PN只有\(O(\sqrt{n})\)个。
PN筛分为以下几个步骤：
- 构造\(g\)
注意构造的\(g\)要是积性函数，且满足\(f(p)=g(p)\)，比如
-
\(f(p^c)=p\)，那么构造\(g(x)=x\)。
-
\(f(p^c)=p^c(p^c-1)\)，那么构造\(g(x)=id(x)\varphi(x)\)。
\(g\)的前缀和\(G\)要好求，注意到需要的是\(\lfloor\frac{n}{i}\rfloor\)处的\(G\)值，所以常用杜教筛。
- 计算\(h(p^c)\)
知道了\(g\)和\(f\)就可以求出\(h\)，由
可得
可以枚举\(p\)和\(c\)递推预处理出来。
- 搜索PN，过程中结果累和起来
直接枚举质因子以及次数搜索，由于PN只有\(O(\sqrt{n})\)个，只用搜索\(O(\sqrt{n})\)次。全部结果加起来就是答案。
如果处理\(G\)的时间复杂度为\(O(n^{\alpha})\)，那么总时间复杂度为\(O(\max(\sqrt{n},n^{\alpha}))\)。
例题
\(f(p^c)=p^c(p^c-1)\)，那么构造\(g(x)=id(x)\varphi(x)\)。
可以发现\((id\cdot\varphi)*id=id^2\)，所以可以使用杜教筛，即
然后各种预处理\(h(p^c)\)，直接搜索PN累加答案即可。
#include <bits/stdc++.h>
#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N)
typedef long long ll;
using namespace std;
/*-----------------------------------------------------------------*/
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f
const int N = 6e6 + 10;
const ll MX = 1e10 + 10;
const int M = 1e9 + 7;
const double eps = 1e-5;
ll pri[N];
bool isnp[N];
int cnt;
inline ll qpow(ll a, ll b, ll m) {
ll res = 1;
while(b) {
if(b & 1) res = (res * a) % m;
a = (a * a) % m;
b = b >> 1;
}
return res;
}
namespace PNS {
map<ll, ll> mpG;
ll g[N], G[N];
ll h[N][35];
ll global_n, ans;
ll r6, r2;
void init() { // 预处理素数和h[p^c]
r2 = qpow(2, M - 2, M);
r6 = qpow(6, M - 2, M);
isnp[1] = 1;
g[1] = 1;
for(int i = 2; i < N; i++) {
if(!isnp[i]) {
pri[cnt++] = i;
g[i] = 1ll * i * (i - 1) % M;
}
for(int j = 0; j < cnt && i * pri[j] < N; j++) {
isnp[i * pri[j]] = 1;
if(i % pri[j] == 0) {
g[i * pri[j]] = g[i] * pri[j] % M * pri[j] % M;
break;
}
g[i * pri[j]] = g[i] * g[pri[j]] % M;
}
}
for(int i = 1; i < N; i++) G[i] = (g[i] + G[i - 1]) % M;
for(int i = 0; i < cnt; i++) {
h[i][0] = 1;
h[i][1] = 0;
}
for(int i = 0; i < cnt; i++) {
ll pp = pri[i] * pri[i], rp = qpow(pri[i], M - 2, M);
for(int c = 2; c < 35; c++) {
ll pw = pp % M;
h[i][c] = pw * (pw - 1) % M;
for(int k = 0; k < c; k++) {
h[i][c] = (h[i][c] - h[i][k] * pw % M * pw % M * (pri[i] - 1) % M * rp % M + M) % M;
pw = pw * rp % M;
}
if(pp > MX / pri[i]) break;
pp = pp * pri[i];
}
}
}
ll sum(ll n) { // 杜教筛预处理
if(n < N) return G[n];
if(mpG.count(n)) return mpG[n];
ll res = n % M * (n % M + 1) % M * (2 * n % M + 1) % M * r6 % M;
ll i = 2;
while(i <= n) {
ll j = n / (n / i);
res = (res - (j - i + 1) % M * ((i + j) % M) % M * r2 % M * sum(n / i) % M + M) % M;
i = j + 1;
}
mpG[n] = res;
return res;
}
void dfs(int p, ll hv, ll v) { // 搜索PN， hv为h函数值，v为当前找到的PN
ans = (ans + hv * sum(global_n / v) % M) % M;
for(int i = p; i < cnt && v <= global_n / pri[i] / pri[i]; i++) {
ll tv = v * pri[i] * pri[i];
for(int c = 2; ; c++) {
dfs(i + 1, hv * h[i][c] % M, tv);
if(tv > global_n / pri[i]) break;
tv = tv * pri[i];
}
}
}
ll solve(ll n) {
global_n = n;
sum(n);
ans = 0;
dfs(0, 1, 1);
return ans;
}
}
int main() {
IOS;
PNS::init();
ll n;
cin >> n;
cout << PNS::solve(n) << endl;
}
