洲阁筛
给定一个积性函数$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一样，分成三段，一段暴力转移，一段需要时前缀和，一段不需要维护。
复杂度也是同样的计算方法
1 /** 2 * SPOJ 3 * Problem#DIVCNT3 4 * Accepted 5 * Time: 35190ms 6 * Memory: 33792k 7 */ 8 #include <bits/stdc++.h> 9 using namespace std; 10 typedef bool boolean; 11 12 template <typename T> 13 void pfill(T* pst, const T* ped, T val) { 14 for ( ; pst != ped; *(pst++) = val); 15 } 16 17 #define ll long long 18 19 typedef class Euler { 20 public: 21 int num; 22 boolean *vis; 23 int *pri, *d3, *mp; 24 25 Euler(int n) : num(0) { 26 d3 = new int[(n + 1)]; 27 mp = new int[(n + 1)]; 28 pri = new int[(n + 1)]; 29 vis = new boolean[(n + 1)]; 30 pfill(vis, vis + n + 1, false); 31 d3[1] = 1; 32 for (int i = 2; i <= n; i++) { 33 if (!vis[i]) { 34 d3[i] = 4, pri[num++] = i, mp[i] = 1; 35 } 36 for (int j = 0, _ = n / i, x; j < num && pri[j] <= _; j++) { 37 vis[x = pri[j] * i] = true; 38 if (!(i % pri[j])) { 39 d3[x] = (d3[i / mp[i]] + 3) * d3[mp[i]]; 40 mp[x] = mp[i]; 41 break; 42 } else { 43 mp[x] = i, d3[x] = d3[i] << 2; 44 } 45 } 46 } 47 } 48 ~Euler() { 49 delete[] mp; 50 delete[] vis; 51 delete[] d3; 52 delete[] pri; 53 } 54 } Euler; 55 56 const int N = 330000; 57 58 Euler *_euler; 59 60 int T; 61 vector<ll> ns; 62 63 ll maxn, D; 64 65 int cnt_prime; 66 int *pri, *d3; 67 68 inline void init() { 69 scanf("%d", &T); 70 for (ll x; T--; ) { 71 scanf("%lld", &x); 72 ns.push_back(x); 73 maxn = max(maxn, x); 74 } 75 D = sqrt(maxn + 0.5) + 1; 76 while (D * 1ll * D <= maxn) 77 D++; 78 _euler = new Euler(D + 23); 79 cnt_prime = _euler->num; 80 pri = _euler->pri, d3 = _euler->d3; 81 } 82 83 ll n; 84 int cnt_P[N]; 85 int range0[N], range1[N]; 86 ll f0[N], f1[N], g0[N], g1[N]; 87 88 void precalc() { 89 for (int i = 1; i < D; i++) { 90 range0[i] = range0[i - 1]; 91 while (pri[range0[i]] * 1ll * pri[range0[i]] <= i) { 92 range0[i]++; 93 } 94 } 95 for (int i = D - 1; i; i--) { 96 ll y = n / i; 97 range1[i] = range1[i + 1]; 98 while (pri[range1[i]] * 1ll * pri[range1[i]] <= y) { 99 range1[i]++; 100 } 101 } 102 pfill(cnt_P, cnt_P + D, 0); 103 for (int i = 0; pri[i] < D; i++) { 104 cnt_P[pri[i]]++; 105 } 106 for (int i = 1; i < D; i++) { 107 cnt_P[i] += cnt_P[i - 1]; 108 } 109 } 110 111 int nump; 112 ll get_value_f(int i, ll j) { 113 if (j >= D) { 114 int rj = n / j; 115 return f1[rj] + (max(0, nump - max(range1[rj], i)) << 2); 116 } 117 return f0[j] + (max(0, cnt_P[j] - max(range0[j], i)) << 2); 118 } 119 120 void calculate_f() { 121 for (nump = 0; pri[nump] < D; nump++); 122 for (int i = 1; i < D; i++) { 123 f0[i] = 1, f1[i] = 1; 124 } 125 for (int i = nump - 1; ~i; i--) { 126 for (int j = 1; j < D && i < range1[j]; j++) { 127 ll m = n / j, coef = 1; 128 while (m /= pri[i]) { 129 f1[j] += (coef += 3) * get_value_f(i + 1, m); 130 } 131 } 132 for (int j = D - 1; j && i < range0[j]; j--) { 133 ll m = j, coef = 1; 134 while (m /= pri[i]) { 135 f0[j] += (coef += 3) * get_value_f(i + 1, m); 136 } 137 } 138 } 139 for (int i = 1; i < D; i++) { 140 f1[i] = get_value_f(0, n / i); 141 } 142 for (int i = 1; i < D; i++) { 143 f0[i] = get_value_f(0, i); 144 } 145 } 146 147 ll get_value_g(int i, ll j) { 148 if (j >= D) { 149 int rj = n / j; 150 return g1[rj] - max(0, i - range1[rj]); 151 } 152 return g0[j] - max(0, min(i, cnt_P[j]) - range0[j]); 153 } 154 155 void calculate_g() { 156 for (int i = 1; i < D; i++) { 157 g0[i] = i, g1[i] = n / i; 158 } 159 int cp = 0; 160 for (int i = 0; pri[i] < D; i++, cp++) { 161 for (int j = 1; j < D && i < range1[j]; j++) { 162 g1[j] -= get_value_g(i, n / j / pri[i]); 163 } 164 for (int j = D - 1; j && i < range0[j]; j--) { 165 g0[j] -= get_value_g(i, j / pri[i]); 166 } 167 } 168 for (int i = 1; i < D; i++) { 169 (g1[i] = get_value_g(cp, n / i) - 1) <<= 2; 170 } 171 for (int i = 1; i < D; i++) { 172 (g0[i] = get_value_g(cp, i) - 1) <<= 2; 173 } 174 } 175 176 void Solve(ll _n) { 177 n = _n; 178 D = sqrt(n + 0.5) + 1; 179 while (D * 1ll * D <= n) 180 D++; 181 precalc(); 182 calculate_g(); 183 calculate_f(); 184 ll ans = f1[1]; 185 for (int i = 1; i < D; i++) { 186 ans += d3[i] * g1[i]; 187 // cerr << d3[i] << " " << g1[i] << '\n'; 188 } 189 printf("%lld\n", ans); 190 } 191 192 inline void solve() { 193 for (int i = 0; i < (signed) ns.size(); i++) { 194 Solve(ns[i]); 195 } 196 } 197 198 int main() { 199 init(); 200 solve(); 201 return 0; 202 }