多項式與點值式
正常\(\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)\)
多項式拉格朗日插值的式子是
那么需要快速求出\(\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
#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);
}