【loj6179】Pyh的求和

Portal -->loj6179

Solution

​　　 这题其实有一个式子一喵一样的版本在bzoj，但是那题是(m)特别大然后只有一组数据

　　 这题多组数据==

​　　

　　 首先根据(varphi(x))的通项(varphi(x)=xprodlimits_{i=1}^{n}(1-frac{1}{p1_i})=prodlimits_{i=1}^{m}(p_i-1)p_i^{a_i-1})（其中(n)是(x)分解质因数之后没有去重的质因数列表(p1)的长度，(m)是去重之后质因数列表(p)的长度，(x=prodlimits_{i=1}^{m} p_i^{a_i})）我们有(varphi(i*j)=varphi(i)*varphi(j)*frac{gcd(i,j)}{varphi(gcd(i,j))})，具体就是因为(varphi(i)*varphi(j))中(gcd)的质因子的部分被算了两次，但是除掉(varphi(gcd(i,j)))之后又没有将(gcd)对(a_i)的贡献算上

​　　 然后我们就可以快乐推式子了，为了让接下来的式子更加简洁，我们默认(n<=m)：

[egin{aligned} &sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}varphi(ij)\ =&sumlimits_{i=1}^n sumlimits_{j=1}^m varphi(i)varphi(j)frac{gcd(i,j)}{varphi(gcd(i,j))}\ =&sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^mvarphi(i)varphi(j)frac{d}{varphi(d)}[gcd(i,j)=d]\ =&sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}varphi(id)varphi(jd)frac{d}{varphi(d)}[gcd(i,j)=1]\ =&sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}varphi(id)varphi(jd)frac{d}{varphi(d)}sumlimits_{k|i,k|j}mu(k)\ =&sumlimits_{d=1}^nfrac{d}{varphi(d)}sumlimits_{k=1}^{lfloorfrac{n}{d} floor}mu(k)sumlimits_{i=1}^{lfloorfrac{n}{dk} floor}varphi(idk)sumlimits_{j=1}^{lfloorfrac{m}{dk} floor}varphi(jdk)\ =&sumlimits_{T=1}^nsumlimits_{k|T}mu(frac{T}{k})frac{k}{varphi(k)}sumlimits_{i=1}^{lfloorfrac{n}{T} floor}varphi(iT)sumlimits_{j=1}^{lfloorfrac{m}{T} floor}varphi(jT)\ end{aligned} ]

​　　 稍微说一下最后一步是相当于枚举(dk)，也就是令(T=dk)然后枚举(T)

​　　 然后我们可以令(g(x)=sumlimits_{k|x}mu(frac{x}{k})frac{k}{varphi(k)})，令(s(i,j)=sumlimits_{k=1}^jvarphi(ik))

　　 那么这个式子就可以写成：

[sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}varphi(ij)=sumlimits_{i=1}^ng(i)cdot s(i,lfloorfrac{n}{i} floor)cdot s(i,lfloorfrac{m}{i} floor) ]

​
　　 
​　　 接下来看起来是怎么化也化不动了qwq但是我们发现(g(i))和(s(i,j))都可以在调和级数的复杂度内预处理出来，但是再接下来我们发现(O(n))求解显然是不现实的

​　　 这时候当然是要大力分段啊，只不过光是普通的操作还是不行（前缀和这个东西很难搞），这里我们还需要一个黑科技，我们手动设定一个阈值(TOP)，然后对于(i<=frac{m}{TOP})的情况我们暴力算，对于(i>frac{m}{TOP})的情况，我们再预处理一个(T[i][j][k])（也就是前缀和）：

[egin{aligned} T[i][j][k]&=sumlimits_{p=1}^kg(p)sumlimits_{t=1}^ivarphi(tp)sumlimits_{t=1}^ivarphi(tp)\ &=sumlimits_{p=1}^kg(p)cdot s(p,i)cdot s(p,j) end{aligned} ]

​　　 然后当(i>frac{m}{TOP})的时候(lfloorfrac{m}{i} floor<=TOP)，所以我们(i,j)只要预处理到(TOP)然后直接用普通的分段操作来搞就好了

　　

　　 代码大概长这个样子

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
using namespace std;
const int N=1e5+10,MOD=998244353,TOP=35,inf=2147483647;
int p[N],g[N],miu[N],phi[N],inv[N];
vector<int>S[N];//s[i][j]=sum_{k=1}^{j}phi[i*k]
vector<int> T[TOP+1][TOP+1];//T[i][j][k]=sum_{p=1}^{k}sum_{t=1}^{i}phi(t*p)sum_{t=1}^{j}phi(t*p)
int vis[N];
int n,m,ans,T1;
void prework(int n){
	int cnt=0;
	miu[1]=1; phi[1]=1;
	for (int i=2;i<=n;++i){
		if (!vis[i])
			miu[i]=-1,p[++cnt]=i,phi[i]=i-1;
		for (int j=1;j<=cnt&&i*p[j]<=n;++j){
			vis[i*p[j]]=true;
			if (i%p[j]==0){
				miu[i*p[j]]=0; phi[i*p[j]]=phi[i]*p[j];
				break;
			}
			else
				miu[i*p[j]]=-miu[i],phi[i*p[j]]=phi[i]*phi[p[j]];
		}
	}
	inv[1]=1;
	for (int i=2;i<=n;++i) inv[i]=1LL*(MOD-MOD/i)*inv[MOD%i]%MOD;
	for (int i=1;i<=n;++i)
		for (int j=i;j<=n;j+=i){
			if (j==9)
				int debug=1;
			g[j]=(1LL*g[j]+1LL*miu[j/i]*(1LL*(i)*inv[phi[i]]%MOD)+MOD)%MOD;
		}
	for (int i=1;i<=n;++i){
		S[i].push_back(0);
		for (int j=1;j<=n/i;++j) S[i].push_back((S[i][j-1]+phi[i*j])%MOD);
	}
	for (int i=1;i<=TOP;++i)
		for (int j=1;j<=TOP;++j){
			T[i][j].push_back(0);
			for (int k=1;k<=n/i&&k<=n/j;++k)
				T[i][j].push_back((1LL*T[i][j][k-1]+1LL*g[k]*S[k][i]%MOD*S[k][j]%MOD)%MOD);
		}
	
}
void solve(){
	if (n>m) swap(n,m);
	ans=0;
	for (int i=1;i<=m/TOP;++i)
		ans=(1LL*ans+1LL*g[i]*S[i][n/i]%MOD*S[i][m/i]%MOD)%MOD;
	for (int i=m/TOP+1,pos=0;i<=n;i=pos+1){
		pos=min(m/(m/i),(n/i)?n/(n/i):inf);
		ans=(1LL*ans+(1LL*T[n/i][m/i][pos]+MOD-T[n/i][m/i][i-1])%MOD)%MOD;
	}
	printf("%lld
",ans);
}

int main(){
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
#endif
	prework(N-10);
	//prework(10);
	scanf("%d",&T1);
	for (int o=1;o<=T1;++o){
		scanf("%d%d",&n,&m);
		solve();
	}
}