【分析】
规定三个求和项：
(displaystyle Sum_1(n)=sum_{i=1}^ni={n(n+1)over 2})
(displaystyle Sum_2(n)=sum_{i=1}^ni^2={n(n+1)(2n+1)over 6})
(displaystyle Sum_3(n)=sum_{i=1}^ni^3=({n(n+1)over 2})^2=Sum^2_1(n))
关于操作 (4) ，是否进行高斯消元并无影响，故放置后文进行分析
高斯消元之前，对于操作 (2) 有：
(displaystylequadsum_{i=1}^xsum_{j=1}^xA_{ij})
(displaystyle=sum_{i=1}^xsum_{j=1}^xijgcd(i,j))
(displaystyle=sum_{i=1}^xsum_{j=1}^xijsum_{dmid iwedge dmid j}oldsymbol varphi(d))
(displaystyle=sum_{d=1}^x(oldsymbol{id}^2cdot oldsymbol varphi)(d)sum_{i=1}^{x/d}isum_{j=1}^{x/d}j)
(displaystyle=sum_{d=1}^x(oldsymbol{id}^2cdot oldsymbol varphi)(d)cdot Sum_3(x/d))
前半部分会达到 (10^8) ，因此需要使用杜教筛：
设 (oldsymbol f=oldsymbol {id}^2cdot oldsymbol varphi,oldsymbol g=oldsymbol {id}^2,oldsymbol h=oldsymbol f*oldsymbol g=oldsymbol {id}^3)
(displaystyle H(n)=sum_{i=1}^noldsymbol h(i)=sum_{i=1}^nsum_{dmid i}oldsymbol f(d)oldsymbol g({iover d})=sum_{d=1}^noldsymbol g(d)sum_{i=1}^{n/d}oldsymbol f(i)=sum_{d=1}^noldsymbol g(d)F(n/d)=sum_{d=2}^noldsymbol g(d)F(n/d)+F(n))
(displaystyle F(n)=H(n)-sum_{d=2}^noldsymbol g(d)F(n/d),H(n)=Sum_3(n),G(n)=Sum_2(n))
现考虑高斯消元后的形式：
先给出结论：({ijgcd(i,j)} o{ijoldsymbol varphi(i)cdot [imid j]})
证明放在后文，现考虑其他部分：
(2) 操作直接分类讨论：(egin{cases} ijoldsymbol varphi(i),imid j \ 0,i mid j end{cases})
(4) 操作由线性代数的知识得到：(displaystyle det B=prod_{i=1}^x[icdot icdot oldsymbol varphi(i)]=prod_{i=1}^xoldsymbol f(i)=Pf(x))
(3) 操作我们可以很容易列出式子：(displaystyle sum_{i=1}^xsum_{j=1}^xijoldsymbol varphi(i)cdot [imid j]=sum_{i=1}^xsum_{imid j}ijoldsymbol varphi(i))
调换顺序，枚举 (j:displaystyle sum_{j=1}^xjsum_{imid j}ioldsymbol varphi(i)=sum_{j=1}^x(oldsymbol {id}cdot oldsymbol varphi*oldsymbol Icdot oldsymbol {id})(j))
所以，记 (oldsymbol g=oldsymbol {id}cdot oldsymbol varphi*oldsymbol Icdot oldsymbol {id}) 则答案为 (G(x))
由于 (oldsymbol g) 是积性函数，考虑如何欧拉筛：
(displaystyle (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(p^k)=sum_{i=0}^k(oldsymbol {id}cdot oldsymbol varphi)(p^i)=sum_{i=1}^kp^icdot p^{i-1}cdot (p-1)+1={pover p+1}cdot (p^{2k}-1)+1)
欧拉筛的时候，除了维护数字 (n) 的最小质因数 (fc_n) ，同时维护最小值因数的 出现次数 次方 (pk_n)
(displaystyle (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(n imes p)=egin{cases} (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(n)cdot (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(p),p mid n \ (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(n)cdot ({fc_nover fc_n+1}cdot (pk_n^2-1)+1)^{-1}cdot ({fc_{n imes p}over fc_{n imes p}+1}cdot (pk_{n imes p}^2-1)+1),pmid n end{cases})
最后 (displaystyle oldsymbol g(n)=displaystyle (oldsymbol {id}cdot oldsymbol varphi*oldsymbol I)(n)cdot n)
【代码】
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=5e6+10,MOD=998244353;
int Fc[MAXN],Prime[MAXN],CntPrime,Phi[MAXN],F[MAXN],PF[MAXN],G[MAXN],PK[MAXN];
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD; return ans; }
inline ll inv(ll a) { return fpow(a,MOD-2); }
inline int add(int a,int b) { return (a+b>=MOD)?(a+b-MOD):(a+b); }
inline int dis(int a,int b) { return (a-b<0)?(a-b+MOD):(a-b); }
inline void sieve(){
Phi[1]=G[1]=1;
for(int i=2;i<=5e6;++i){
if(!Fc[i]){
Fc[i]=Prime[++CntPrime]=i;
Phi[i]=i-1;
PK[i]=i;
G[i]=((ll)i*i-i+1+MOD)%MOD;
}
for(int j=1;j<=CntPrime;++j)
if(Prime[j]*i>5e6||Prime[j]>Fc[i]) break;
else if(Prime[j]==Fc[i]) {
PK[Prime[j]*i]=PK[i]*Prime[j];
Fc[Prime[j]*i]=Prime[j];
Phi[Prime[j]*i]=Phi[i]*Prime[j];
ll tmp=Prime[j]*inv(Prime[j]+1)%MOD;
ll tmp0=tmp*((ll)PK[i]*PK[i]%MOD-1)%MOD+1;
ll tmp1=tmp*((ll)PK[Prime[j]*i]*PK[Prime[j]*i]%MOD-1)%MOD+1;
G[Prime[j]*i]=G[i]*inv(tmp0)%MOD*tmp1%MOD;
}
else{
PK[Prime[j]*i]=Prime[j];
Fc[Prime[j]*i]=Prime[j];
Phi[Prime[j]*i]=Phi[i]*(Prime[j]-1);
G[Prime[j]*i]=(ll)G[Prime[j]]*G[i]%MOD;
}
}
PF[0]=1;
for(int i=1;i<=5e6;++i){
G[i]=((ll)G[i]*i+G[i-1])%MOD;
F[i]=(ll)Phi[i]*i%MOD*i%MOD;
PF[i]=(ll)PF[i-1]*F[i]%MOD;
F[i]=add(F[i],F[i-1]);
}
}
bool app;
unordered_map<ll,int> sumf;
const int inv2=(MOD+1>>1),inv6=(MOD+1)/6;
inline int sum1(int n) { return (ll)inv2*n%MOD*(n+1)%MOD; }
inline int sum2(int n) { return (ll)inv6*n%MOD*(n+1)%MOD*(n<<1|1)%MOD; }
inline int sum3(int n) { int res=(ll)inv2*n%MOD*(n+1)%MOD; return (ll)res*res%MOD; }
inline int getf(int n){
if(n<=5e6) return F[n];
if(sumf.count(n)) return sumf[n];
int res=sum3(n);
for(register int l=2,r;l<=n;l=r+1){
r=n/(n/l);
res=dis( res , (sum2(r)-sum2(l-1)+(ll)MOD)*getf(n/l)%MOD );
}
return sumf[n]=res;
}
inline int ans(int x){
if(app) return G[x];
int res=0;
for(register int l=1,r;l<=x;l=r+1){
r=x/(x/l);
res=add( res , sum3(x/l)*(getf(r)-getf(l-1)+(ll)MOD)%MOD );
}
return res;
}
int main(){
sieve();
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
ll M,opt,x,y; cin>>M;
while(M--&&cin>>opt){
if(opt!=1) cin>>x;
if(opt==2) cin>>y;
if(opt==1) app=1;
else if(opt==2) cout<<( app? (y%x?0:x%MOD*(y%MOD)%MOD*Phi[x]%MOD) : (x%MOD*(y%MOD)%MOD*(__gcd(x,y)%MOD)%MOD) )<<"
";
else if(opt==3) cout<<ans(x)<<"
";
else if(opt==4) cout<<PF[x]<<"
";
}
cout.flush();
return 0;
}
【证明】
证明分为两步：
-
证明矩阵 (A,B) 若满足 (b_{ij}=ija_{ij}) ，则分别对其高斯消元后得到 (A',B') ，仍满足 (b'_{ij}=ija'_{ij})
-
证明矩阵 ({gcd(i,j)}) 高斯消元后得到 ({oldsymbol varphi(i)cdot [imid j]})
关于 (1) ，记 (Lambda = ext{diag}(1,2,3,cdots )) 则 (Lambda^{-1}= ext{diag}(1,2^{-1},3^{-1},cdots ))
则 (B=Lambda ALambda)
设高斯消元的矩阵为 (G) ，即 (A'=GA,B'=GB)
不难发现，(G) 为一系列初等矩阵相乘得出的，记为 (G=prod_i E_i)
因为 ((Lambdacdot E_i)^{-1}=Lambda^{-1}cdot E_i^{-1})
故 (Lambda) 与各个 (E_i) 都是可交换的
(displaystyle herefore Lambda G=Lambda prod_iE_i=E_1Lambdaprod_{i>1}E_i=cdots =prod_i E_iLambda=GLambda)
因此得出：(GB=GLambda ALambda=Lambda GALambda)
这表明，先将 (A) 去除 (ij) ，再高斯消元，最后再乘回去；和直接高斯消元是等价的
由此证明结论 (1)
现假设前 ((i-1)) 行高斯消元后满足该结论
则第 (i) 行减去 (d) 行 ((dmid i)) 得到：
(displaystyle a'_{ij}=gcd(i,j)-sum_{dmid gcd(i,j)wedge d<i}oldsymbol varphi(i)=gcd(i,j)-sum_{dmid gcd(i,j)}oldsymbol varphi(i)+[gcd(i,j)=i]cdot oldsymbol varphi(i)=gcd(i,j)-gcd(i,j)+[imid j]cdot oldsymbol varphi(i)=oldsymbol varphi(i)cdot [imid j])
则 (j<i) 时 (i mid j o a'_{ij}=0)
并且 (j=i) 时 (imid j o a'_{ij}=oldsymbol varphi(i))
成功消出 (i) 行的上三角矩阵，且仍然满足结论
由数学归纳法，显然可证：矩阵 ({gcd(i,j)}) 高斯消元后得到 ({oldsymbol varphi(i)cdot [imid j]})
所以先取出 (ij) ，高斯消元得到 (a'_{ij}=oldsymbol varphi(i)cdot [imid j]) ，再算回去得到 (b'_{ij}=ijoldsymbol varphi(i)cdot [imid j])
因此，原结论得证