P13645 Totient with Divisors

Luogu 链接

题意

求

\[\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\sigma(ij) \]

多测，\(1\le n,m,T\le N=10^5\)。

思路

首先有一个不知道哪里冒出来的结论：

\[\sigma(ij)=\sum_{d_1\mid i}\sum_{d_2\mid j}\dfrac{id_2}{d_1}[\gcd(d_1,d_2)=1] \]

先给出这个结论的证明：

对于质数 \(p\)，我们设它在 \(i\) 中的指数为 \(x\)，在 \(j\) 中的指数为 \(y\)，在 \(d_1\) 中的指数为 \(x'\)，在 \(d_2\) 中的指数为 \(y'\)，所以 \(p\) 在 \(ij\) 中的指数为 \(x+y\)。（显然有 \(0\le x'\le x\)，\(0\le y'\le y\)）

则 \(p\) 在 \(\dfrac{id_2}{d_1}\) 的指数为 \(x+y'-x'\)，并且由于 \(\gcd(d_1,d_2)=1\)，所以有 \(\min(x',y')=0\)。

• 当 \(x'=0\) 时，有 \(0\le y'\le y\)，因此 \(x+y'-x'\) 可以取遍 \([x,x+y]\) 中的所有整数。
• 当 \(x'\neq0\) 时，有 \(y'=0\)，以及 \(0<x'\le x\)，所以 \(x+y'-x'\) 可以取遍 \([0,x)\) 中的所有整数。

综上 \(x+y'-x'\) 可以取遍 \([0,x+y]\) 中的所有整数，故原式成立。

然后我们就可以来推式子了（不妨设 \(n\le m\)）。

\[\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\sigma(ij)\\ =&\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\sum_{d_1\mid i}\sum_{d_2\mid j}\dfrac{id_2}{d_1}[\gcd(d_1,d_2)=1]\\ =&\sum_{i=1}^n\varphi(i)i\sum_{j=1}^m\varphi(j)\sum_{d_1\mid i}\sum_{d_2\mid j}\dfrac{d_2}{d_1}[\gcd(d_1,d_2)=1]\\ =&\sum_{d_1=1}^n(\sum_{i=1}^{ n/d_1}\varphi(id_1)i)\sum_{d_2=1}^md_2(\sum_{j=1}^{m/d_2}\varphi(jd_2))[\gcd(d_1,d_2)=1]\\ =&\sum_{t=1}^n\mu(t)t(\sum_{d_1=1}^{n/t}\sum_{i=1}^{n/d_1t}\varphi(id_1t)i)(\sum_{d_2=1}^{m/t}d_2\sum_{i=1}^{m/d_2t}\varphi(id_2t))\\ =&\sum_{t=1}^n\mu(t)t(\sum_{d=1}^{n/t}\sum_{i=1}^{n/dt}\varphi(idt)i)(\sum_{d=1}^{m/t}\sum_{i=1}^{m/dt}\varphi(idt)i)\\ =&\sum_{t=1}^n\mu(t)tf(\lfloor\dfrac{n}{t}\rfloor,t)f(\lfloor\dfrac{m}{t}\rfloor,t) \end{aligned}\]

---

我们有：

\[\begin{aligned} f(x,y) &=\sum_{i=1}^x\sum_{j=1}^{x/i}\varphi(ijy)j\\ &=\sum_{t=1}^x\varphi(ty)\sigma(t)\\ &=(\sum_{t=1}^{x-1}\varphi(ty)\sigma(t))+\varphi(xy)+\sigma(x)\\ &=f(x-1,y)+\varphi(xy)+\sigma(x) \end{aligned}\]

显然有 \(f(1,y)=\varphi(y)\)

注意到，\(\lfloor\dfrac{n}{t}\rfloor\times t\le n\le N\)，\(\lfloor\dfrac{m}{t}\rfloor\times t\le m\le N\)。
因此我们需要用到的 \(f(x,y)\) 都有 \(1\le xy\le N\)。

所以我们需要计算的 \(f(x,y)\) 的个数为 \(\displaystyle O(\sum_{i=1}^N\dfrac{N}{i})=O(N\log N)\)

并且，在我们线性筛出 \([1,N]\) 处 \(\varphi(x)\)、\(\sigma(x)\) 的值后，我们可以 \(O(1)\) 递推算出每个 \(f(x,y)\)。

因此预处理所有 \(f(x,y)\) 的复杂度是 \(O(N\log N)\) 的。

---

注意到我们显然需要用数论分块计算答案。
也就是说，我们需要处理 \(\mu(t)tf(\lfloor\dfrac{n}{t}\rfloor,t)f(\lfloor\dfrac{m}{t}\rfloor,t)\) 在 \(\lfloor\dfrac{n}{t}\rfloor\) 和 \(\lfloor\dfrac{m}{t}\rfloor\) 相同时关于 \(t\) 的前缀和。

但是总共有 \(\displaystyle O(\sum_{t=1}^N\frac{N^2}{t^2})=O(N^2)\) 个需要计算的值，全部预处理复杂度无法接受。

但是我们肯定也要预处理一部分，不然计算的时间复杂度就无法接受。

考虑设置一个阈值 \(S\)，当 \(\lfloor\dfrac{m}{t}\rfloor\le S\) 时将值预处理。
复杂度是 \(\displaystyle O(\sum_{i=1}^Si\times\dfrac{N}{i})=O(NS)\)。

当 \(\lfloor\dfrac{m}{t}\rfloor>S\) 时，此时由于 \(t\le \lfloor\dfrac{m}{S+1}\rfloor\)，我们可以对于每个 \(t\) 处的值 \(O(1)\) 计算。

当 \(t>\lfloor\dfrac{m}{S+1}\rfloor\) 时，我们就可以用预处理的值进行数论分块了。

时间复杂度是 \(O(NS+T\sqrt{N}+\dfrac{NT}{S})\)。
当 \(S\) 取 \(\sqrt{T}\) 时可以达到最优时间复杂度 \(O(N\sqrt{T}+T\sqrt{N})\)。

对于具体的实现，我们可以令 \(\displaystyle g(x,y,t)=\sum_{i=1}^t\mu(i)if(x,i)f(y,i)\)，直接递推计算即可。
显然递推是 \(O(1)\) 的，因为 \(f(x,y)\) 的值我们都计算过了。

---

现在我们来理一下思路。

预处理：

• 使用线性筛筛出 \(\mu(n)\)、\(\varphi(n)\)、\(\sigma(n)\) 在 \([1,N]\) 的值，复杂度 \(O(N)\)。
• \(O(1)\) 递推出每个 \(f(x,y)\) 的值，复杂度 \(O(N\log N)\)。
• \(O(1)\) 递推出我们需要预处理的 \(g(x,y,t)\) 的值，复杂度 \(O(N\sqrt{T})\)。

计算答案：

• 对于 \(t\le \lfloor\dfrac{m}{S+1}\rfloor\)，我们 \(O(1)\) 计算每个值，复杂度 \(O(N\sqrt{T})\)。
• 对于 \(t> \lfloor\dfrac{m}{S+1}\rfloor\)，我们可以利用预处理的 \(g(x,y,t)\) 使用数论分块计算，复杂度 \(O(T\sqrt{N})\)。

程序

AC 记录

#include<bits/stdc++.h>
#define forUp(i,a,b) for(int i=(a);i<=(b);++i)
#define forUP(i,a,b) for(int i=(a);i<(b);++i)
#define forDown(i,a,b) for(int i=(a);i>=(b);--i)
#define forG(i,u,v) for(int i=head[u],v=to[i];i;i=nxt[i],v=to[i])
#define pb emplace_back
using ll=long long;using ull=unsigned long long;using uint=unsigned int;using db=double;using ld=long db;using pii=std::pair<int,int>;using pdi=std::pair<db,int>;using vl=__int128;using uvl=unsigned __int128;
constexpr int INF=0x3f3f3f3f,MINF=0xcfcfcfcf;constexpr long long INFLL=0x3f3f3f3f3f3f3f3f,MINFLL=0xcfcfcfcfcfcfcfcf;constexpr double INFDB=1e50,eps=1e-9;
template<class _Tp>void chkMax(_Tp &x,const _Tp &y){x<y?x=y:0;}template<class _Tp>void chkMin(_Tp &x,const _Tp &y){x>y?x=y:0;}
constexpr int N=1e5+10,B=sqrt(1e5)+10,mod=998244353;int __test_num=1,__test_id;using namespace std;void __init();

int cnt,P[N],mu[N],phi[N],minp[N],sigma[N];bool np[N];
int S;vector<int> f[N],g[B][B];

int n,m,ans;

void __solve(int __test_id){
	scanf("%d%d",&n,&m);if(n>m)swap(n,m);ans=0;
	for(int i=1;i<=n&&i<=m/(S+1);++i)ans=(ans+(1ll*mu[i]*i*f[n/i][i]%mod*f[m/i][i]%mod+mod)%mod)%mod;
	for(int l=m/(S+1)+1,r;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans=(ans+(g[n/l][m/l][r]-g[n/l][m/l][l-1]+mod)%mod)%mod;
	}
	printf("%d\n",ans);
}
signed main(){
	__init();
	forUp(i,1,__test_num)__solve(i);
	return 0;
}
void __init(){
	//const string __file_name="test";freopen((__file_name+".in").c_str(),"r",stdin);freopen((__file_name+".out").c_str(),"w",stdout);
	scanf("%d",&__test_num);
	int n=N-10;S=sqrt(__test_num);
	mu[1]=phi[1]=sigma[1]=minp[1]=1;
	forUp(i,2,n){
		if(!np[i])P[++cnt]=i,mu[i]=-1,phi[i]=i-1,sigma[i]=i+1,minp[i]=i;
		for(int j=1;j<=cnt&&i*P[j]<=n;++j){
			int p=P[j],x=i*p;
			np[x]=true;
			if(i%p==0){
				minp[x]=minp[i]*p;
				phi[x]=phi[i]*p;
				if(minp[x]==x)sigma[x]=p*sigma[i]+1;
				else sigma[x]=sigma[x/minp[x]]*sigma[minp[x]];
				break;
			}
			mu[x]=-mu[i],phi[x]=phi[i]*phi[p],sigma[x]=sigma[i]*sigma[p],minp[x]=p;
		}
	}
	f[1].pb(0);forUp(j,1,n)f[1].pb(phi[j]);
	forUp(i,2,n){
		f[i].pb(0);
		forUp(j,1,n/i){
			int tmp=(f[i-1][j]+1ll*phi[i*j]*sigma[i]%mod)%mod;
			f[i].pb(tmp);
		}
	}
	forUp(j,1,S)forUp(i,1,j){
		g[i][j].pb(0);
		forUp(k,1,n/j){
			int tmp=(g[i][j][k-1]+(1ll*mu[k]*k*f[i][k]%mod*f[j][k]%mod+mod)%mod)%mod;
			g[i][j].pb(tmp);
		}
	}
}