P3768 简单的数学题 解题报告（莫比乌斯反演+杜教筛）

题目描述

题面简单，我打下：$$\sum _{i=1}^n\sum _{j=1}^{n}ij\gcd (i,j)\mathrm{m}\mathrm{o}\mathrm{d}p$$
n最大1e10，p为质数

思路分析

拿到这种题，肯定是先推式子，暴力跑不过去的。
    ${\sum }_{i=1}^n{\sum }_{j=1}^{n}ij\gcd (i,j)\mathrm{m}\mathrm{o}\mathrm{d}p$
$={\sum }_{d=1}^n{d}^3{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }[\gcd (i,j)=1]$(改变枚举顺序)
$={\sum }_{d=1}^n{d}^3{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }i{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }j{\sum }_{x|\gcd (i,j)}\mu (x)$(反演结论)
$={\sum }_{d=1}^n{d}^3{\sum }_{x=1}^{\lfloor \frac{n}{d}\rfloor }{x}^2\mu (x){\sum }_{i=1}^{\lfloor \frac{n}{xd}\rfloor }i{\sum }_{i=1}^{\lfloor \frac{n}{xd}\rfloor }j$(改变枚举顺序)

令sum(i)=$\frac{(1+i)i}{2}$，则，原式化为

$={\sum }_{d=1}^n{d}^3{\sum }_{x=1}^{\lfloor \frac{n}{d}\rfloor }\mu (x)x^{2}su{m}^2(\lfloor \frac{n}{xd}\rfloor )$
令T=xd
$={\sum }_{T=1}^n{T}^{2}su{m}^2(\lfloor \frac{n}{T}\rfloor ){\sum }_{d|T}\mu (d)\frac{T}{d}$
到此为止光从式子就推不下去了，但是，联想莫比乌斯函数的性质$$(\mu \ast id)=\varphi$$
也就是说，莫比乌斯函数和单位函数的迪利克雷卷积为欧拉函数。因此，原式可写为
$$\sum _{T=1}^{n}su{m}^2(\lfloor \frac{n}{T}\rfloor )T^2\varphi (T)$$
前面的sum函数部分可以用分块整除，那么，就需要考虑求出后面两个因式的部分和。
由于函数f(T)=T和欧拉函数都为积性函数，则，后面两个因式的乘积也为积性函数，再看看数据范围，故考虑杜教筛。
拿出杜教筛的套路式子：
$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )（答案式）$$
其中
$h(i)=(f\ast g)(i)$
$S(n)={\sum }_{i=1}^{n}f(i)$
若能得出上式就能递归求解了。
重点在于g(x)的选择，我们看下我们要处理的式子：$f(x)=x^2\varphi (x)$，带入上式得到：
$$g(1)S(n)=\sum _{i=1}^n\sum _{d|i}g(i)f(\frac{i}{d})-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )（1）$$
单独看左边的一项：
    ${\sum }_{i=1}^n{\sum }_{d|i}g(d)f(\frac{i}{d})$
$={\sum }_{i=1}^n{\sum }_{d|i}g(d)(\frac{i}{d}{)}^2\varphi (\frac{i}{d})（2）$
联想欧拉函数性质：${\sum }_{d|n}\varphi (d)=n$
因此我们考虑凑出上面的式子，也就是取$g(x)=x^2$
我们回代：
对于（2）式：
原式=${\sum }_{i=1}^n{i}^2{\sum }_{d|i}\varphi (d)={\sum }_{i=1}^n{i}^3$
带入（1）式：
$$S(n)=\sum _{i=1}^n{i}^3-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )（g(1)=1）$$
根据立方前缀和公式：
$$\sum _{i=1}^n{i}^3=(1+2+\dots \dots +n{)}^2=su{m}^2(n)$$
就能杜教筛求解上面的式子了。$S(\lfloor \frac{n}{d}\rfloor )$整除分块+递归求解+线性筛预处理解决，g(d)的部分和，由于$g(x)=x^2$，可以通过平方前缀和公式：
$$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$
转化为前缀和之差。
然后回代到上面的答案式就行了。
有点复杂的一题Orz。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int maxn = 2e6 + 10;
template<typename T>
void read(T&x){
	x=0;
	int f=1;
	char ch=getchar();
	while(!isdigit(ch)){
		if(ch=='-')f*=-1;
		ch=getchar();
	}
	while(isdigit(ch)){
		x=x*10+ch-'0';
		ch=getchar();
	}
	x*=f;
}

