洛谷P3312 [SDOI2014]数表 莫比乌斯反演+线段树+约数和筛

洛谷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;
}