Min_25筛有关求解次小质因子

#188. 【UR #13】Sanrd

题意化简就是求次小质因子，这一步我们可以在Min_25筛的ans计算中得到，

S(n, j)表示的是最小质因子大于等于primejprime_jprimej的加上质数的答案贡献，

要满足次小质因子，一定有除去这个数之后只剩下质数了，显然当我们计算到S(n, j)时，我们可以计算primej−1prime_{j - 1}primej−1对答案的贡献，

接着我们还需要加上质数的次方项对答案的贡献，最后我们就可以得到答案了。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1using namespace std;typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;inline ll read() {ll f = 1, x = 0;char c = getchar();while(c < '0' || c > '9') {if(c == '-')    f = -1;c = getchar();}while(c >= '0' && c <= '9') {x = (x << 1) + (x << 3) + (c ^ 48);c = getchar();}return f * x;
}const int N = 1e6 + 10;namespace Min_25 {int prime[N], id1[N], id2[N], m, cnt, T;ll a[N], g[N], sum[N], f[N], n;bool st[N];int ID(ll x) {return x <= T ? id1[x] : id2[n / x];}void init() {cnt = m = 0;T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum[cnt] = sum[cnt - 1] + 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if(a[m] <= T) id1[a[m]] = m;else id2[n / a[m]] = m;g[m] = a[m] - 1;}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];}}for(int i = 1; i <= T; i++) {st[i] = 0;}}ll solve(ll n, int m) {if(n < prime[m]) return 0;ll ans = 1ll * prime[m - 1] * (g[ID(n)] - sum[m - 1]);for(int j = m; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {for(ll i = prime[j]; i * prime[j] <= n; i *= prime[j]) {ans += solve(n / i, j + 1) + prime[j];}}return ans;}ll solve(ll x) {if(x <= 1) return 0;n = x;init();return solve(x, 1);}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);ll l = read(), r = read();cout << Min_25::solve(r) - Min_25::solve(l - 1) << endl;return 0;
}

#572. 「LibreOJ Round #11」Misaka Network 与求和

推式子

∑i=1n∑j=1nf(gcd(i,j))k∑d=1nf(d)k∑i=1nd∑j=1nd[gcd(i,j)=1]∑d=1nf(d)k∑K=1ndμ(k)(nKd)2t=Kd∑t=1n(nt)2∑d∣tf(d)kμ(td)我们记f(x)k=F(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I，有F∗μ∗I=F∗ϵ=F得到后半部分的前缀和S(n)=∑i=1nF(i)−∑i=2nS(ni)\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i, j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^k \sum_{K = 1} ^{\frac{n}{d}} \mu(k) \left( \frac{n}{Kd} \right) ^2\\ t = Kd\\ \sum_{t = 1} ^{n} \left(\frac{n}{t} \right) ^ 2 \sum_{d \mid t} f(d) ^ k \mu(\frac{t}{d})\\ 我们记f(x) ^ k = F(x)\\ 上面式子后半部分是一个迪利克雷卷积形式:F * \mu\\ 所以我们卷上一个I，有F * \mu * I = F * \epsilon = F\\ 得到后半部分的前缀和S(n) = \sum_{i = 1} ^{n} F(i) - \sum_{i = 2} ^{n} S(\frac{n}{i})\\ i=1∑nj=1∑nf(gcd(i,j))kd=1∑nf(d)ki=1∑dnj=1∑dn[gcd(i,j)=1]d=1∑nf(d)kK=1∑dnμ(k)(Kdn)2t=Kdt=1∑n(tn)2d∣t∑f(d)kμ(dt)我们记f(x)k=F(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I，有F∗μ∗I=F∗ϵ=F得到后半部分的前缀和S(n)=i=1∑nF(i)−i=2∑nS(in)

化简到这里只需要跟上面一题类似用Min_25求∑i=1nF(i)\sum\limits_{i = 1} ^{n} F(i)i=1∑nF(i)，然后用杜教筛求S(n)S(n)S(n)即可得到答案。

代码

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1using namespace std;typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;inline ll read() {ll f = 1, x = 0;char c = getchar();while(c < '0' || c > '9') {if(c == '-')    f = -1;c = getchar();}while(c >= '0' && c <= '9') {x = (x << 1) + (x << 3) + (c ^ 48);c = getchar();}return f * x;
}#define uint unsigned intconst int N = 1e6 + 10;uint quick_pow(uint a, int n) {uint ans = 1;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}return ans;
}namespace Min_25 {uint prime[N], g[N], sum[N], f[N];int a[N], id1[N], id2[N], n, m, k, cnt, T;bool st[N];int ID(int x) {return x <= T ? id1[x] : id2[n / x];}void init() {cnt = m = 0;T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;f[cnt] = quick_pow(i, k);sum[cnt] = sum[cnt - 1] + 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if(a[m] <= T) id1[a[m]] = m;else id2[n / a[m]] = m;g[m] = a[m] - 1;}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];}}for(int i = 1; i <= T; i++) {st[i] = 0;}}uint solve(int n, int m) {if(n < prime[m] || n <= 1) return 0;uint ans = f[m - 1] * (g[ID(n)] - sum[m - 1]);for(int j = m; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {for(ll i = prime[j]; i * prime[j] <= n; i *= prime[j]) {ans += solve(n / i, j + 1) + f[j];}}return ans;}uint solve(int n) {if(n <= 1) return 0;return solve(n, 1) + g[ID(n)];}
}unordered_map<int, uint> ans_s;uint S(int n) {if(ans_s.count(n)) return ans_s[n];uint ans = Min_25::solve(n);for(uint l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans -= (r - l + 1) * S(n / l);}return ans_s[n] = ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);uint n = read(), k = read();Min_25::n = n, Min_25::k = k;Min_25::init();uint ans = 0;for(uint l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans += (n / l) * (n / l) * (S(r) - S(l - 1));}cout << ans << endl;return 0;
}