template <typename T>
void write(T x)
{
	if (x < 0)
		putchar('-'), x = -x;
	if (x > 9)
		write(x / 10);
	putchar(x % 10 + '0');
}
//=================================================
ll n, mod;
bitset<maxn> vis;
unordered_map<ll, int> dp;
int pri[maxn], tot;
int phi[maxn];
int pre_Sum[maxn];
int inv2, inv6;
inline int ksm(int a, int n)
{
	int res = 1;
	while (n)
	{
		if (n & 1)
			res = 1ll * res * a % mod;
		a = 1ll * a * a % mod;
		n >>= 1;
	}
	return res % mod;
}
inline void init()
{
	inv2 = ksm(2, mod - 2);
	inv6 = ksm(6, mod - 2);
	phi[1] = 1;
	for (int i = 2; i < maxn; i++)
	{
		if (!vis[i])
			pri[tot++] = i, phi[i] = i - 1;
		for (int j = 0; j < tot && i * pri[j] < maxn; j++)
		{
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0)
			{
				phi[i * pri[j]] = (1ll * phi[i]) * pri[j] % mod;
				break;
			}
			phi[i * pri[j]] = (1ll * phi[i]) * (pri[j] - 1) % mod;
		}
	}
	for (int i = 1; i < maxn; i++)
	{
		pre_Sum[i] = (((1ll * (i) * (i)) % mod) * phi[i]) % mod;
	}
	for (int i = 1; i < maxn; i++)
	{
		pre_Sum[i] = (pre_Sum[i] + pre_Sum[i - 1]) % mod;
	}
}

inline int sum1(ll n)
{
	n %= mod;
	return (n * (n + 1) % mod * inv2 % mod) % mod;
}

inline int sum3(ll n)
{
	n %= mod;
	return (n * (n + 1) % mod * (2 * n + 1) % mod * inv6) % mod;
}
int L;
inline int sum2(ll n)
{
	//if(L==1957302059)cerr<<"ok"<<endl;
	//if(L==1957302059)cerr<<(n<maxn)<<endl,cerr<<n<<endl;
	if (n < maxn)
		return pre_Sum[n];
	//if(L==1957302059)cerr<<"ok"<<endl;
	if (dp.count(n))
		return dp[n];
	int tt = sum1(n) % mod;
	int res = 1ll * tt * tt % mod;
	for (ll l = 2, r; l <= n; l = r + 1)
	{
		//if(L==1957302059)cerr<<"**"<<l<<endl;
		ll t = n / l;
		r = n / t;
		int p1 = ((sum3(r) - sum3(l - 1) % mod) + mod) % mod;
		int p2 = sum2(t) % mod;
		res = (1ll * (res - 1ll*p1 * p2 % mod) % mod);
	}
	return dp[n] = (res % mod + mod) % mod;
}

int solve()
{
	int res = 0;
	for (ll l = 1, r; l <= n; l = r + 1)
	{
		ll t = n / l;
		r = n / t;
		//cerr<<l<<" "<<r<<endl;
		int p1 = sum1(t) % mod;
		//if(l==1957302059) cerr<<l<<endl;
		//cerr<<l<<" "<<r<<endl;
		p1 = 1ll * p1 * p1 % mod;
		int p2 = ((sum2(r) - sum2(l - 1)) % mod + mod) % mod;
		//if(l==1957302059) cerr<<l<<endl;
		//cerr<<l<<" "<<r<<endl;
		res = (res % mod + 1ll * p1 * p2 % mod) % mod;
	}
	return (res % mod + mod) % mod;
}

signed main()
{
	//freopen("in.txt", "r", stdin);
	read(mod), read(n);
	init();
	write(solve());
	return 0;
}