洛谷P3768 简单的数学题 杜教筛+莫比乌斯反演

---

洛谷P3768 简单的数学题

---

标签

• 莫比乌斯反演
• 狄利克雷卷积
• 杜教筛

  ---

前言

• 我的csdn和博客园是同步的，欢迎来访danzh-博客园~
• 很简单很好推~

  ---

简明题意

• 求
  $$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)(模p意义下)$$

  ---

思路

• 很简单鸭。顺着推一遍就出来了~
  $$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$$
• 改为枚举gcd：
  $$\sum _{x=1}^{n}x\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)==x]$$

都是套路啊！倍加式中出现$gcd(i,j)$，则改为枚举

• 见到$[gcd(i,j)==x]$再熟悉不过了，但我还是一步步的说吧。先更换枚举上限：
  $$\sum _{x=1}^n{x}^3\sum _{i=1}^{[\frac{n}{x}]}\sum _{j=1}^{[\frac{n}{x}]}ij[gcd(i,j)==1]$$

这里简单说一下为什么是$x^3$。因为更换枚举上限后，枚举的是就是原来的$\frac{1}{x}$，所以ij都变成了原来的$\frac{1}{x}$，所以应该乘上$x^2$，然后再移到前面就是$x^3$了

• 然后反演鸭，用${\sum }_{d|n}\mu (d)$替换$[gcd(i,j)==1]$，并且将$d|gcd(i,j)$改为$d|i且d|j$得到：
  $$\sum _{x=1}^n{x}^3\sum _{i=1}^{[\frac{n}{x}]}\sum _{j=1}^{[\frac{n}{x}]}ij\sum _{d|i且d|j}\mu (d)$$

记得关于gcd的性质：$d|gcd(i,j)⟺d|i且d|j$

• 然后就又是套路了，替换后就改为枚举d，那么就是：
  $$\sum _{x=1}^n{x}^3\sum _{d=1}^{[\frac{n}{x}]}\mu (d)\sum _{i=1}^{[\frac{n}{x}]}\sum _{j=1}^{[\frac{n}{x}]}ij[d|i且d|j]$$
• 然后看到$[d|i且d|j]$我们就又忍不住跟换枚举上限了，记得换了上限后，枚举的数成了原来的$\frac{1}{d}$，所以应该乘回去：
  $$\sum _{x=1}^n{x}^3\sum _{d=1}^{[\frac{n}{x}]}\mu (d)\sum _{i=1}^{[\frac{n}{dx}]}\sum _{j=1}^{[\frac{n}{dx}]}ij{d}^2$$
• 小学奥数告诉我们$g(x)=\sum _{i=1}^x\sum _{j=1}^{x}ij=x^2\ast (x+1{)}^2/4$，所以原式实际是求：
  $$\sum _{x=1}^n\sum _{d=1}^{[\frac{n}{x}]}\mu (d)g([\frac{n}{dx}])d^2{x}^3$$

上面的小学奥数式子是开玩笑的。实际上是这样推导的：$g(x)=\sum _{i=1}^x\sum _{j=1}^{x}ij=\sum _{i=1}^{x}i\sum _{i=1}^{x}i=(n\ast (n+1)/2{)}^2$

• 这个时候一看，x和d的上限的关系，就发现可以更换枚举项，我们改为枚举$K=dx$，然后原式就是：
  $$\sum _{k=1}^n{k}^{2}g([\frac{n}{k}])\sum _{d|k}\mu (\frac{k}{d})d$$

这一步如何改为枚举$K=dx$的，不清楚的同学请看我的博客《数论公式总结》戳这里

• 很显然$\sum _{d|k}\mu (\frac{k}{d})d$就是$\mu 和id$的狄利克雷卷积鸭！人尽皆知$\mu \ast id=\varphi$，所以我们实际求：
  $$\sum _{k=1}^n{k}^2\varphi (k)g([\frac{n}{k}])$$

• 这里一看g里面的$\frac{n}{k}$，就知道要分块了~所以哦我们需要求出$k^2\mu (k)$的前缀和，暴力显然是不行的，因此我们杜教筛呗，所以现在问题就是如何杜教筛

• 现在要求的是$x^2\varphi (x)$的前缀和，我们令$f(x)=x^2\varphi (x)$，再找一个函数$g(x)$，于是杜教筛一下：
  $$S(n)g(1)=\sum _{i=1}^{n}f\ast g-\sum _{i=2}^{n}g(i)S([\frac{n}{i}])$$

• 其实，当杜教筛找$g(x)$时，往往更注重的是$f\ast g$的狄利克雷卷积的前缀和好求，通常两者卷积是$ϵ$或者$id$或者$\varphi (n)$。因此现在需要找到一个$g(x)$，使得$f\ast g$是$ϵ$或者$id$。我们写出两者的卷积：
  $$(f\ast g)(n)=\sum _{d|n}f(d)g([\frac{n}{d}])$$

