Description
Solution
如题，有
[sum_{j = 1} ^ n gcd(i, j) ^ c imes lcm(i, j) ^ d imes x_j equiv b_i pmod p
]
首先先把(lcm(i, j))用(frac{i imes j}{gcd(i, j)})替换，得到
[sum_{j = 1}^n gcd(i, j)^{c - d} imes i^d imes j^d imes x_j equiv b_i pmod p
]
设(h(x) = x ^ d)，(f(x) = x ^ {c - d})，那么将它们带入，得到
[sum_{j = 1}^n f(gcd(i, j)) imes h(i) imes h(j) imes x_j equiv b_i pmod p
]
式子中的(f(gcd(i, j)))看起来不大好处理，那么假设有(f(n) = sum_{d | n} fr(d))，就可以将(fr)带入替换(f)，得到
[sum_{j = 1}^n sum_d [d | i] [d | j] fr(d) imes h(i) imes h(j) imes x_j equiv b_i pmod p
]
我们发现这样就好处理很多，所以只要求(fr)就行，那么求(fr)就可以直接莫反。
交换求和号，将([d|i])扔进(sum)里，得到
[sum_{d | i} fr(d) sum_{j = 1}^n [d | j] h(i) imes h(j) imes x_j equiv b_i pmod p
]
设(z(d) = sum_{d | j, j <= n} h(j) imes x_j)，将(h(i))移到等号右边，得到
[sum_{d | i} fr(d) imes z(d) equiv frac{b_i}{h(i)} pmod p
]
那么发现这也是莫反的形式，可以莫反一次求出(fr(i) imes z(d))的值。
然后都(div fr(i))就得到了(z(i))的值，发现(z)的形式也是莫反，就可以求出(h(i) * x_i)的值，再(div h(i))就能得到答案。
注意因为不能有正数除(0)，所以要判断无解的情况。
复杂度的瓶颈在于莫反，最快大概能做到(O(n log log n))，懒得写了，直接(O(n log n))滚粗。
注意(c - d)可能很大，因为模数是质数所以可以将(c - d mod (MOD - 1))
Code
#include <iostream>
#include <cstdio>
#include <cstring>
const int N = 100000;
const int MOD = 998244353;
int n, c, d, q, b[N + 50], tmp[N + 50], h[N + 50], p, f[N + 50], fr[N + 50], invfr[N + 50], invh[N + 50];
void Read(int &x)
{
x = 0; int p = 0; char st = getchar();
while (st < '0' || st > '9') p = (st == '-'), st = getchar();
while (st >= '0' && st <= '9') x = (x << 1) + (x << 3) + st - '0', st = getchar();
x = p ? -x : x;
return;
}
int Ksm(int a, int b)
{
int tmp = 1;
while (b)
{
if (b & 1) tmp = 1LL * tmp * a % MOD;
a = 1LL * a * a % MOD;
b >>= 1;
}
return tmp;
}
int main()
{
Read(n); Read(c); Read(d); Read(q); p = (c - d + MOD - 1) % (MOD - 1);
for (int i = 1; i <= n; i++) h[i] = Ksm(i, d), invh[i] = Ksm(h[i], MOD - 2), f[i] = Ksm(i, p);
for (int i = 1; i <= n; i++) fr[i] = f[i];
for (int i = 1; i <= n; i++)
for (int j = i + i; j <= n; j += i)
fr[j] = (fr[j] - fr[i] + MOD) % MOD;
for (int i = 1; i <= n; i++) invfr[i] = Ksm(fr[i], MOD - 2);
while (q--)
{
int flag = 0;
for (int i = 1; i <= n; i++) Read(b[i]);
for (int i = 1; i <= n; i++)
{
if (b[i] && !h[i]) flag = 1;
b[i] = 1LL * b[i] * invh[i] % MOD;
}
for (int i = 1; i <= n; i++) tmp[i] = b[i];
for (int i = 1; i <= n; i++)
for (int j = i + i; j <= n; j += i)
tmp[j] = (tmp[j] - tmp[i] + MOD) % MOD;
for (int i = 1; i <= n; i++)
{
if (tmp[i] && !fr[i]) flag = 1;
tmp[i] = 1LL * tmp[i] * invfr[i] % MOD;
b[i] = 0;
}
for (int i = 1; i <= n; i++) b[i] = tmp[i];
for (int i = n; i >= 1; i--)
for (int j = i + i; j <= n; j += i)
b[i] = (b[i] - b[j] + MOD) % MOD;
for (int i = 1; i <= n; i++)
{
if (b[i] && !h[i]) flag = 1;
b[i] = 1LL * b[i] * invh[i] % MOD;
}
if (flag) { puts("-1"); continue; }
for (int i = 1; i <= n; i++) printf("%d ", b[i]);
printf("
");
}
return 0;
}