多项式与点值式

正常$\text{DFT/IDFT}$是构造一个特殊的点值式，即$x_i=\omega_{n}^i$

如果能通过题目条件构造出来这样的点值，就可以直接$\text{DFT/IDFT}$

那如果不能的话。。。。。

多项式多点求值

一个多项式$F(x)$我们求它在$x_0,x_0,\cdots x_{m-1}$上的点值

核心是分治+多项式取模，因此常数很大

对于当前分治区间$[l,r]\in[0,m-1]$

需要快速构造一个长度为$\frac{r-l+1} {2}$的等价多项式进入分治区间

令$G_{l,r}(x)=\prod_l^r(1-x_i)$

由于$G_{l,r(x_l)}=\cdots=G_{l,r}(x_r)=0$

所以可以将$F(x)$对于$G_{l,mid}(x)$和$G_{mid+1,r}(x)$分别取模之后得到两个等价式

递归到$[l=r]$时，$F(x)$只剩下常数项

需要被访问的$G(x)$可以预先跑一遍分治NTT求出

那么复杂度就是$O(n\log ^2n)$

这种做法代码实现困难，而且常数非常大

多项式快速插值

对于点对$(x_i,y_i)$

多项式拉格朗日插值的式子是

$$\begin{aligned}F(x) = \sum_{i=0}^{n-1} y_i \prod_{i\ne j} \frac{x-x_j} {x_i-x_j}\end{aligned}$$

那么需要快速求出$\prod \frac{1} {x_i-x_j}$

构造多项式$G(x)=\prod (x-x_i)$

那么$\prod (x_i-x_j)=\frac{G} {x-x_i}(x_i)$

由于$G(x),x-x_i$在$x_i$上的点值均为$0$

我们要求的多项式就是$\begin{aligned} \prod_{i\ne j} (x_i-x_j) \end{aligned}=\frac{G(x)} {x-x_i}$

即求出$\frac{G} {x-x_i}(x_i)$

分母分子均为$0$，所以带入洛必达法则$\begin{aligned}\frac{G} {x-x_i}(x_i)=\frac{G'} {(x-x_i)'}(x_i)=G'(x_i)\end{aligned}$

那么求出$G’(x)$，然后多项式多点求值即可

剩下那一部分的答案，可以简单地分治合并上来，$[l=r]$时，多项式是一个常数

合并上来时

$[l,mid]$的答案补上$\prod_{mid+1}^r (x-x_i)$

$[mid+1,r]$的答案补上$\prod_{l}^{mid} (x-x_i)$

即复杂度为$O(n\log ^2n)$

垃圾模板题卡常

应用转置原理对于多点求值的优化

由于这个东西实在是太新了，所以没有什么文献可以看

关于转置原理的前置定义

矩阵的转置:

对于$n\cdot m$的矩阵$M$，它的转置$M^{T}$为交换行列坐标后得到的$m\cdot n$的矩阵

其满足运算性质:

1.逆: ${(A^T)}^T=A$

2.和:$(A+B)^T=A^T+B^T$

3.反积:$(AB)^T=B^TA^T$

初等矩阵:

初等矩阵是指单位矩阵通过初等变换(交换行列，某一行(列)乘上$k$加到另一行(列)上，类似高斯消元)得到的矩阵

对于计算$b=A\cdot a$，其中$A$为矩阵，$a,b$为列向量

考虑先计算$b’=A^T\cdot a$

出计算$b’$的过程，这可以分解成若干步操作（或说是初等矩阵）$E_1,E_2,\cdots E_k$

即$b’=E_1\cdot E_2\cdot E_3\cdots E_k\cdot a$

将$E_i$倒序执行，并且每一步都换成原先操作的转置$E_i^T$，就能得到$A\cdot a$

即$b=E^T_k\cdot E^T_{k-1}\cdots E^T_1\cdot a$


应用转置原理的优化核心

如果把多项式系数视为列向量$F$，则可以把多项式多点求值的过程视为一个矩阵运算$M$