• 我们把$f(x)=x^2\varphi (x)$带进去，得出$f$和$g$的卷积是：
  $$(f\ast g)(n)=\sum _{d|n}{d}^2\varphi (d)g([\frac{n}{d}])$$

• 到了这里，是不是很容易确定$g$呢？如果令$g=i{d}^2$，会发现$d^2$被消掉了耶。那么$f$和$g$的卷积是：
  $$(f\ast g)(n)=\sum _{d|n}{d}^2\varphi (d)\frac{n}{d}\frac{n}{d}=\sum _{d|n}\varphi (d)n^2=n^2\sum _{d|n}\varphi (d)=n^3$$

$\sum _{d|n}\varphi (d)=n$，n得约数得欧拉函数之和是n，这是常用的结论！！！

• 把上面的卷积带入杜教筛中，可以得：
  $$S(n)=\sum _{i=1}^n{i}^3-\sum _{i=2}^{n}g(i)S([\frac{n}{i}])$$

• 由小学奥数可以$O(1)$求出$i^3$的前缀和，然后杜教筛写出来就A啦！！

---

注意事项

• 在计算$i^3,i^2$前缀和是会溢出，要算一下乘法逆元。而且，$i$可能大于模数，应该先取模再运算

  ---

总结

• 其实，当杜教筛找$g(x)$时，往往更注重的是$f\ast g$的狄利克雷卷积的前缀和好求，通常两者卷积是$ϵ$或者$id$。因此现在需要找到一个$g(x)$，使得$f\ast g$是$ϵ$或者$id$。常见的狄利克雷卷积见我的另一篇博客戳这里
• 总结几个求和公式：
  $$\sum _{i=1}^n{i}^3=\sum _{i=1}^n\sum _{j=1}^{n}ij={\left(\frac{n(n-1)}{2}\right)}^2$$
  $$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$

---

AC代码

#include<cstdio>
#include<unordered_map>
using namespace std;

const int maxn = 1e7 + 10;

long long n;
int mod;
unordered_map<long long, int> rec;

int ksm(int a, int b)
{
	int ans = 1, base = a % mod;
	while (b)
	{
		if (b & 1)
			ans = 1ll * ans * base % mod;
		base = 1ll * base * base % mod;
		b >>= 1;
	}
	return ans;
}

bool no_prime[maxn];
int prime[maxn], phi[maxn], pre_phi[maxn];
int shai(int n)
{
	int cnt = 0;
	phi[1] = 1;

	for (int i = 2; i <= n; i++)
	{
		if (!no_prime[i])
			prime[++cnt] = i, phi[i] = i - 1;

		for (int j = 1; j <= cnt && prime[j] * i <= n; j++)
		{
			no_prime[prime[j] * i] = 1;
			phi[prime[j] * i] = i % prime[j] == 0 ? phi[i] * prime[j] : phi[i] * (prime[j] - 1);
			if (i % prime[j] == 0) break;
		}
	}

	for (int i = 1; i <= n; i++)
		pre_phi[i] = (pre_phi[i - 1] + 1ll * phi[i] * i % mod * i) % mod;

	return cnt;
}

int inv2;
int inv6;

int pre_x2(long long n)
{
	n %= mod;
	return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}

int pre_x3(long long n)
{
	n %= mod;
	return (n * (n + 1) % mod * inv2 % mod) * (n * (n + 1) % mod * inv2 % mod) % mod;
}

int S(long long n)
{
	if (n <= maxn - 10) return pre_phi[n];
	if (rec[n]) return rec[n];
	int ans = pre_x3(n);
	
	long long l = 2, r;
	while (l <= n)
	{
		r = n / (n / l);
		ans = (ans - (pre_x2(r) - pre_x2(l - 1)) * 1ll * S(n / l)) % mod;
		l = r + 1;
	}
	return rec[n] = ans;
}

void solve()
{
	scanf("%d%lld", &mod, &n);
	shai(maxn - 10);
	inv2 = ksm(2, mod - 2);
	inv6 = ksm(6, mod - 2);

	long long l = 1, r;
	int ans = 0;
	while (l <= n)
	{
		r = n / (n / l);
		int gg = (S(r) - S(l - 1));
		ans = (ans + 1ll * (S(r) - S(l - 1)) * 1ll * pre_x3(n / l)) % mod;
		l = r + 1;
		//printf("%d\n", gg);
	}
	
	printf("%d", (ans % mod + mod) % mod);
}

int main()
{
	solve();
	return 0;
}