洛谷P3312 [SDOI2014]数表
标签
- 莫比乌斯反演
- 线段树
- 约数和筛
前言
- 用数据结构维护的数论题..有点难，不过做完后收获很大~
简明题意
- 有一个\(n *m\)的数表，每个点\((i,j)\)的值是能同时整除\(i,j\)的自然数之和。再给定\(a\)，需要你求表中所有值\(>a\)的项之和。
思路
- 首先用数学公式表示出来：
\[\sum_{i=1}^n\sum_{j=1}^m(\sum_{d|i且d|j}d)*[\sum_{d|i且d|j}d<=a]\] - 我们知道关于\(gcd\)有一个很常用的性质：\(d|gcd(i,j) \iff d|i且d|j\)（这个太常用了，不会的同学一定要自行搞清楚）。然后我们替换掉，就成了：
\[\sum_{i=1}^n\sum_{j=1}^m(\sum_{d|gcd(i,j)}d)*[\sum_{d|gcd(i,j)}d<=a]\] 然后就是改为枚举\(gcd(i,j)\)，显然上限是\(n\)
\[\sum_{x=1}^n\sum_{i=1}^n\sum_{j=1}^m\sum_{d|x}d*[\sum_{d|x}d<=a][gcd(i,j)==x]\]式子反而变复杂了？不要担心，通常改为枚举会使得式子最前面多一个\(\sum\)，被加式中多一个条件。但后面的步骤就会把他们都简化
- 接下来移项，成为：
\[\sum_{x=1}^n\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==x]\sum_{d|x}d*[\sum_{d|x}d<=a]\] 然后由莫比乌斯反演我们又有这样一个式子：
\[\sum_{i=1}^{n}\sum_{j=1}^m[gcd(i,j)==x]=\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\]这个式子很常用。不清楚的同学参见我的另一篇博客戳这里
- 把这个式子代进去，原式就是：
\[\sum_{x=1}^n\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\sum_{d|x}d*[\sum_{d|x}d<=a]\] - 然后观察这个式子：
\[\sum_{d|x}d\] - 数论比较好的同学可以立马发现这就是\(\sigma\)函数，方便起见，我们就用\(\sigma\)替换
\[\sum_{x=1}^n\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\sigma(x)*[\sigma(x)<=a]\] - 然后式子比较长，我们令\(f(x)=\sigma(x)*[\sigma(x)<=a]\)：
\[\sum_{x=1}^nf(x)\sum_{d=1}^{[\frac nx]}\left( \mu(d)*[\frac {n}{dx}]*[\frac {m}{dx}]\right)\] - 到了这里，我们可以用线性筛预处理出\(f,\mu\)函数（不会约数个数和筛的同学请看看别的博客弄明白）,再前缀和一下，后面的分以下块，整体复杂度就会是\(O(n+Tn\sqrt{n})\)，可以拿到60分。我们进一步优化。
这里我们令\(k=dx\)，然后改为枚举\(k\)：
\[\sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]\sum_{d|k}f(d)\mu([\frac kd])\]这一步一定很多同学会产生疑惑。我在这里解释一下，我们有这样一个结论：
\[\sum_{i=1}^n\sum_{j=1}^{[\frac ni]}f(ij)*g(i)*h(j)=\sum_{k=1}^nf(k)*\sum_{d|k}g([\frac kd])h(d)\]
只要是二重枚举式，第一重枚举随意，且第二重枚举枚举的是第一重的上限除以第一重的枚举项,就可以将被加式写成关于\(i，j，i*j\)的三个函数，所以可以令\(k=ij\)。这个式子要牢记。- 现在令\(g(k) =\sum\limits_{d|k}f(d)\mu([\frac kd])\),原式就是：
\[\sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]g(k)\] - 这个式子，我们十分想对它进行整除分块，如果整除分块，显然应该预处理出\(g(k)\)的前缀和。然而，每组样例\(a\)不同，会影响\(g(k)\)的取值，因此无法对\(g(k)\)前缀和。怎么办呢？所以现在的关键就是如何处理对不同的\(a\)计算。现在先把\(f(d)\)代回去：
\[g(k) =\sum\limits_{d|k}\sigma(d)*[\sigma(d)<=a]\mu([\frac kd]\] - 我们先不考虑\(g(k)\)的前缀和，而是考虑单个的\(g(k)\)。实际上，\(g(k)\)的取值并不是只和\(k\)有关，还与\(a\)有关，严谨一点，\(g(k)\)应该是一个二元函数\(g(k,a)\)。但是实际上\(k\)并不重要，因为每一次分块都要求\(g(k)\)的前缀和，所以实际上对于给定的\(a\)，必须求出所有的\(g(k)\)，而不是只用求出单独一项\(g(k)\)。所以，关键的是\(a\)的取值。这是我们换一个角度思考这个问题
- 对于给定的\(k\)，显然\(a\)越大，就会有更多的\(\sigma\)贡献到\(g(k)\)。如果\(a\)从\(3\)变化到\(10\)，对于同一个\(k\)，会有\(3<\sigma(d)<=10的\sigma(d)*\mu(\frac kd)(d|k)\)被加入到了原先的\(g(k)\)。（多读几遍黑色的部分），是不是注意到了\(a\)越大，\(g(k)\)越大，而且更大的\(a\)得到的\(g(k)\)可以由更小的\(a\)推过来。
- 所以这里，我们先离线一下，把\(a\)从小到大排序。然后对\(g(k)\)建一棵线段树，每一次\(a\)增大，更新一下线段树。具体更新的方法是：把位于\((上一次的a，当前a]\)范围内的\(\sigma(d)\)找出来，符合要求的每一个\(\sigma(d)\)都会对\(d|k\)的\(g(k)\)产生贡献。就对\(d\)枚举它的所有倍数\(k\)，让其值加上\(\sigma(d)*\mu(\frac kd)\)。每一次修改操作进行完，\(g(k)\)就是正确的，所以回想刚刚要求的式子：
\[\sum_{k=1}^n[\frac {n}{k}][\frac {m}{k}]g(k)\] 我们直接数论分块，每一块求出区间\([l,r]\)的\(g(k)\)之和乘以\([\frac nl][\frac mk]\)，全部累加起来就是这一次的答案。
注意事项
- 无
总结
- 对\(2^{31}\)取模等价于最终的答案&\(2^{31}-1\)
- 前缀和如果动态变化，可以尝试用线段树维护。
- \[\sum_{i=1}^n\sum_{j=1}^{[\frac ni]}f(ij)*g(i)*h(j)=\sum_{k=1}^nf(k)*\sum_{d|k}g([\frac kd])h(d)\]
这个式子很常用！ - 注意约数和筛的写法！
AC代码
#include<cstdio>
#include<algorithm>
using namespace std;
const int maxn = 1e5 + 10;
struct Sig
{
int sigma, id;
bool operator < (const Sig& a)const
{
return sigma < a.sigma;
}
};
Sig sigma[maxn];
bool no_prime[maxn];
int prime[maxn], s[maxn], ssum[maxn], mu[maxn], pre_mu[maxn];
int shai(int n)
{
int cnt = 0;
mu[1] = 1, s[1] = 1;
for (int i = 2; i <= n; i++)
{
if (!no_prime[i])
prime[++cnt] = i, mu[i] = -1, s[i] = i + 1, ssum[i] = i + 1;
for (int j = 1; j <= cnt && prime[j] * i <= n; j++)
{
no_prime[prime[j] * i] = 1;
mu[prime[j] * i] = i % prime[j] == 0 ? 0 : -mu[i];
s[prime[j] * i] = i % prime[j] == 0 ? s[i] / ssum[i] + s[i] * prime[j] : s[i] * (1 + prime[j]);
ssum[prime[j] * i] = i % prime[j] == 0 ? ssum[i] * prime[j] + 1 : prime[j] + 1;
if (i % prime[j] == 0) break;
}
}
for (int i = 1; i <= n; i++)
pre_mu[i] = pre_mu[i - 1] + mu[i], sigma[i].sigma = s[i], sigma[i].id = i;
sort(sigma + 1, sigma + 1 + n);
return cnt;
}
int n, m, a;
int f(int x)
{
if (s[x] <= a) return s[x];
return 0;
}
struct Query
{
int n, m, a, id;
bool operator < (const Query& b)const
{
return a < b.a;
}
};
Query query[maxn];
struct Node
{
int l, r, sum;
};
Node tree[maxn * 4];
void build(int o, int l, int r)
{
tree[o].l = l, tree[o].r = r;
if (l == r)
return;
int mid = (l + r) / 2;
build(o * 2, l, mid);
build(o * 2 + 1, mid + 1, r);
}
void change(int o, int x, int k)
{
if (tree[o].l == tree[o].r)
{
tree[o].sum += k;
return;
}
int mid = (tree[o].l + tree[o].r) / 2;
if (x <= mid)
change(o * 2, x, k);
else
change(o * 2 + 1, x, k);
tree[o].sum = (tree[o * 2].sum + tree[o * 2 + 1].sum);
}
int ask(int o, int l, int r)
{
if (tree[o].l == l && tree[o].r == r)
return tree[o].sum;
int mid = (tree[o].l + tree[o].r) / 2;
if (r <= mid)
return ask(o * 2, l, r);
else if (l > mid)
return ask(o * 2 + 1, l, r);
return ask(o * 2, l, mid) + ask(o * 2 + 1, mid + 1, r);
}
int cal(int n, int m)
{
int l = 1, r, ans = 0;
while (l <= n)
{
r = min(n / (n / l), m / (m / l));
ans += (n / l) * (m / l) * ask(1, l, r);
l = r + 1;
}
return ans;
}
int ans[maxn];
void solve()
{
shai(maxn - 10);
build(1, 1, maxn - 10);
int t;
scanf("%d", &t);
for (int i = 1; i <= t; i++)
{
scanf("%d%d%d", &query[i].n, &query[i].m, &query[i].a);
query[i].id = i;
if (query[i].n > query[i].m) swap(query[i].n, query[i].m);
}
sort(query + 1, query + 1 + t);
int pt = 1;
for (int i = 1; i <= t; i++)//i表示第i组样例
{
for (; sigma[pt].sigma <= query[i].a; pt++)//处理每一个sigma
{
for (int k = sigma[pt].id; k <= maxn - 10; k += sigma[pt].id)
change(1, k, sigma[pt].sigma * mu[k / sigma[pt].id]);
}
ans[query[i].id] = cal(query[i].n, query[i].m);
}
for (int i = 1; i <= t; i++)
printf("%d\n", ans[i] & 2147483647);
}
int main()
{
solve();
return 0;
}
扫一扫
5万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
