题意
给定 (n,c,d,q)，(q) 次给定 (b_1,dots b_n)，解方程组：
[sum_{j = 1}^{n} gcd(i, j)^c cdot operatorname{lcm}(i, j)^d cdot x_j equiv b_i pmod{p}
]
(nleq 3 imes 10^5)，(qleq 3)。
题解
(operatorname{lcm}(x,y)=dfrac{xy}{gcd(x,y)})
原式相当于是
[sum_{j = 1}^{n} f(gcd(i,j))g(i)h(j) x_iequiv b_i pmod{p}
]
稍微化一下：
[sum_{j = 1}^{n} f(gcd(i,j))x'_jequiv b'_i pmod{p}
]
[sum_{j = 1}^{n} sum_{d|i,d|j} (f*mu)(d)x'_jequiv b'_i pmod{p}
]
[sum_{d|i} (f*mu)(d) sum_{j=1}^{n} [d|j]x'_jequiv b'_i pmod{p}
]
令 (z(d)=sum_{j=1}^{n} [d|j]x'_j)。
[sum_{d|i} (f*mucdot z)(d)equiv b'_i pmod{p}
]
[((f*mucdot z)*1)(i)equiv b'_i pmod{p}
]
[(f*mucdot z)equiv b'*mu pmod{p}
]
[zequiv dfrac{b'*mu}{f*mu} pmod{p}
]
求出 (z) 后再莫反得 (x_i)，一共 3 个莫反。(b'*mu) 和 (f*mu) 点除时可能有无解或多解。
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int qpow(int x,int y){
int ans=1;
while(y){
if(y&1)ans=ans*1ll*x%mod;
x=x*1ll*x%mod;
y>>=1;
}
return ans;
}
const int N=1e5+10;
int b[N],x[N],n,c,d,q;
int ipw[N];
int f[N],i_f[N],fz[N],z[N];
int main(){
scanf("%d%d%d%d",&n,&c,&d,&q);
c%=mod-1;d%=mod-1;
for(int i=1;i<=n;i++)ipw[i]=qpow(i,mod-1-d);
for(int i=1;i<=n;i++)f[i]=qpow(i,mod-1+c-d);
for(int i=1;i<=n;i++)
for(int j=i+i;j<=n;j+=i)
f[j]=(f[j]+mod-f[i])%mod;
for(int i=1;i<=n;i++)i_f[i]=qpow(f[i],mod-2);
while(q--){
for(int i=1;i<=n;i++)scanf("%d",b+i);
for(int i=1;i<=n;i++)fz[i]=b[i]*1ll*ipw[i]%mod;
for(int i=1;i<=n;i++)
for(int j=i+i;j<=n;j+=i)
fz[j]=(fz[j]+mod-fz[i])%mod;
bool ok=1;
for(int i=1;i<=n;i++){
if(!f[i]){
if(fz[i]){
ok=0;
break;
}else x[i]=0;
}else x[i]=fz[i]*1ll*i_f[i]%mod;
}
if(!ok){
puts("-1");
continue;
}
for(int i=n;i>=1;i--)
for(int j=i+i;j<=n;j+=i)
x[i]=(x[i]+mod-x[j])%mod;
for(int i=1;i<=n;i++)x[i]=x[i]*1ll*ipw[i]%mod;
for(int i=1;i<=n;i++)printf("%d ",x[i]);printf("
");
}
}
