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

洛谷P3768 简单的数学题

标签

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

前言

• 很简单很好推~

简明题意

• 求

\[\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)(模p意义下) \]

思路

• 很简单鸭。顺着推一遍就出来了~

\[\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j) \]

• 改为枚举gcd：

\[\sum_{x=1}^nx\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)==x] \]

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

• 见到\([gcd(i,j)==x]\)再熟悉不过了，但我还是一步步的说吧。先更换枚举上限：

\[\sum_{x=1}^nx^3\sum_{i=1}^{[\frac nx]}\sum_{j=1}^{[\frac nx]}ij[gcd(i,j)==1] \]

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

• 然后反演鸭，用\(\sum_{d|n}\mu(d)\)替换\([gcd(i,j)==1]\)，并且将\(d|gcd(i,j)\)改为\(d|i且d|j\)得到：

\[\sum_{x=1}^nx^3\sum_{i=1}^{[\frac nx]}\sum_{j=1}^{[\frac nx]}ij\sum_{d|i且d|j}\mu(d) \]

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

• 然后就又是套路了，替换后就改为枚举d，那么就是：

\[\sum_{x=1}^nx^3\sum_{d=1}^{[\frac nx]}\mu(d)\sum_{i=1}^{[\frac nx]}\sum_{j=1}^{[\frac nx]}ij[d|i且d|j] \]

• 然后看到\([d|i且d|j]\)我们就又忍不住跟换枚举上限了，记得换了上限后，枚举的数成了原来的\(\frac 1d\)，所以应该乘回去：

\[\sum_{x=1}^nx^3\sum_{d=1}^{[\frac nx]}\mu(d)\sum_{i=1}^{[\frac n{dx}]}\sum_{j=1}^{[\frac n{dx}]}ijd^2 \]

• 小学奥数告诉我们\(g(x)=\sum\limits_{i=1}^x\sum\limits_{j=1}^xij=x^2*(x+1)^2/4\)，所以原式实际是求：

\[\sum_{x=1}^n\sum_{d=1}^{[\frac nx]}\mu(d)g([\frac n {dx}])d^2x^3 \]

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

• 这个时候一看，x和d的上限的关系，就发现可以更换枚举项，我们改为枚举\(K=dx\)，然后原式就是：

\[\sum_{k=1}^nk^2g([\frac nk])\sum_{d|k}\mu(\frac kd)d \]

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

• 很显然\(\sum\limits_{d|k}\mu(\frac kd)d\)就是\(\mu和id\)的狄利克雷卷积鸭！人尽皆知\(\mu*id=\phi\)，所以我们实际求：

\[\sum_{k=1}^nk^2\phi(k)g([\frac nk]) \]

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

• 现在要求的是\(x^2\phi(x)\)的前缀和，我们令\(f(x)=x^2\phi(x)\)，再找一个函数\(g(x)\)，于是杜教筛一下：

\[S(n)g(1)=\sum_{i=1}^nf*g-\sum_{i=2}^ng(i)S([\frac ni]) \]

• 其实，当杜教筛找\(g(x)\)时，往往更注重的是\(f*g\)的狄利克雷卷积的前缀和好求，通常两者卷积是\(\epsilon\)或者\(id\)或者\(\phi(n)\)。因此现在需要找到一个\(g(x)\)，使得\(f*g\)是\(\epsilon\)或者\(id\)。我们写出两者的卷积：

\[(f*g)(n)=\sum_{d|n}f(d)g([\frac nd]) \]

• 我们把\(f(x)=x^2\phi(x)\)带进去，得出\(f\)和\(g\)的卷积是：

\[(f*g)(n)=\sum_{d|n}d^2\phi(d)g([\frac nd]) \]

• 到了这里，是不是很容易确定\(g\)呢？如果令\(g=id\)，会发现\(d^2\)被消掉了耶。那么\(f\)和\(g\)的卷积是：

\[(f*g)(n)=\sum_{d|n}d^2\phi(d)\frac nd\frac nd=\sum_{d|n}\phi(d)n^2=n^2\sum_{d|n}\phi(d)=n^3 \]

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

• 把上面的卷积带入杜教筛中，可以得：

\[S(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ng(i)S([\frac ni]) \]

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

注意事项

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

总结

• 其实，当杜教筛找\(g(x)\)时，往往更注重的是\(f*g\)的狄利克雷卷积的前缀和好求，通常两者卷积是\(\epsilon\)或者\(id\)。因此现在需要找到一个\(g(x)\)，使得\(f*g\)是\(\epsilon\)或者\(id\)。常见的狄利克雷卷积见我的另一篇博客戳这里
• 总结几个求和公式：

\[\sum_{i=1}^ni^3=\sum_{i=1}^n\sum_{j=1}^nij=\left( \frac {n(n-1)}2 \right)^2 \]

\[\sum_{i=1}^ni^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;
}