P4240 毒瘤之神的考验

题目

P4240 毒瘤之神的考验

神仙题(emmm)

前置

首先有一个很神奇的性质：

(varphi(ij)=dfrac{varphi(i)varphi(j)gcd(i,j)}{varphi(gcd(i,j))})

证明：

[egin{aligned} varphi(i)varphi(j) &= iprodlimits_{p|i}frac{p-1}{p}jprodlimits_{p|j}frac{p-1}{p}\ &= ijprodlimits_{p|ij}frac{p-1}{p}prodlimits_{q|gcd(i,j)}frac{p-1}{p} end{aligned}]

[ herefore varphi(i)varphi(j)gcd(i,j)=varphi(ij)varphi(gcd(i,j)) ]

具体做法

求(sum_{i=1}^nsum_{j=1}^mvarphi(i,j))

我们最终是能得到：$$egin{aligned}
sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}varphi(ij)=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}$$

网上好像很多人的证明都有些错误

证明：

[egin{aligned}\ Ans=&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}]

看我们最终得到的式子，多次查询完全在线做，参数这么多时间肯定承受不住

(sumlimits_{k|T}mu(frac{T}{k})frac{k}{varphi(k)})设有函数(F(x)=sumlimits_{k|x}mu(frac{x}{k})frac{k}{varphi(k)})再普通不过就不多讲了

(sumlimits_{i=1}^{lfloorfrac{n}{T} floor}varphi(iT))参数有两个，设有函数(G(y,x)=sumlimits_{i=1}^{x}varphi(iy))

显然(G(y,x)=G(y,x-1)+varphi(xy))

整个式子设为函数(S(y,z,x)=sumlimits_{T=1}^xsumlimits_{k|T}mu(frac{T}{k})frac{k}{varphi(k)}sumlimits_{i=1}^{y}varphi(iT)sumlimits_{j=1}^{z}varphi(jT))

显然(S(y,z,x)=S(y,z,x-1)+F[x]*G(x,y)*G(x,z))

(F(x)(x<=n),G(y,x)(x<=n,y<=B),S(y,z,x)(x<=n,y,z<=B))

B是自己取的任意参数，或许你会想(frac{n}{T}>B)怎么办，(Rightarrow T<frac{n}{B})，这部分暴力算，后面分块

#include<cstring>
#include<iostream>
#include<cstdio>
#include<string>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long LL;
const int B=35;
const int maxn=1e5+9;
const LL p=998244353;
inline int Read(){
    int x(0),f(1); char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-') f=-1; c=getchar();
    }
    while(c>='0'&&c<='9')
        x=(x<<3)+(x<<1)+c-'0',c=getchar();
    return x*f;
}
bool visit[maxn];
int mu[maxn],phi[maxn],prime[maxn];
int *G[maxn],*S[B+1][B+1],F[maxn];
int inv[maxn];

inline void First(LL N){
    mu[1]=phi[1]=inv[1]=1;
    int tot(0);
    for(int i=2;i<=N;++i){
        if(!visit[i]){
            prime[++tot]=i,
            mu[i]=-1,
            phi[i]=i-1;
        }
        for(int j=1;j<=tot&&i*prime[j]<=N;++j){
            visit[i*prime[j]]=true;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }else{
                phi[i*prime[j]]=phi[i]*phi[prime[j]],
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
    for(int i=2;i<=N;++i)
        inv[i]=p-1ll*p/i*inv[p%i]%p;
    for(int i=1;i<=N;++i)
        for(int j=1;j*i<=N;++j)
            F[i*j]=1ll*(1ll*F[i*j]+1ll*i*inv[phi[i]]%p*mu[j]%p)%p;
    for(int i=1;i<=N;++i){
        G[i]=new int [N/i+1];
        G[i][0]=0;
        for(int j=1;j<=N/i;++j)
            G[i][j]=1ll*(1ll*G[i][j-1]+1ll*phi[j*i])%p;
    }
    for(int j=1;j<=B;++j)
        for(int k=1;k<=B;++k){
        	int len(N/(max(j,k)));
        	S[j][k]=new int [len+1];
        	S[j][k][0]=0;
        	for(int i=1;i<=len;++i)
        		S[j][k][i]=1ll*(1ll*S[j][k][i-1]+1ll*F[i]*G[i][j]%p*G[i][k]%p)%p;
        }
}
inline LL Solve(int n,int m){
    if(n>m)
        swap(n,m);
    LL ans(0);
    for(int i=1;i<=m/B;++i)
        ans=(ans+1ll*F[i]*G[i][n/i]%p*G[i][m/i]%p)%p;
    for(int l=m/B+1,r;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=(ans+1ll*(S[n/l][m/l][r]-S[n/l][m/l][l-1]+p)%p)%p;
    }
    return ans;
}
int main(){
    int T=Read();
    First(100000);
    while(T--){
        int n(Read()),m(Read());
        printf("%lld
",Solve(n,m));
    }
    return 0;
}/*
3
1 1
2 2
3 3

1 5 19
*/