bzoj [SDOI2014]数表 莫比乌斯反演 BIT
链接
思路
[sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}a*[f[gcd(i,j)]<=a]
]
[f[]可以O(n)预处理出来
]
[sumlimits_{k=1}^{n}f[k]*sumlimits_{i=1}^{m}sumlimits_{j=1}^{m}[gcd(i,j)==k]
]
[sumlimits_{k=1}^{n}f[k]*sumlimits_{i=1}^{frac{n}{k}}{frac{n}{ki}}{frac{m}{ki}}mu(i)
]
[sumlimits_{i=1}^{n}f[i]*sumlimits_{k|d}{frac{n}{d}frac{m}{d}}mu(frac{d}{i})
]
d替换k*i
[sumlimits_{d=1}^{n} frac{n}{d} frac{m}{d} sumlimits_{k|d} mu(frac{d}{k})f(k)
]
[sumlimits_{d=1}^{n} frac{n}{d} frac{m}{d} g(d)
]
[g(d)=sumlimits_{k|d} mu(frac{d}{k})f(k)
]
询问按照a排序，每次加入f(k)时候影响的只是能k|d的g(d)
每次修改就是(O(sqrt{n}logn))
查询也是一样，
总的就是(O(Tsqrt{n}logn))
注意，ll+取模的话，loj会超时，用int的自然溢出就快了三倍(300ms)，是int，不是unsigned int。
其他
线性筛约数和
[x=p_{1}^{w_1}p_{2}^{w_2}…p_{k}^{w_k}
]
那么
[SD(x)=约数和=(1+p_1^1+p_1^2+…+p_1^{w_1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k})
]
0x00 是个素数
显然(SD(pri)=pri+1)
0x01 两两互质
是个积性函数
因为x,y两两互质，所以他们质因子互不相交，所以显然脑补公式
[SD(x*y)=SD(x)*SD(y)(gcd(x,y)==1)
]
0x02 两两不互质（i%pri[j]!=0）
再开个数组tmp，记录最小质因子因子的贡献((1+p_1^1+p_1^2+…+p_1^{w_1}))
因为pri[j]是他的最小质因子（因为这是线性筛）
我们之前求出的i的
[SD(i)=(1+p_1^1+p_1^2+…+p_1^{w_1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k})
]
现在的(i*pri[j])的SD显然就是
[SD(i*pri[j])=(1+p_1^1+p_1^2+…+p_1^{w_1}+p_1^{w_1+1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k})
]
改变的只有最小因子的贡献，tmp的作用就来了
tmp[i]我们已经求出来了，那么
[tmp[i*pri[j]]=tmp[i]*pri[j]+1（这是个等差数列）
]
[SD[i*pri[j]]=SD[i]/tmp[i]*tmp[i*pri[j]]
]
代码
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 6;
int read() {
int x = 0, f = 1; char s = getchar();
for (; s > '9' || s < '0'; s = getchar()) if (s == '-') f = -1;
for (; s >= '0' && s <= '9'; s = getchar()) x = x * 10 + s - '0';
return x * f;
}
int pri[N], tot, vis[N], mu[N];
int f[N], tmp[N];
struct node {
int n, m, a, id;
int ans;
bool operator < (const node &b) const {
return a < b.a;
}
} Q[N];
bool cmp(node a, node b) {
return a.id < b.id;
}
pair<int,int> F[N];
void Euler(int limit) {
f[1] = tmp[1] = mu[1] = 1;
for (int i = 2; i <= limit; ++i) {
if (!vis[i]) {
mu[i] = -1;
pri[++tot] = i;
f[i] = i + 1;
tmp[i] = i + 1;
}
for (int j = 1; j <= tot && i * pri[j] <= limit; ++j) {
vis[i * pri[j]] = 1;
if (i % pri[j] == 0) {
tmp[i * pri[j]] = tmp[i] * pri[j] + 1;
f[i * pri[j]] = f[i] / tmp[i] * tmp[i * pri[j]];
mu[i * pri[j]] = 0;
break;
}
mu[i * pri[j]] = -mu[i];
f[i * pri[j]] = f[i] * f[pri[j]];
tmp[i * pri[j]] = pri[j] + 1;
}
}
for (int i = 1; i <= limit; ++i) {
F[i].first = f[i], F[i].second = i;
}
sort(F + 1, F + 1 + limit);
}
namespace BIT {
int sum[N];
int lowbit(int x) {return x & (-x);}
void add(int x, int ad) {
for (int i = x; i <= 100000; i += lowbit(i)) sum[i] = (sum[i] + ad);
}
int query(int x) {
int ans = 0;
for (int i = x; i >= 1; i -= lowbit(i)) ans = (ans + sum[i]);
return ans;
}
}
int main() {
Euler(100000);
int 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;
}
sort(Q + 1, Q + 1 + T);
int now = 0;
for (int i = 1; i <= T; ++i) {
while (now + 1 <= 100000 && Q[i].a >= F[now + 1].first) {
now++;
for (int j = F[now].second; j <= 100000; j += F[now].second) {
BIT::add(j, mu[j / F[now].second] * F[now].first);
}
}
int ans = 0;
if (Q[i].n > Q[i].m) swap(Q[i].n, Q[i].m);
for (int l = 1, r; l <= Q[i].n; l = r + 1) {
r = min(Q[i].n / (Q[i].n / l), Q[i].m / (Q[i].m / l));
ans += 1LL * (Q[i].n / l) * (Q[i].m / l) * (BIT::query(r) - BIT::query(l - 1));
}
Q[i].ans = ans;
}
sort(Q + 1, Q + 1 + T, cmp);
for (int i = 1; i <= T; ++i) printf("%u
", Q[i].ans < 0 ? Q[i].ans + 2147483648 : Q[i].ans);
return 0;
}