Residual Polynomial
写出所有的 f i ( x ) f_i(x) fi(x)出来, f i , j f_{i, j} fi,j表示 f i ( x ) f_i(x) fi(x)的第 j j j项系数
{
f
1
,
0
f
1
,
1
f
1
,
2
…
f
1
,
n
−
1
f
1
,
n
f
2
,
0
f
2
,
1
f
2
,
2
…
f
2
,
n
−
1
f
2
,
n
f
3
,
0
f
3
,
1
f
3
,
2
…
f
3
,
n
−
1
f
3
,
n
⋮
⋮
⋮
⋱
⋮
⋮
f
n
−
1
,
0
f
n
−
1
,
1
f
n
−
1
,
2
…
f
n
−
1
,
n
−
1
f
n
−
1
,
n
f
n
,
0
f
n
,
1
f
n
,
2
…
f
n
,
n
−
1
f
n
,
n
}
\left\{ \begin{matrix} f_{1, 0} & f_{1, 1} & f_{1, 2} & \dots & f_{1, n - 1} & f_{1, n}\\ f_{2, 0} & f_{2, 1} & f_{2, 2} & \dots & f_{2, n - 1} & f_{2, n}\\ f_{3, 0} & f_{3, 1} & f_{3, 2} & \dots & f_{3, n - 1} & f_{3, n}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ f_{n - 1, 0} & f_{n - 1, 1} & f_{n - 1, 2} & \dots & f_{n - 1, n - 1} & f_{n - 1, n}\\ f_{n, 0} & f_{n, 1} & f_{n, 2} & \dots & f_{n, n - 1} & f_{n, n}\\ \end{matrix} \right\}
⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧f1,0f2,0f3,0⋮fn−1,0fn,0f1,1f2,1f3,1⋮fn−1,1fn,1f1,2f2,2f3,2⋮fn−1,2fn,2………⋱……f1,n−1f2,n−1f3,n−1⋮fn−1,n−1fn,n−1f1,nf2,nf3,n⋮fn−1,nfn,n⎭⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎫
第一行是
a
0
,
a
1
,
a
2
,
…
,
a
n
−
1
,
a
n
a_0, a_1, a_2, \dots, a_{n - 1}, a_{n}
a0,a1,a2,…,an−1,an。
考虑下一行如何从上一行转移,分两种:
一、 f i , j × c i + 1 = f i + 1 , j f_{i, j} \times c_{i + 1} = f_{i + 1, j} fi,j×ci+1=fi+1,j
二、 f i , j × j × b i + 1 = f i + 1 , j − 1 f_{i, j} \times j \times b_{i + 1} = f_{i + 1, j - 1} fi,j×j×bi+1=fi+1,j−1
可以想想成每个点有一条向下的连边 v a l u e = c value = c value=c,同时有一条向左下的连边 v a l u e = b × j value = b \times j value=b×j。
我们从 f 1 , i f_{1, i} f1,i出发,走向 f n , j f_{n, j} fn,j,显然选择了 i − j i- j i−j个 b × j b \times j b×j,所以还有 n − 1 − ( i − j ) n - 1 - (i - j) n−1−(i−j)个 c c c。
先不考虑 j j j,只算 b , c b, c b,c的贡献,设选 b b b的生成函数为 F ( x ) F(x) F(x),对于某个点有 F ( x ) = c + b x F(x) = c + bx F(x)=c+bx,因为上面的转移是具有区间性的。
所以 F ( l , r , x ) = F ( l , m i d , x ) ∗ F ( m i d , r , x ) F(l, r, x) = F(l, mid, x) * F(mid, r, x) F(l,r,x)=F(l,mid,x)∗F(mid,r,x),也就是在左区间选几个,右区间选几个,然后组合一下,这个组合刚好满足卷积的形式。
这一步就可以递归 O ( n log n log n ) O(n \log n \log n) O(nlognlogn)求解了。
接下来考虑 j j j的影响:
f n , k = ∑ i − k = j F ( j ) f 1 , i i ! k ! f_{n, k} = \sum\limits_{i - k = j} F(j)f_{1, i} \frac{i!}{k!} fn,k=i−k=j∑F(j)f1,ik!i!
这里可以简单理解一下,每次只能向左向下移动,k是我们最后在的列,只能从初始列 i ≥ k i \geq k i≥k的列转移过来。
( k ! × f n , k ) = ∑ i − k = j F ( j ) ( f 1 , i × i ! ) \left(k! \times f_{n, k}\right) = \sum\limits_{i - k = j} F(j) \left(f_{1, i} \times i!\right) (k!×fn,k)=i−k=j∑F(j)(f1,i×i!)
设 H ( k ) = k ! × f n , k , G ( i ) = f 1 , i × i ! 设H(k) = k! \times f_{n, k}, G(i) = f_{1, i} \times i ! 设H(k)=k!×fn,k,G(i)=f1,i×i!, 考 虑 把 G 翻 转 , G ( n − i ) = G ( i ) 考虑把G翻转,G(n - i) = G(i) 考虑把G翻转,G(n−i)=G(i),
H ( k ) = ∑ i + j = n − k F ( j ) G ( i ) H(k) = \sum_{i + j = n - k} F(j) G(i) H(k)=∑i+j=n−kF(j)G(i),于是再做一次卷积即可求得答案。
#include <bits/stdc++.h>
using namespace std;
const int mod = 998244353, inv2 = mod + 1 >> 1;
namespace Quadratic_residue {
struct Complex {
int r, i;
Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
};
int I2;
Complex operator * (const Complex &a, Complex &b) {
return Complex((1ll * a.r * b.r % mod + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
}
Complex quick_pow(Complex a, int n) {
Complex ans = Complex(1, 0);
while (n) {
if (n & 1) {
ans = ans * a;
}
a = a * a;
n >>= 1;
}
return ans;
}
int get_residue(int n) {
mt19937 e(233);
if (n == 0) {
return 0;
}
if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
return -1;
}
uniform_int_distribution<int> r(0, mod - 1);
int a = r(e);
while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
a = r(e);
}
I2 = (1ll * a * a % mod - n + mod) % mod;
int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
if(x > y) swap(x, y);
return x;
}
}
const int N = 1e6 + 10;
int r[N], inv[N], b[N], c[N], d[N], e[N], t[N];
int quick_pow(int a, int n) {
int ans = 1;
while (n) {
if (n & 1) {
ans = 1ll * a * ans % mod;
}
a = 1ll * a * a % mod;
n >>= 1;
}
return ans;
}
void get_r(int lim) {
for (int i = 0; i < lim; i++) {
r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
}
}
void get_inv(int n) {
inv[1] = 1;
for (int i = 2; i <= n; i++) {
inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
}
}
void NTT(int *f, int lim, int rev) {
for (int i = 0; i < lim; i++) {
if (i < r[i]) {
swap(f[i], f[r[i]]);
}
}
for (int mid = 1; mid < lim; mid <<= 1) {
int wn = quick_pow(3, (mod - 1) / (mid << 1));
for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
int w = 1;
for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
}
}
}
if (rev == -1) {
int inv = quick_pow(lim, mod - 2);
reverse(f + 1, f + lim);
for (int i = 0; i < lim; i++) {
f[i] = 1ll * f[i] * inv % mod;
}
}
}
void polyinv(int *f, int *g, int n) {
if (n == 1) {
g[0] = quick_pow(f[0], mod - 2);
return ;
}
polyinv(f, g, n + 1 >> 1);
for (int i = 0; i < n; i++) {
t[i] = f[i];
}
int lim = 1;
while (lim < 2 * n) {
lim <<= 1;
}
get_r(lim);
NTT(t, lim, 1);
NTT(g, lim, 1);
for (int i = 0; i < lim; i++) {
int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
g[i] = 1ll * g[i] * cur % mod;
t[i] = 0;
}
NTT(g, lim, -1);
for (int i = n; i < lim; i++) {
g[i] = 0;
}
}
void polysqrt(int *f, int *g, int n) {
if (n == 1) {
g[0] = Quadratic_residue::get_residue(f[0]);
return ;
}
polysqrt(f, g, n + 1 >> 1);
polyinv(g, b, n);
int lim = 1;
while (lim < 2 * n) {
lim <<= 1;
}
get_r(lim);
for (int i = 0; i < n; i++) {
t[i] = f[i];
}
NTT(g, lim, 1);
NTT(b, lim, 1);
NTT(t, lim, 1);
for (int i = 0; i < lim; i++) {
g[i] = (1ll * inv2 * g[i] % mod + 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;
b[i] = t[i] = 0;
}
NTT(g, lim, -1);
for (int i = n; i < lim; i++) {
g[i] = 0;
}
}
void derivative(int *a, int *b, int n) {
for (int i = 0; i < n; i++) {
b[i] = 1ll * a[i + 1] * (i + 1) % mod;
}
}
void integrate(int *a, int n) {
for (int i = n - 1; i >= 1; i--) {
a[i] = 1ll * a[i - 1] * inv[i] % mod;
}
a[0] = 0;
}
void polyln(int *f, int *g, int n) {
polyinv(f, b, n);
derivative(f, g, n);
int lim = 1;
while (lim < 2 * n) {
lim <<= 1;
}
get_r(lim);
NTT(g, lim, 1);
NTT(b, lim, 1);
for (int i = 0; i < lim; i++) {
g[i] = 1ll * g[i] * b[i] % mod;
b[i] = 0;
}
NTT(g, lim, -1);
for (int i = n; i < lim; i++) {
g[i] = 0;
}
integrate(g, n);
}
void polyexp(int *f, int *g, int n) {
if (n == 1) {
g[0] = 1;
return ;
}
polyexp(f, g, n + 1 >> 1);
int lim = 1;
while (lim < 2 * n) {
lim <<= 1;
}
polyln(g, d, n);
for (int i = 0; i < n; i++) {
t[i] = (f[i] - d[i] + mod) % mod;
}
t[0] = (t[0] + 1) % mod;
get_r(lim);
NTT(g, lim, 1);
NTT(t, lim, 1);
for (int i = 0; i < lim; i++) {
g[i] = 1ll * g[i] * t[i] % mod;
t[i] = d[i] = 0;
}
NTT(g, lim, -1);
for (int i = n; i < lim; i++) {
g[i] = 0;
}
}
/*
b存放多项式逆,
c存放多项式开根,
d存放多项式对数ln,
e存放多项式指数exp,
t作为中间转移数组,
如果要用到polyln,得提前调用get_inv(n)先预先得到我们想要得到的逆元范围。
*/
int B[N], C[N], f[N], fac[N], inv1[N], n;
vector<int> F[N];
int A1[N], B1[N];
void init() {
fac[0] = 1;
for (int i = 1; i < N; i++) {
fac[i] = 1ll * fac[i - 1] * i % mod;
}
inv1[N - 1] = quick_pow(fac[N - 1], mod - 2);
for (int i = N - 2; i >= 0; i--) {
inv1[i] = 1ll * inv1[i + 1] * (i + 1) % mod;
}
}
void merge(int rt) {
int ls = rt << 1, rs = rt << 1 | 1;
int n = F[ls].size(), m = F[rs].size();
F[rt].resize(n + m);
for (int i = 0; i < n + m; i++) {
F[rt][i] = 0;
}
if (n <= 50 && m <= 50) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
F[rt][i + j] = (F[rt][i + j] + 1ll * F[ls][i] * F[rs][j] % mod) % mod;
}
}
return ;
}
for (int i = 0; i < n; i++) {
A1[i] = F[ls][i];
}
for (int i = 0; i < m; i++) {
B1[i] = F[rs][i];
}
int lim = 1;
while (lim < n + m) {
lim <<= 1;
}
get_r(lim);
NTT(A1, lim, 1);
NTT(B1, lim, 1);
for (int i = 0; i < lim; i++) {
A1[i] = 1ll * A1[i] * B1[i] % mod;
}
NTT(A1, lim, -1);
for (int i = 0; i < n + m; i++) {
F[rt][i] = A1[i];
}
for (int i = 0; i < lim; i++) {
A1[i] = B1[i] = 0;
}
}
void divide(int rt, int l, int r) {
if (l == r) {
F[rt].push_back(C[l]);
F[rt].push_back(B[l]);
return ;
}
int mid = l + r >> 1;
divide(rt << 1, l, mid), divide(rt << 1 | 1, mid + 1, r);
merge(rt);
F[rt << 1].clear(), F[rt << 1 | 1].clear();
}
int main() {
// freopen("in.txt", "r", stdin);
// freopen("out.txt", "w", stdout);
// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
int T;
init();
scanf("%d", &T);
while (T--) {
scanf("%d", &n);
for (int i = 0; i <= n; i++) {
scanf("%d", &f[i]);
f[i] = 1ll * f[i] * fac[i] % mod;
}
for (int i = 0; i <= n - 2; i++) {
scanf("%d", &B[i]);
}
for (int i = 0; i <= n - 2; i++) {
scanf("%d", &C[i]);
}
divide(1, 0, n - 2);
int lim = 1;
while (lim <= 2 * n) {
lim <<= 1;
}
get_r(lim);
for (int i = 0; i <= n; i++) {
A1[i] = f[n - i];
}
for (int i = 0; i < F[1].size(); i++) {
B1[i] = F[1][i];
}
F[1].clear();
NTT(A1, lim, 1);
NTT(B1, lim, 1);
for (int i = 0; i < lim; i++) {
A1[i] = 1ll * A1[i] * B1[i] % mod;
}
NTT(A1, lim, -1);
for (int i = 0; i <= n; i++) {
printf("%lld%c", 1ll * inv1[i] * A1[n - i] % mod, i == n ? '\n' : ' ');
}
for (int i = 0; i < lim; i++) {
A1[i] = B1[i] = 0;
}
}
return 0;
}