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

推式子

$$\begin{gathered} \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??=F \\ 得到后半部分的前綴和S(n)=\sum _{i=1}^{n}F(i)?\sum _{i=2}^{n}S(\frac{n}{i}) \end{gathered}$$

化簡到這里只需要跟上面一題類似用Min_25求$\sum _{i=1}^{n}F(i)$，然后用杜教篩求$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], calc[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;}/*非遞歸版本*/// for(int j = cnt; j >= 1; j--) {// for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {// for(ll k = prime[j]; k * prime[j] <= a[i]; k *= prime[j]) {// calc[i] += calc[ID(a[i] / k)] + (g[ID(a[i] / k)] - sum[j - 1]) * f[j];// }// }// }}// uint solve(int x) {// if(x <= 1) return 0;// return calc[ID(x)] + g[ID(x)];// }/*下面是遞歸版本*/// uint solve(int n, int m) {// if(n < prime[m] || n <= 1) return 0;// uint ans = 0;// 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) + (g[ID(n / i)] - sum[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; }

寫法二

$$\begin{gathered} \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(2\sum _{i=1}^{\frac{n}{d}}?(i)?1) \end{gathered}$$
充分利用上面遞推Min_25，得到一個復雜度更優的算法。

/*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], calc[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 = 1000000;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;}for(int j = cnt; j >= 1; j--) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {for(ll k = prime[j]; k * prime[j] <= a[i]; k *= prime[j]) {calc[i] += calc[ID(a[i] / k)] + (g[ID(a[i] / k)] - sum[j - 1]) * f[j];}}}}uint solve(int x) {if(x <= 1) return 0;return calc[ID(x)] + g[ID(x)];} }namespace Djs {uint prime[N], phi[N], cnt;bool st[N];void init() {phi[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[++cnt] = i;phi[i] = i - 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {phi[i * prime[j]] = phi[i] * prime[j]; break;}phi[i * prime[j]] = phi[i] * (prime[j] - 1);}}for(int i = 1; i < N; i++) {phi[i] += phi[i - 1];}}unordered_map<int, uint> ans_s;uint S(int n) {if(n < N) return phi[n];if(ans_s.count(n)) return ans_s[n];uint ans = 1ll * n * (n + 1) / 2;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();Djs::init();uint ans = 0;for(uint l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans += (2 * Djs::S(n / l) - 1) * (Min_25::solve(r) - Min_25::solve(l - 1));}cout << ans << endl;return 0; } 創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的「LibreOJ Round #11」Misaka Network 与求和（杜教筛 + Min_25）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。