所谓的亚线性筛，是指一类用低于线性复杂度求出积性函数前缀和的方法。
杜教筛
杜教筛可以在(O(n^{frac{2}{3}}))的时间复杂度求出(sqrt{n}) 个点值，原理和实现都比较简单。
原理与实现
对于数论函数(f)，要求计算(S(n)=sumlimits_{i=1}^{n}{f(i)})。对于任意的数论函数(g)，有
如果可以快速求出(sumlimits_{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)=sumlimits_{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)要好求，注意到需要的是(lfloorfrac{n}{i} floor)处的(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))。
可以发现((idcdotvarphi)*id=id^2)，所以可以使用杜教筛，即
然后各种预处理(h(p^c))，直接搜索PN累加答案即可。
#include <bits/stdc++.h>
#define endl '
'
#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;
}