网站开发小组分工,常州网站建设公司效果,长沙网页制作开发公司,网络管理系统的基本组成和功能#572. 「LibreOJ Round #11」Misaka Network 与求和
推式子 ∑i1n∑j1nf(gcd(i,j))k∑d1nf(d)k∑i1nd∑j1nd[gcd(i,j)1]∑d1nf(d)k∑K1ndμ(k)(nKd)2tKd∑t1n(nt)2∑d∣tf(d)kμ(td)我们记f(x)kF(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I#xff0c…#572. 「LibreOJ Round #11」Misaka Network 与求和
推式子
∑i1n∑j1nf(gcd(i,j))k∑d1nf(d)k∑i1nd∑j1nd[gcd(i,j)1]∑d1nf(d)k∑K1ndμ(k)(nKd)2tKd∑t1n(nt)2∑d∣tf(d)kμ(td)我们记f(x)kF(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I有F∗μ∗IF∗ϵF得到后半部分的前缀和S(n)∑i1nF(i)−∑i2nS(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})\\ i1∑nj1∑nf(gcd(i,j))kd1∑nf(d)ki1∑dnj1∑dn[gcd(i,j)1]d1∑nf(d)kK1∑dnμ(k)(Kdn)2tKdt1∑n(tn)2d∣t∑f(d)kμ(dt)我们记f(x)kF(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I有F∗μ∗IF∗ϵF得到后半部分的前缀和S(n)i1∑nF(i)−i2∑nS(in)
化简到这里只需要跟上面一题类似用Min_25求∑i1nF(i)\sum\limits_{i 1} ^{n} F(i)i1∑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 pairint, 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_mapint, 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;
}
写法二
∑i1n∑j1nf(gcd(i,j))k∑d1nf(d)k∑i1nd∑j1nd[gcd(i,j)1]∑d1nf(d)k(2∑i1ndϕ(i)−1)\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}} \phi(i) - 1)\\ i1∑nj1∑nf(gcd(i,j))kd1∑nf(d)ki1∑dnj1∑dn[gcd(i,j)1]d1∑nf(d)k(2i1∑dnϕ(i)−1) 充分利用上面递推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 pairint, 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_mapint, 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;
}