为了便于描述，设要求的点值和多项式项数均为$n$

设要求的点横坐标为$x_i$，则$M$是范德蒙德矩阵

$1$ $x_0^1$ $x_0^2$ …
1 $x_1^1$ $x_1^2$ …
1 $x_2^1$ $x_2^2$ …
…

分析会发现我们要求的实际上是$b=M\cdot F$（到底是谁对矩阵乘法有误解？）

现在来将问题转置，先假装求$b’=M^T\cdot F$

$1$ 1 1 …
$x_0^1$ $x_1^1$ $x_2^1$ …
$x_0^2$ $x_1^2$ $x_2^2$ …
…

实际$M^T\cdot F$得到的结果用形式幂级数表示是

$\displaystyle\sum F_i\sum_{j=0}^{n-1}x_i^j\equiv \sum \frac{F_i} {1-x_ix}\pmod {x^n}$

求$\displaystyle M^T\cdot F= \sum \frac{F_i} {1-x_ix}\pmod {x^n}$

可以用两次分治 $\text{NTT}$ 解决，大致过程可以描述为

1.将问题转化为求$\begin{aligned} \frac{\sum F_i\prod _{i\ne j} {(1-x_jx)} } {\prod (1-x_ix)}\end{aligned} $

2.对于分治节点$[L,R]$，求得$T(L,R)=\prod_{i=L}^R{(1-x_i)}$

3.从下往上合并，每次合并答案为$A(L,R)=A(L,mid)\cdot T(mid+1,R)+A(mid+1,R)\cdot T(L,mid)$

4.最后将答案$A(0,n-1)$除以$\prod(1-x_ix)$

然后我们考虑把所有的操作都反过来并且换成转置，求得$M\cdot F$

因为过程中涉及到多项式卷积，设其转置运算为$\oplus$

我们知道普通的多项式卷积为$F(x)\cdot G(x)=\sum_i\sum_j [x^i]F(x)[x^j]G(x)x^{i+j}$

则其转置为$mul^T(F(x),G(x))=F(x)\oplus G(x)=\sum_i\sum_{j\leq i} [x^i]F(x)[x^j]G(x)x^{i-j}$

可以看到这个操作会导致多项式项数降低，若原先$F(x),G(x)$最高项为$n,m$，则转置卷积后最高项为$n-m$

那么给出整个转置后的过程为

1.在$F(x)$后面加上若干个$0$，求出$\displaystyle A(0,n-1)=F(x) \oplus \frac{1} {\prod(1-x_ix)}$的前$n$项

2.对于每个分治节点，依然预处理$\displaystyle T(L,R)=\prod_{i=L}^R{(1-x_ix)}$

3.从顶向下递归，向子节点下传

$A(L,mid)= A(L,R)\oplus T(mid+1,R)$

$A(mid+1,R)= A(L,R)\oplus T(L,mid)$

递归到子节点时，只剩一项，即是每一个点值

关于这个优化的效果:

1.不需要写多项式除法和取模了!

2.第二次分治的过程中调用的$mul^T$长度短一倍

下面这份代码是优化过的版本，能快一倍左右，但关键还是代码短听说可以被卡常卡到1s跑1e6

Luogu Submission

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
#define pb push_back
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=1<<17,P=998244353;
typedef vector <int> V;

int n,m;
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int w[N],Inv[N+1],rev[N];
void Init(){
	w[N/2]=1;
	for(int t=qpow(3,(P-1)/N),i=N/2+1;i<N;++i) w[i]=1ll*w[i-1]*t%P;
	drep(i,N/2-1,1) w[i]=w[i<<1];
	Inv[0]=Inv[1]=1;
	rep(i,2,N) Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
}
int Init(int n) {
	int R=1,cc=-1;
	while(R<n) R<<=1,cc++;
	rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
	return R;
}

