Description：

给出 $n,c,d$ ，多次询问，每次给出 $b_1,b_2,\dots,b_n$ ，求一组 $x_1,x_2,\dots,x_n$ 满足： $$ \sum_{j=1}^n\gcd(i,j)^c\times \text{lcm}(i,j)^d\times x_j\equiv b_i\mod p $$ $1\leqslant n\leqslant 10^5$

Solution：

首先把 $\text{lcm}$ 拆开，就是： $$ \sum_{j=1}^n\gcd(i,j)^{c-d}i^dj^d\times x_j\equiv b_i $$ 我们可以把这个式子抽象一下： $$ \sum_{j=1}^nf(\gcd(i,j))g(i)h(j)\times x^j=b_i $$ 那么我们构造 $f_r(i)$ 满足： $$ \sum_{d|n}f_r(d)=f(n) $$ 这个 $f_r$ 可以像 $O(n\ln n)$ 求狄利克雷卷积那样倒着求，也就是说： $$ f_r(n)=f(n)-\sum_{d|n,d\ne n}f_r(d) $$ 于是递推就可以了。

那么式子变成了： $$ \sum_{j=1}^n\sum_{d|i,d|j}f_r(d)h(j)\times x^j=b_i/g(i) $$ 套路交换求和号： $$ \sum_{d|i}f_r(d)\sum_{j=1,d|j}^nh(j)\times x_j=b_i/g(i) $$ 设 $v(d)=f_r(d)\sum_{j=1,d|j}^nh(j)\times x_j$ ，显然有 $v\times 1=b_i/g(i)$ ，于是莫比乌斯反演即可得到 $v$ ，具体做法是像上面一样 $O(n\ln n)$ 求，再除 $f_r$ 就可以得到 $z(d)=\sum_{j=1,d|j}^nh(j)\times x_j$ ，再莫比乌斯反演即可得到 $x_j$ ，有一些细节看代码。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,q;

#define MAXN 100010

typedef long long ll;

ll b[MAXN];

ll fr[MAXN];

ll z[MAXN];

ll c,d;

#define MOD 998244353

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % MOD;

		a = a * a % MOD;

		b = b >> 1;

	}

	return res;

}

int main()

{

	scanf("%d%lld%lld%d",&n,&c,&d,&q);

	for(int i = 1;i < MAXN;++i)fr[i] = power(i,(c - d + MOD - 1) % (MOD - 1));

	for(int i = 1;i < MAXN;++i)for(int j = i + i;j < MAXN;j += i)fr[j] = (fr[j] - fr[i] + MOD) % MOD;

	while(q--)

	{

		for(int i = 1;i <= n;++i)scanf("%lld",&b[i]);

		for(int i = 1;i <= n;++i)z[i] = b[i] * power(power(i,d),MOD - 2) % MOD;

		for(int i = 1;i <= n;++i)for(int j = i + i;j <= n;j += i)z[j] = (z[j] - z[i] + MOD) % MOD;

		bool tag = false;

		for(int i = 1;i <= n;++i)if(fr[i] == 0 && z[i] != 0)tag = true;

		if(tag){puts("-1");continue;}

		for(int i = 1;i <= n;++i)z[i] = z[i] * power(fr[i],MOD - 2) % MOD;

		for(int i = n;i >= 1;--i)for(int j = i + i;j <= n;j += i)z[i] = (z[i] - z[j] + MOD) % MOD;

		for(int i = 1;i <= n;++i)z[i] = z[i] * power(power(i,d),MOD - 2) % MOD;

		for(int i = 1;i <= n;++i)printf("%lld ",(z[i] + MOD) % MOD);

		puts("");

	}

	return 0;

}