<!--
.mycontent {
border-left-style: solid;
border-left-width: 10px;
border-left-color: rgb(70, 131, 255);
padding-left: 10px;
}
h3.subtitle {
padding-top: 3px;
padding-bottom: 3px;
background-color: rgba(70, 131, 255, 0.5);
}
.warning {
padding-left: 30px;
padding-top: 3px;
padding-bottom: 3px;
background-color: rgba(255, 255, 83, 0.5);
}
.theory {
padding-left: 30px;
padding-top: 3px;
padding-bottom: 3px;
background-color: rgba(0, 255, 64, 0.5);
}
.methon {
padding-left: 30px;
padding-top: 3px;
padding-bottom: 3px;
background-color: rgba(135, 206, 235, 0.5);
}
-->
洲阁筛
给定一个积性函数$F(n)$，求$\sum_{i = 1}^{n}F(n)$。并且$F(n)$满足在素数和素数次幂的时候易于计算。
显然有：
$\sum_{i = 1}^{n} F(n) = \sum_{i = 1}^{\sqrt{n}}F(i) \left(\sum_{\sqrt{n} < p\leqslant n/i, p\ is\ a\ prime} F(p) \right) + \sum_{i = 1, i\ has\ no\ prime\ factor\ greater\ than\ \sqrt{n}}^{n} F(i)$。
Part 1
第一部分是要求出$\sum_{\sqrt{n}<p \leqslant n/i, p\ is\ a\ prime} F(p)$
假设小于等于$\sqrt{n}$的素数从小到大分别是$P_1, P_2, \cdots, P_{m}$。
设$f_{i, j}$表示$[1, j]$中和前$i$个素数互质的数的函数值和。
当$i = 0$时候想办法计算。
考虑当$i > 0$的时候的转移。考虑从$f_{i - 1, j}$中需要删掉最小质因子等于$P_i$的数的贡献。
所以有转移式子$f_{i, j} = f_{i - 1, j} - F(P_i) f_{i - 1, \left \lfloor j / P_i \right \rfloor}$
如果$F(n)$本身不是完全积性函数，需要把它拆成若干积性函数的和。注意这一部分只用保证素数处取值相同即可。
注意到第一维约有$O\left (\frac{\sqrt{n}}{\log n} \right)$种取值，第二维有$O\left (\sqrt{n}\right)$种取值。暴力计算的复杂度是$O\left (\frac{n}{\log n} \right )$
优化
考虑当$j < P_{i+ 1}$的时候，$f_{i, j} = 1$。考虑把条件放缩成$j < P_i$。
因此，当$P_i \leqslant j < P_{i}^2$的时候$f_{i, j} = f_{i - 1, j} - F(P_i)$。
所以做法是：
维护一个$F(P_i)$的前缀和
考虑转移时的$f_{i - 1, \left \lfloor j / P_i \right \rfloor}$
当$j \geqslant P_{i}^2$的时候暴力计算
当$P_i \leqslant j < P_{i}^2$的时候，需要这样的状态的时候减去一段和就行了。
当$j < P_i$的时候，由于$f_{i, 0} = 0$，所以这一部分的值不会改变。
所以对于每个$j$，有效的状态只有$O\left (\frac{\sqrt{j}}{\log j} \right )$。
时间复杂度为：$\sum_{i = 1}^{\sqrt{n}} O \left (\frac{\sqrt{i}}{\log i} \right ) + \sum_{i = 1}^{\sqrt{n}} O\left ( \frac{\sqrt{n/i}}{\log {n/i}}\right)$
由于第一部分显然小于第二部分，所以只用考虑第二部分的和。
所以有：$\sum_{i = 1}^{\sqrt{n}}O \left (\frac{\sqrt{n/i}}{\log n - \log i}\right)$
由于$\log i \leqslant \log \sqrt{n} = \frac{1}{2}\log n$
所以：$\sum_{i = 1}^{\sqrt{n}}O \left (\frac{\sqrt{n/i}}{\frac{1}{2}\log n}\right) = \sum_{i = 1}^{\sqrt{n}}O\left (\frac{\sqrt{n / i}}{\log n}\right) = O\left ( \frac{n^{3/4}}{\log n}\right)$
Part 2
这一部分是要求出$\sum_{i = 1, i\ has\ no\ prime\ factor\ greater\ than\ \sqrt{n}}^{n} F(i)$
假设小于等于$\sqrt{n}$的素数从小到大分别是$P_1, P_2, \cdots, P_{m}$。
设$g_{i, j}$表示在$[1, j]$中，质因子只包含$P_{i}, P_{i + 1}, \cdots, P_{m}$的数的函数值的和。
转移考虑枚举当前质数$P_i$的指数。所以有：$g_{i, j} = g_{i + 1, j} + \sum_{e = 1}^{\log_{P_i} j} F(P_{i}^{e}) g_{i + 1, \lfloor j / P_i^e\rfloor}$
对于时间复杂度，考虑对于每一个$j$的转移的枚举量，可以得到大概是$O(\frac{\sqrt{j}}{\log j})$。(不会积分，告辞)
和Part 1一样，得到暴力计算这一部分的复杂度是$O(\frac{n}{\log n})$。
优化
因为质因子是从大到小枚举的，所以当$j < P_{i}$的时候，$g_{i, j} = 1$。
所以当满足$P_i \leqslant j < P_{i}^2$，满足$g_{i, j} = g_{i + 1, j} + F(P_i)$。
和Part 1一样，分成三段，一段暴力转移，一段需要时前缀和，一段不需要维护。
复杂度也是同样的计算方法
/**
* SPOJ
* Problem#DIVCNT3
* Accepted
* Time: 35190ms
* Memory: 33792k
*/
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean; template <typename T>
void pfill(T* pst, const T* ped, T val) {
for ( ; pst != ped; *(pst++) = val);
} #define ll long long typedef class Euler {
public:
int num;
boolean *vis;
int *pri, *d3, *mp; Euler(int n) : num() {
d3 = new int[(n + )];
mp = new int[(n + )];
pri = new int[(n + )];
vis = new boolean[(n + )];
pfill(vis, vis + n + , false);
d3[] = ;
for (int i = ; i <= n; i++) {
if (!vis[i]) {
d3[i] = , pri[num++] = i, mp[i] = ;
}
for (int j = , _ = n / i, x; j < num && pri[j] <= _; j++) {
vis[x = pri[j] * i] = true;
if (!(i % pri[j])) {
d3[x] = (d3[i / mp[i]] + ) * d3[mp[i]];
mp[x] = mp[i];
break;
} else {
mp[x] = i, d3[x] = d3[i] << ;
}
}
}
}
~Euler() {
delete[] mp;
delete[] vis;
delete[] d3;
delete[] pri;
}
} Euler; const int N = ; Euler *_euler; int T;
vector<ll> ns; ll maxn, D; int cnt_prime;
int *pri, *d3; inline void init() {
scanf("%d", &T);
for (ll x; T--; ) {
scanf("%lld", &x);
ns.push_back(x);
maxn = max(maxn, x);
}
D = sqrt(maxn + 0.5) + ;
while (D * 1ll * D <= maxn)
D++;
_euler = new Euler(D + );
cnt_prime = _euler->num;
pri = _euler->pri, d3 = _euler->d3;
} ll n;
int cnt_P[N];
int range0[N], range1[N];
ll f0[N], f1[N], g0[N], g1[N]; void precalc() {
for (int i = ; i < D; i++) {
range0[i] = range0[i - ];
while (pri[range0[i]] * 1ll * pri[range0[i]] <= i) {
range0[i]++;
}
}
for (int i = D - ; i; i--) {
ll y = n / i;
range1[i] = range1[i + ];
while (pri[range1[i]] * 1ll * pri[range1[i]] <= y) {
range1[i]++;
}
}
pfill(cnt_P, cnt_P + D, );
for (int i = ; pri[i] < D; i++) {
cnt_P[pri[i]]++;
}
for (int i = ; i < D; i++) {
cnt_P[i] += cnt_P[i - ];
}
} int nump;
ll get_value_f(int i, ll j) {
if (j >= D) {
int rj = n / j;
return f1[rj] + (max(, nump - max(range1[rj], i)) << );
}
return f0[j] + (max(, cnt_P[j] - max(range0[j], i)) << );
} void calculate_f() {
for (nump = ; pri[nump] < D; nump++);
for (int i = ; i < D; i++) {
f0[i] = , f1[i] = ;
}
for (int i = nump - ; ~i; i--) {
for (int j = ; j < D && i < range1[j]; j++) {
ll m = n / j, coef = ;
while (m /= pri[i]) {
f1[j] += (coef += ) * get_value_f(i + , m);
}
}
for (int j = D - ; j && i < range0[j]; j--) {
ll m = j, coef = ;
while (m /= pri[i]) {
f0[j] += (coef += ) * get_value_f(i + , m);
}
}
}
for (int i = ; i < D; i++) {
f1[i] = get_value_f(, n / i);
}
for (int i = ; i < D; i++) {
f0[i] = get_value_f(, i);
}
} ll get_value_g(int i, ll j) {
if (j >= D) {
int rj = n / j;
return g1[rj] - max(, i - range1[rj]);
}
return g0[j] - max(, min(i, cnt_P[j]) - range0[j]);
} void calculate_g() {
for (int i = ; i < D; i++) {
g0[i] = i, g1[i] = n / i;
}
int cp = ;
for (int i = ; pri[i] < D; i++, cp++) {
for (int j = ; j < D && i < range1[j]; j++) {
g1[j] -= get_value_g(i, n / j / pri[i]);
}
for (int j = D - ; j && i < range0[j]; j--) {
g0[j] -= get_value_g(i, j / pri[i]);
}
}
for (int i = ; i < D; i++) {
(g1[i] = get_value_g(cp, n / i) - ) <<= ;
}
for (int i = ; i < D; i++) {
(g0[i] = get_value_g(cp, i) - ) <<= ;
}
} void Solve(ll _n) {
n = _n;
D = sqrt(n + 0.5) + ;
while (D * 1ll * D <= n)
D++;
precalc();
calculate_g();
calculate_f();
ll ans = f1[];
for (int i = ; i < D; i++) {
ans += d3[i] * g1[i];
// cerr << d3[i] << " " << g1[i] << '\n';
}
printf("%lld\n", ans);
} inline void solve() {
for (int i = ; i < (signed) ns.size(); i++) {
Solve(ns[i]);
}
} int main() {
init();
solve();
return ;
}
洲阁筛
Min25筛
Part 1
假设小于等于$\sqrt{n}$的素数从小到大分别是$P_1, P_2, \cdots, P_{m}$。
设$g_{i, j} = \sum_{k = 2}^{j} [k的最小质因子大于P_i或者k是质数] F(k)$
转移式子：
$$g_{i, j} =
\begin{cases}g_{i - 1, j} - F(P_i)(g_{i - 1, \lfloor j / P_i\rfloor} - g_{i - 1, P_i - 1}) & (j \geqslant P_i^2) \\
g_{i - 1, j} & (1 \leqslant j < P_i^2)
\end{cases}
$$
和洲阁筛类似，不过第一种情况还需要抠掉素数的贡献。
这里也要求$F$是一个完全积性函数。
我们考虑求出$1, 2, \cdots, \sqrt{n}, \lfloor n / \sqrt{n}\rfloor, \cdots, \lfloor n / 2\rfloor, n$处的值。
Part 2
这里考虑计算$S(n, i) = \sum_{i = 2}^{n}[n的最小质因子大于等于P_i]F(n)$。
那么显然答案等于$S(n, 1) + F(1)$。
这里简单粗暴一点，如果是素数，那么可以用第一部分计算的答案再减去某个前缀和得到，否则枚举这个合数的最小质因子。不难得到这样一个式子：
$$
S(n, i) = g_{m, n} - \sum_{j = 1}^{i - 1}F(P_j) + \sum_{j = i \wedge P_j^2 \leqslant n}\sum_{e = 1\wedge P_{j}^{e + 1}\leqslant n} F(P_{j}^{e}) S(\lfloor n/P_{j}^e\rfloor, j + 1) + F(P_{j}^{e + 1})
$$
前一半是计算素数的贡献。
考虑计算合数的贡献，后一半是先枚举最小质因子，因为它是一个合数最小质因子，所以$p^2 \leqslant n$
$p^{e + 1} \leqslant n$是类似的原因。
时间复杂度被证明为是$O(\frac{n}{Poly(\log n)})$。
（似乎目前min_25在$10^{13}$范围内吊打其他筛）
/**
* loj
* Problem#6053
* Accepted
* Time: 332ms
* Memory: 2684k
*/
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean; #define ll long long void exgcd(int a, int b, int& x, int& y) {
if (!b) {
x = 1, y = 0;
} else {
exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
} int inv(int a, int n) {
int x, y;
exgcd(a, n, x, y);
return (x < 0) ? (x + n) : (x);
} const int Mod = 1e9 + 7; template <const int Mod = :: Mod>
class Z {
public:
int v; Z() : v(0) { }
Z(int x) : v(x){ }
Z(ll x) : v(x % Mod) { } friend Z operator + (const Z& a, const Z& b) {
int x;
return Z(((x = a.v + b.v) >= Mod) ? (x - Mod) : (x));
}
friend Z operator - (const Z& a, const Z& b) {
int x;
return Z(((x = a.v - b.v) < 0) ? (x + Mod) : (x));
}
friend Z operator * (const Z& a, const Z& b) {
return Z(a.v * 1ll * b.v);
}
friend Z operator ~(const Z& a) {
return inv(a.v, Mod);
}
friend Z operator - (const Z& a) {
return Z(0) - a;
}
Z& operator += (Z b) {
return *this = *this + b;
}
Z& operator -= (Z b) {
return *this = *this - b;
}
Z& operator *= (Z b) {
return *this = *this * b;
}
friend boolean operator == (const Z& a, const Z& b) {
return a.v == b.v;
}
}; Z<> qpow(Z<> a, int p) {
Z<> rt = Z<>(1), pa = a;
for ( ; p; p >>= 1, pa = pa * pa) {
if (p & 1) {
rt = rt * pa;
}
}
return rt;
} typedef Z<> Zi; const Zi inv2 ((Mod + 1) >>1); const int Dmx = 1e5 + 5; ll n;
int D;
int pnum;
int pri[Dmx], sf[Dmx]; void Euler(int n) {
static bitset<Dmx + 23> vis;
for (int i = 2; i <= n; i++) {
if (!vis.test(i)) {
pri[pnum++] = i;
}
for (int *p = pri, *_ = pri + pnum, x; p != _ && (x = *p * i) <= n; p++) {
vis.set(x);
if (!(i % *p)) {
break;
}
}
}
for (int i = 1; i <= pnum; i++) {
sf[i] = sf[i - 1] + (pri[i - 1] ^ 1);
}
} template <typename T>
class SieveArray {
public:
ll n;
vector<T> a0;
vector<T> a1; void init(ll n, int D) {
this->n = n;
a0.assign(D, 0);
a1.assign(D, 0);
} T& operator [] (ll p) {
return (p >= (signed) a0.size()) ? a1[n / p] : a0[p];
}
}; SieveArray<int> range;
SieveArray<Zi> f0, f1; Zi S(ll n, int m) {
if (pri[m] > n) {
return 0;
}
Zi rt = f1[n] - sf[m];
for (int j = m; j < pnum && 1ll * pri[j] * pri[j] <= n; j++) {
int p = pri[j];
ll nn = n, p2 = 1ll * p * p;
for (int e = 1; p2 <= nn; e++) {
rt += S(nn /= p , j + 1) * (p ^ e) + (p ^ (e + 1));
}
}
return rt;
} int main() {
scanf("%lld", &n);
D = sqrt(D + 0.5) + 1;
while (1ll * D * D <= n)
D++;
Euler(D + 3);
range.init(n, D);
f0.init(n, D);
f1.init(n, D);
ll t = 0;
range[0] = 0;
for (int i = 1; i < D; i++) {
while (t * t <= i)
t++;
range.a0[i] = t;
}
t = 0;
for (int i = D; --i; ) {
ll v = n / i;
while (t * t <= v)
t++;
range.a1[i] = t;
}
for (int i = 1; i < D; i++) {
f0.a0[i] = i - 1;
f1.a0[i] = ((Zi(i) * (i + 1)) * inv2) - 1;
ll v = n / i;
f0.a1[i] = v - 1;
f1.a1[i] = ((Zi(v) * (v + 1)) * inv2) - 1;
}
for (int i = 0; i < pnum; i++) {
int p = pri[i];
if (p >= D) {
break;
}
ll v;
for (int j = 1; j < D && p < range.a1[j]; j++) {
v = n / (1ll * j * p);
f0.a1[j] -= f0[v] - f0.a0[p - 1];
f1.a1[j] -= (f1[v] - f1.a0[p - 1]) * p;
}
for (int j = D; p < range.a0[--j]; ) {
f0.a0[j] -= f0.a0[j / p] - f0.a0[p - 1];
f1.a0[j] -= (f1.a0[j / p] - f1.a0[p - 1]) * p;
}
}
for (int i = 1; i < D; i++) {
f1.a0[i] -= f0.a0[i];
f1.a0[i] += (i >= 2) * 2;
f1.a1[i] -= f0.a1[i];
f1.a1[i] += (2 * i <= n) * 2;
}
Zi ans = S(n, 0) + 1;
printf("%d\n", ans.v);
return 0;
}