void NTT(int n,V &A,int f){
    ull a[N];
	if((int)A.size()<n) A.resize(n);
	rep(i,0,n-1) a[i]=A[rev[i]];
	for(int i=1;i<n;i<<=1) {
		int *e=w+i;
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;++j) {
				int t=a[j+i]*e[j-l]%P;
				a[j+i]=a[j]+P-t;
				a[j]+=t;
			}
		}
	}
	rep(i,0,n-1) A[i]=a[i]%P;
	if(f==-1) {
		reverse(A.begin()+1,A.end());
		rep(i,0,n-1) A[i]=1ll*A[i]*Inv[n]%P;
	}
}

V operator ~ (V F) {
	int n=F.size();
	if(n==1) return V{(int)qpow(F[0])};
	V G=F; G.resize((n+1)/2),G=~G;
	int R=Init(n*2);
	NTT(R,F,1),NTT(R,G,1);
	rep(i,0,R-1) F[i]=(2-1ll*F[i]*G[i]%P+P)*G[i]%P;
	NTT(R,F,-1),F.resize(n);
	return F;
}

V operator * (V A,V B) {
	int n=A.size()+B.size()-1,R=Init(n);
	NTT(R,A,1),NTT(R,B,1);
	rep(i,0,R-1) A[i]=1ll*A[i]*B[i]%P;
	NTT(R,A,-1),A.resize(n);
	return A;
}

V Evaluate(V F,V X){
	static int ls[N<<1],rs[N<<1],cnt;
	static V T[N<<1];
	static auto TMul=[&](V F,V G){
		reverse(G.begin(),G.end());
		int n=F.size(),m=G.size(),R=Init(n);
		NTT(R,F,1),NTT(R,G,1);
		rep(i,0,R-1) F[i]=1ll*F[i]*G[i]%P;
		NTT(R,F,-1); V T(n-m+1);
		rep(i,0,n-m) T[i]=F[i+m-1];
		return T;
	};
	static function <int(int,int)> Build=[&](int l,int r) {
		int u=++cnt; ls[u]=rs[u]=0;
		if(l==r) {
			T[u]=V{1,P-X[l]};
			return u;
		}
		int mid=(l+r)>>1;
		ls[u]=Build(l,mid),rs[u]=Build(mid+1,r);
		T[u]=T[ls[u]]*T[rs[u]];
		return u;
	};

	int n=F.size(),m=X.size();
	cmax(n,m),F.resize(n),X.resize(n);
	cnt=0,Build(0,n-1);

	F.resize(n*2+1),T[1]=TMul(F,~T[1]);
	int p=0;
	rep(i,1,cnt) if(ls[i]) {
		swap(T[ls[i]],T[rs[i]]);
		int R=Init(T[i].size()),n=T[i].size(),m1=T[ls[i]].size(),m2=T[rs[i]].size();
		NTT(R,T[i],1);
		reverse(T[ls[i]].begin(),T[ls[i]].end()); reverse(T[rs[i]].begin(),T[rs[i]].end());
		NTT(R,T[ls[i]],1); NTT(R,T[rs[i]],1);
		rep(j,0,R-1) {
			T[ls[i]][j]=1ll*T[ls[i]][j]*T[i][j]%P;
			T[rs[i]][j]=1ll*T[rs[i]][j]*T[i][j]%P;
		}
		NTT(R,T[ls[i]],-1); NTT(R,T[rs[i]],-1);
		rep(j,0,n-m1) T[ls[i]][j]=T[ls[i]][j+m1-1];
		T[ls[i]].resize(n-m1+1);
		rep(j,0,n-m2) T[rs[i]][j]=T[rs[i]][j+m2-1];
		T[rs[i]].resize(n-m2+1);

		//T[ls[i]]=TMul(T[i],T[ls[i]]); T[rs[i]]=TMul(T[i],T[rs[i]]);
	} else X[p++]=T[i][0];
	X.resize(m);
	return X;
}

int main(){
	Init(),n=rd(),m=rd();
	V F(n+1),X(m);
	rep(i,0,n) F[i]=rd();
	rep(i,0,m-1) X[i]=rd();
	V Res=Evaluate(F,X);
	for(int i:Res) printf("%d\n",i);
}