先扯些别的。
2021 年 7 月的某一天,我和 ycx 对话:
- tzc:你做过哪些名字里带“毒瘤”的题目,我做过一道名副其实的毒瘤题就叫毒瘤,是个虚树+dp
- ycx:还有毒瘤之神的考验
- tzc:???那是个啥?
- ycx:一道数论水题
然后我便做到了这个题,然后却发现它一点也不水……
跑题了跑题了
首先我们显然不可能硬着头皮算 (varphi(ij)),肯定要想办法将 (varphi(ij)) 中的 (i,j) 独立开来。通过这题的套路可知 (varphi(ij)=dfrac{varphi(i)varphi(j)gcd(i,j)}{varphi(gcd(i,j))}),于是下面我们就可以开始推式子了:
如果我们记 (f(n)=sumlimits_{dmid n}dfrac{d}{varphi(d)}·mu(dfrac{n}{d})),(S(n,m)=sumlimits_{i=1}^mvarphi(ni)),那么
(f) 显然可以调和级数地求出,有用的 (S(n,m)) 的个数也是 (nln n) 级别的,因此暴力计算上式可实现 (Tn+nln n) 的复杂度,无法通过。直接整除分块看起来不太容易的亚子,这就使我们陷入了一个比较尴尬的局面。但是一个非常直观的 observation 是当 (T) 比较大时,可能的 ((lfloordfrac{n}{T} floor,lfloordfrac{m}{T} floor)) 组成的二元组并不是特别多。因此考虑从这个角度入手。记 (B=316),对于 (Tle B),我们暴力计算即可,对于 (T>B),写个程序算一下可以发现 (sumlimits_{i=317}^{10^5}lfloordfrac{10^5}{i} floor^2) 只有大概 (3 imes 10^7) 的样子,因此这部分我们就预处理时,暴力枚举所有可能的 ((lfloordfrac{n}{T} floor,lfloordfrac{m}{T} floor)),然后前缀和算一下贡献,这样每次询问整除分块即可求出这部分 (T) 的贡献。
时间复杂度 (Tsqrt{n}+nsqrt{n}+nln n),可以通过。
const int MAXN=1e5;
const int B=316;
int inv[MAXN+5],pr[MAXN/7+5],prcnt=0,phi[MAXN+5],vis[MAXN+5],f[MAXN+5],mu[MAXN+5];
vector<int> s[MAXN+5],ss[B+4][B+4];
void init(){
for(int i=(inv[0]=inv[1]=1)+1;i<=MAXN;i++) inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
phi[1]=mu[1]=1;
for(int i=2;i<=MAXN;i++){
if(!vis[i]) phi[i]=i-1,mu[i]=-1,pr[++prcnt]=i;
for(int j=1;j<=prcnt&&pr[j]*i<=MAXN;j++){
vis[pr[j]*i]=1;
if(i%pr[j]==0){phi[i*pr[j]]=phi[i]*pr[j];break;}
phi[i*pr[j]]=phi[i]*phi[pr[j]];mu[i*pr[j]]=-mu[i];
}
}
for(int i=1;i<=MAXN;i++) for(int j=i;j<=MAXN;j+=i)
f[j]=(0ll+f[j]+1ll*i*inv[phi[i]]%MOD*mu[j/i]%MOD+MOD)%MOD;
for(int i=1;i<=MAXN;i++){
s[i].resize(MAXN/i+1);
for(int j=1;j<=MAXN/i;j++) s[i][j]=(s[i][j-1]+phi[i*j])%MOD;
}
for(int i=1;i<=B+1;i++) for(int j=1;j<=B+1;j++){
ss[i][j].pb(0);int sum=0;
for(int k=B+1;k<=MAXN;k++){
if(k*i>MAXN||k*j>MAXN) break;
sum=(sum+1ll*f[k]*s[k][i]%MOD*s[k][j])%MOD;
ss[i][j].pb(sum);
}
}
}
int main(){
init();int qu;scanf("%d",&qu);
while(qu--){
int n,m;scanf("%d%d",&n,&m);int res=0;
for(int i=1;i<=B;i++) res=(res+1ll*f[i]*s[i][n/i]%MOD*s[i][m/i])%MOD;
for(int l=B+1,r;l<=min(n,m);l=r+1){
r=INF;
if(n/l) chkmin(r,n/(n/l));
if(m/l) chkmin(r,m/(m/l));
res=(0ll+res+ss[n/l][m/l][r-B]-ss[n/l][m/l][l-B-1]+MOD)%MOD;
}
printf("%d
",res);
}
return 0;
}