【loj6179】Pyh的求和
Solution
这题其实有一个式子一喵一样的版本在bzoj，但是那题是\(m\)特别大然后只有一组数据
这题多组数据==
首先根据\(\varphi(x)\)的通项\(\varphi(x)=x\prod\limits_{i=1}^{n}(1-\frac{1}{p1_i})=\prod\limits_{i=1}^{m}(p_i-1)p_i^{a_i-1}\)（其中\(n\)是\(x\)分解质因数之后没有去重的质因数列表\(p1\)的长度，\(m\)是去重之后质因数列表\(p\)的长度，\(x=\prod\limits_{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\)：
稍微说一下最后一步是相当于枚举\(dk\)，也就是令\(T=dk\)然后枚举\(T\)
然后我们可以令\(g(x)=\sum\limits_{k|x}\mu(\frac{x}{k})\frac{k}{\varphi(k)}\)，令\(s(i,j)=\sum\limits_{k=1}^j\varphi(ik)\)
那么这个式子就可以写成：
接下来看起来是怎么化也化不动了qwq但是我们发现\(g(i)\)和\(s(i,j)\)都可以在调和级数的复杂度内预处理出来，但是再接下来我们发现\(O(n)\)求解显然是不现实的
这时候当然是要大力分段啊，只不过光是普通的操作还是不行（前缀和这个东西很难搞），这里我们还需要一个黑科技，我们手动设定一个阈值\(TOP\)，然后对于\(i<=\frac{m}{TOP}\)的情况我们暴力算，对于\(i>\frac{m}{TOP}\)的情况，我们再预处理一个\(T[i][j][k]\)（也就是前缀和）：
然后当\(i>\frac{m}{TOP}\)的时候\(\lfloor\frac{m}{i}\rfloor<=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\n",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();
}
}
