bzoj3529 - 数表 (莫比乌斯函数,树状数组)
题目
有一张 n×m 的数表,其第 i 行第 j 列(1 <= i <= n, 1 <= j <= m)的数值为
能同时整除 i 和 j 的所有自然数之和。给定 a , 计算数表中不大于 a 的数之和。
题解
先不考虑不大于a这一限制条件。
令\(f(x)\)为x的约数和,那么数表(i,j)的值即为\(f(\gcd(i, j))\)。
开始的处理思路类似YY的gcd
\[\sum\limits^n_{i=1}{\sum\limits^m_{j=1}{f(\gcd(i, j))}}
\]
\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits^{\frac{n}{d_1}}_{i=1}{\sum\limits^{\frac{m}{d_1}}_{j=1}{[gcd(i,j)=1]f(d_1)}}}
\]
\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits^{\frac{n}{d_1}}_{i=1}{\sum\limits^{\frac{m}{d_1}}_{j=1}{\sum\limits_{d_2|\gcd(i,j)}{\mu(d_2)}f(d_1)}}}
\]
\[\sum\limits^{\min(n,m)}_{d_1}{\sum\limits_{d_2}{\mu(d_2)f(d_1) \lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor}}
\]
令\(T=d_1d_2\)
\[\sum\limits^{min(n, m)}_{T=1}{\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}}
\]
\[\sum\limits^{min(n, m)}_{T=1}{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)}}
\]
令\(F(T)=\sum\limits_{d|T}{\mu(\frac{T}{d})f(d)}\)。注意,式子中的f(d)若大于a则为零。如果我们能快速求出在限制下,\(F(T)\)的前缀和,那么就可以以\(O(\sqrt{n})\)的复杂度处理每个询问。
这里将询问按a从小到大排序,将f(d)从小到大排序。若当前f(i)小于等于a,则将其贡献加入(将其加入每个i的倍数的F中)。由于a的值单调,一共加n次,每次加\(\frac{n}{i}\),复杂度\(O(nln(n))\)。然后用树状数组维护前缀和。时间复杂度\(O(nln(n)+Q*\sqrt{n})\)。
#include <iostream>
#include <map>
#include <algorithm>
#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N)
typedef long long ll;
using namespace std;
/*-----------------------------------------------------------------*/
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f
const int N = 3e5 + 10;
const double eps = 1e-5;
const ll M = (1ll<<31) - 1;
bool vis[N];
ll prime[N];
ll cnt;
ll g[N];
ll sum[N];
struct node {
ll v;
int p;
bool operator<(const node&rhs) const {
return v < rhs.v;
}
} f[N];
node ask[N];
node nm[N];
ll tarr[N];
int mu[N];
ll ans[N];
void init() {
mu[1] = sum[1] = 1;
for(ll i = 2; i < N; i++) {
if(!vis[i]) {
mu[i] = -1;
g[i] = i + 1;
sum[i] = g[i];
vis[i] = 1;
prime[cnt++] = i;
}
for(ll j = 0; j < cnt && prime[j] * i < N; j++) {
ll p = prime[j];
vis[p * i] = 1;
if(i % p == 0) {
g[p * i] = g[i] * p + 1;
sum[p * i] = sum[i] / g[i] * g[p * i];
mu[p * i] = 0;
break;
} else {
sum[p * i] = sum[p] * sum[i];
g[p * i] = p + 1;
mu[p * i] = -mu[i];
}
}
}
for(int i = 1; i < N; i++) {
f[i].p = i;
f[i].v = sum[i];
}
sort(f + 1, f + N);
}
int lowbit(int x) {
return x&-x;
}
void add(int p, ll val) {
while(p < N) {
tarr[p] += val;
p += lowbit(p);
}
}
ll get(int p) {
ll res = 0;
while(p) {
res += tarr[p];
p -= lowbit(p);
}
return res;
}
int main() {
IOS;
init();
int q;
cin >> q;
for(int i = 1; i <= q; i++) {
ll n, m, a;
cin >> n >> m >> a;
ask[i] = node{a, i};
nm[i] = node{n, m};
}
sort(ask + 1, ask + 1 + q);
int p = 1;
for(int i = 1; i <= q; i++) {
while(p < N && f[p].v <= ask[i].v) {
for(int i = f[p].p; i < N; i += f[p].p) {
add(i, f[p].v * mu[i / f[p].p]);
}
p++;
}
ll res = 0;
int cur = 1;
int n = nm[ask[i].p].v, m = nm[ask[i].p].p;
while(cur <= min(n, m)) {
int p1 = n / (n / cur);
int p2 = m / (m / cur);
int nt = min(min(p1, p2), min(n, m));
res += 1ll * (get(nt) - get(cur - 1)) * (n / cur) * (m / cur);
cur = nt + 1;
}
ans[ask[i].p] = (res & M);
}
for(int i = 1; i <= q; i++) cout << ans[i] << endl;
}