BZOJ3529: [Sdoi2014]数表(莫比乌斯反演 树状数组)
程序员文章站
2022-04-14 10:59:40
题意 "题目链接" Sol 首先不考虑$a$的限制 我们要求的是 $$\sum_{i = 1}^n \sum_{j = 1}^m \sigma(gcd(i, j))$$ 用常规的套路可以化到这个形式 $$\sum_{d = 1}^n \sigma (d) \sum_{k = 1}^{\frac{n} ......
题意
sol
首先不考虑\(a\)的限制
我们要求的是
\[\sum_{i = 1}^n \sum_{j = 1}^m \sigma(gcd(i, j))\]
用常规的套路可以化到这个形式
\[\sum_{d = 1}^n \sigma (d) \sum_{k = 1}^{\frac{n}{d}} \mu(k) \frac{n}{kd} \frac{m}{kd}\]
设\(kd = t\)
那么
\(\sum_{t = 1}^n \left\lfloor \frac{n}{t} \right\rfloor \left\lfloor \frac{m}{t} \right\rfloor \sum_{d \ | t} \sigma(d) \mu(\frac{t}{d})\)
然后按照套路筛出后面的卷积。
但是现在有\(a\)的限制了。。
那么可以直接离线之后树状数组维护一下。
时间复杂度:\(o(t\sqrt{n}logn)\)
约数个数和用埃氏筛两行就筛出来了
#include<bits/stdc++.h> #define chmax(a, b) (a = (a > b ? a : b)) #define chmin(a, b) (a = (a < b ? a : b)) //#define int long long using namespace std; const int maxn = 3e5 + 10; inline int read() { char c = getchar(); int x = 0, f = 1; while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();} while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar(); return x * f; } int t, mx, mu[maxn], phi[maxn], prime[maxn], tot, vis[maxn], si[maxn], out[maxn], p[maxn]; struct query{ int n, m, a, id; bool operator < (const query &rhs) const { return a < rhs.a; } }q[maxn]; int comp(const query &a, const query &b) { return a.id < b.id; } int comp2(const int &a, const int &b) { return si[a] < si[b]; } void get(int n) { for(int i = 1; i <= n; i++) for(int j = i; j <= n; j += i) si[j] += i; vis[1] = 1; mu[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) prime[++tot] = i, mu[i] = -1; for(int j = 1; j <= tot && i * prime[j] <= n; j++) { vis[i * prime[j]] = 1; if(!(i % prime[j])) {mu[i * prime[j]] = 0; break;} else mu[i * prime[j]] = -mu[i]; } } } #define lb(x) (x & (-x)) int tr[maxn]; void add(int x, int val) { while(x <= mx) tr[x] += val, x += lb(x); } int sum(int x) { int ans = 0; while(x) ans += tr[x], x -= lb(x); return ans; } int solve(int n, int m) { int rt = 0, las = 0, now; for(int i = 1, nxt; i <= n; i = nxt + 1) { nxt = min(m / (m / i), n / (n / i)); now = sum(nxt); rt += (n / i) * (m / i) * (now - las); las = now; } return rt; } signed main() { t = read(); for(int i = 1; i <= t; i++) { q[i].n = read(), q[i].m = read(), q[i].a = read(), q[i].id = i; if(q[i].n > q[i].m) swap(q[i].n, q[i].m);; chmax(mx, q[i].n); } get(mx); for(int i = 1; i <= mx; i++) p[i] = i; sort(p + 1, p + mx + 1, comp2); sort(q + 1, q + t + 1); for(int t = 1, d = 1; t <= t; t++) { for(; d <= mx && si[p[d]] <= q[t].a; d++) { for(int k = p[d]; k <= mx; k += p[d]) { if(!mu[k / p[d]]) continue; add(k, si[p[d]] * mu[k / p[d]]); } } out[q[t].id] = solve(q[t].n, q[t].m); } for(int i = 1; i <= t; i++) printf("%d\n", out[i] < 0 ? out[i] + 2147483648 : out[i]); //printf("%d\n", out[i] <); return 0; }