好像又很久沒寫博客了
總結成一個小札還是香啊 awa
把 @tly 的博客學了一遍然后翻譯成自己看得懂的備忘
# 目錄
目錄- # 目錄
- # 基礎多項式版子
- # 多項式多點求值
- - 問題1
- - 新的解法
- - 源代碼1
- # 多項式快速插值
- - 問題2
- - 拉格朗日插值
- - 快速插值
- - 源代碼2
- THE END
- Thank for reading!
# 基礎多項式版子
你可能會在下面看到 BASIC_POLY 這么一個命名空間,代碼如下,
點擊展開/折疊代碼
#define con(type) const type &
const int N=1<<20,MOD=998244353;
inline int add(con(int)a,con(int)b){return a+b>=MOD? a+b-MOD:a+b;}
inline int sub(con(int)a,con(int)b){return a-b<0? a-b+MOD:a-b;}
inline int mul(con(int)a,con(int)b){return int(1ll*a*b%MOD);}
inline int ina_pow(con(int)a,con(int)b){return b?mul(ina_pow(mul(a,a),b>>1),(b&1)?a:1):1;}
namespace BASIC_POLY{
int w[N+10],eta_lg[N+10],rev[N+10];
void init(){
w[0]=eta_lg[1]=0;w[1]=ina_pow(3,(MOD-1)>>20);
for(int i=2;i<N;i++){
w[i]=mul(w[i-1],w[1]);
eta_lg[i]=eta_lg[(i+1)>>1]+1;
}
}
void ntt(int *ary,con(int)n,con(int)typ){
for(int i=1;i<n;i++){
rev[i]=rev[i>>1]>>1|((i&1)*(n>>1));
if(rev[i]<i) swap(ary[i],ary[rev[i]]);
}
for(int i=1,ii=2;i<n;i<<=1,ii<<=1){
int u=N/ii;
for(int j=0;j<n;j+=ii){
int *a=ary+j,*b=ary+j+i,*p=w,q=*b;
for(int k=0;k<i;k++,a++,b++,p+=u,q=mul(*b,*p))
*b=sub(*a,q),*a=add(*a,q);
}
}
if(typ==-1){
reverse(ary+1,ary+n);
for(int i=0,ivn=MOD-(MOD-1)/n;i<n;i++)
ary[i]=mul(ary[i],ivn);
}
}
int modelize(con(int)l){return 1<<eta_lg[l];}
void polyInv(int *a,int *b,con(int)n){
if(n==1){b[0]=ina_pow(a[0],MOD-2);return;}
int m=(n+1)>>1,l=modelize(2*(m-1)+(n-1)+1);
polyInv(a,b,m);
static int tmp1[N+10],tmp2[N+10];
for(int i=0;i<n;i++) tmp1[i]=a[i];fill(tmp1+n,tmp1+l,0);
for(int i=0;i<m;i++) tmp2[i]=b[i];fill(tmp2+m,tmp2+l,0);
ntt(tmp1,l,1),ntt(tmp2,l,1);
for(int i=0;i<l;i++) tmp1[i]=mul(tmp2[i],sub(2,mul(tmp1[i],tmp2[i])));
ntt(tmp1,l,-1);
for(int i=0;i<n;i++) b[i]=tmp1[i];
}
void polyDet(int *a,int *b,con(int)n){
for(int i=1;i<n;i++) b[i-1]=mul(a[i],i);
b[n-1]=0;
}
}
# 多項式多點求值
- 問題1
給定一個最高次項為 \(n-1\) 的多項式 \(f(x)\),求 \(f(a_0),f(a_1)\cdots f(a_{m-1})\),
- 新的解法
之前的做法要用到多項式取模這種大常數計算,而且代碼還很復雜……
先定義一種所謂的“差卷積”,記作
\[f\otimes g=\sum_{i=0}\sum_{j=0}^if_ig_jx^{i-j} \]我們可以構造出 \(g(x)=\frac{1}{1-ax}\) 使得 \(h=f\otimes g\) 滿足 \([x^0]h(x)=f(a)\),
點擊展開/折疊證明
把 $g(x)$ 的閉形式展開,則 $(f\otimes g)$ 相當于
$$(f_0+f_1x+f_2x^2+\cdots)\otimes(1+ax+a^2x^2+\cdots)$$觀察其常數項,一定是 $f_ix^i\otimes a^ix^i$,也就是 $\sum a^if_i$,
這樣我們可以把一個點的求值問題轉化成卷積問題,那怎么擴展到多點求值呢?
關于差卷積,我們還需要一個性質——
\[(f\otimes g)\otimes h=f\otimes(g\cdot h) \]根據定義很容易理解,即 \(f_ig_jh_kx^{i-j-k}=f_ix^i(g_jh_k)x^{-(j+k)}\),
有了這個性質,就可以考慮用線段樹的方式分治,我們希望能在線段樹的第 \(i\) 個葉子處得到 \([x^0](f\otimes\frac{1}{1-a_ix})\),那么線段樹上區間 \([l,r]\) 就維護
\[S_{l,r}=f\otimes\Big(\prod_{i=l}^r\frac{1}{1-a_ix}\Big) \]要從 \([l,r]\) 遞推到子區間 \([l,m]\),我們需要計算:
\[S_{l,m}=S_{l,r}\otimes\Big(\prod_{i=m+1}^r{1-a_ix}\Big) \]代入即可證明,于是我們還需要對線段樹的每個節點 \([l,r]\) 先預處理出
\[T_{l,r}=\prod_{i=l}^r1-a_ix \]\(T_{l,r}\) 的次數和節點的大小是一致的,因為線段樹上每一層節點大小減半,可以直接 \(O(n\log n)\) 預處理,
似乎現在就可以解決多點求值了?但還剩了一個非常棘手的問題——\(S_{l,r}\) 是一個形式冪級數,我們應該保留多少項才能計算出葉子處的常數項?
實際上我們只需要保留 \(S_{l,r}\) 的前 \(r-l+1\) 項,歸納證明如下:
- 歸納邊界:葉子處只需要保留常數項;
- 考慮 \([l,r]\) 的子節點 \([l,m]\),\(S_{l,m}\) 需要保留前 \(m-l+1\) 項(最高項為 \(m-l\) 次),而 \(S_{l,m}\) 的計算方法是:\[S_{l,m}=S_{l,r}\otimes T_{m+1,r} \]注意到 \(T_{m+1,r}\) 最高項是 \(r-m-1\) 次的,根據差卷積的計算方式,\(S_{l,r}\) 的最高項應為 \(r-l\) 次,即只需要保留前 \(r-l+1\) 項,
這樣的話復雜度就有保證了,每層項數減半,仍然可以 \(O(n\log n)\) 解決,
- 源代碼1
點擊展開/折疊代碼
//BASIC_POLY 中主要是NTT這些東西
namespace UPPER_POLY{
typedef vector<int> vint;
//普通卷積
vint polyMul(con(vint)a,con(vint)b){
static int tmp1[N+10],tmp2[N+10];
int na=(int)a.size(),nb=(int)b.size(),nc=na+nb-1;
vint c(nc);
int l=BASIC_POLY::modelize(nc);
for(int i=0;i<na;i++) tmp1[i]=a[i];fill(tmp1+na,tmp1+l,0);
for(int i=0;i<nb;i++) tmp2[i]=b[i];fill(tmp2+nb,tmp2+l,0);
BASIC_POLY::ntt(tmp1,l,1),BASIC_POLY::ntt(tmp2,l,1);
for(int i=0;i<l;i++) tmp1[i]=mul(tmp1[i],tmp2[i]);
BASIC_POLY::ntt(tmp1,l,-1);
for(int i=0;i<nc;i++) c[i]=tmp1[i];
return c;
}
//差卷積,卷積結果保留前 nc項
vint polySubMul(con(vint)a,con(vint)b,con(int)nc){
vint c=a;reverse(c.begin(),c.end());
c=polyMul(c,b);
c.resize(a.size()),reverse(c.begin(),c.end());
c.resize(nc);
return c;
}
vint seg[N];
#define idx(l,r) (((l)+(r))|((l)!=(r)))
//預處理出 T[l,r]
void build(int *a,con(int)le,con(int)ri){
if(le==ri){
seg[idx(le,ri)].clear();
seg[idx(le,ri)].push_back(1);
seg[idx(le,ri)].push_back(sub(0,a[le]));
return;
}
int mi=(le+ri)>>1;
build(a,le,mi),build(a,mi+1,ri);
seg[idx(le,ri)]=polyMul(seg[idx(le,mi)],seg[idx(mi+1,ri)]);
}
//p即 S[l,r]
void solve(int *res,con(int)le,con(int)ri,con(vint)p){
if(le==ri){res[le]=p[0];return;}
int mi=(le+ri)>>1;
solve(res,le,mi,polySubMul(p,seg[idx(mi+1,ri)],mi-le+1));
solve(res,mi+1,ri,polySubMul(p,seg[idx(le,mi)],ri-mi));
}
void multiVal(int *f,con(int)n,int *pos,con(int)m,int *r){
static int tmp1[N+10],tmp2[N+10];
build(pos,0,m-1);
int rt=idx(0,m-1);
for(int i=0,ii=min(n,(int)seg[rt].size());i<ii;i++)
tmp1[i]=seg[rt][i];
BASIC_POLY::polyInv(tmp1,tmp2,n);
//T[0,n-1] 求個逆再和 f差卷積得到 S[0,n-1]
vint v1(n),v2(n);
for(int i=0;i<n;i++) v1[i]=f[i],v2[i]=tmp2[i];
solve(r,0,m-1,polySubMul(v1,v2,m));
}
}
# 多項式快速插值
- 問題2
給定 \(n\) 組 \((x_i,y_i)\),\(n-1\) 次多項式 \(f(x)\) 滿足 \(\forall i,f(x_i)=y_i\),求 \(f(x)\),
- 拉格朗日插值
由若干個 \(n-1\) 次多項式疊加得到 \(f(x)\),具體構造如下:
\[g_i(x)=y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \]這樣構造的目的是 \(g_i(x_j)\) 當且僅當 \(i=j\) 時 \(g_i(x_j)\neq0\),于是可以將 \(g_0(x),g_1(x)\cdots g_{n-1}(x)\) 直接相加得到 \(f(x)\),
- 快速插值
考慮對拉格朗日插值進行優化(優化計算方式以加速),觀察 \(f(x)\) 的運算式:
\[f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}=\sum_{i=0}^{n-1}y_i\frac{\prod(x-x_j)}{\prod(x_i-x_j)} \]先看看怎么計算分式的分母部分,該部分對于固定的 \(i\) 來說是一個常數,記為 \(k_i\),
\[k_i=\prod_{j\neq i}(x_i-x_j) \]根據極限的相關知識,不難證明下面這個式子是成立的:
\[k_i=\lim_{x\to x_i}\frac{\prod\limits_{j=0}^{n-1}(x_i-x_j)}{x-x_i} \]便于書寫,記 \(g(x)=\prod\limits_{i=0}^{n-1}(x-x_i)\),分式上下都趨近于 \(0\),符合洛必達法則的適用條件——由此可得 \(k_i=g'(x_i)\)
,
可以先用類似線段樹的分治(本質是啟發式合并)在 \(O(n\log n)\) 的復雜度內算出 \(g(x)\),然后 \(O(n)\) 多項式求導得到 \(g'(x)\),再多點求值就可以得到 \(k_i\) 了,
再看 \(f(x)\),我們現在可以把它寫成
\[f(x)=\sum_{i=0}^{n-1}\frac{y_i}{k_i}\prod_{j\neq i}(x-x_j) \]則需要解決后面這個式子,這其實可以用到多點求值的舊方法中的一個技巧——仍然是分治:
分治
記 $h_{l,r}(x)$:
$$h_{l,r}(x)=\sum_{i=l}^r\frac{y_i}{k_i}\prod_{j\in[l,r]}^{j\neq i}(x-x_j)$$
仍然用類似線段樹的方法分治,考慮如何從 $h_{l,m}(x)$ 和 $h_{m+1,r}(x)$ 合并到 $h_{l,r}(x)$,
$$ \begin{aligned} S_{l,r}&=\prod_{i=l}^r(x-x_i)\\ h_{l,r}(x)&=h_{l,m}S_{m+1,r}+S_{l,m}h_{m+1,r} \end{aligned} $$
這樣就可以 $O(n\log n)$ 求出答案,
總的復雜度就是 \(O(n^2+n\log n)\)
- 源代碼2
點擊展開/折疊代碼
namespace UPPER_POLY{
typedef vector<int> vint;
vint polyAdd(con(vint)a,con(vint)b){
int na,nb;
vint c(max(na=(int)a.size(),nb=(int)b.size()));
for(int i=0;i<na;i++) c[i]=a[i];
for(int i=0;i<nb;i++) c[i]=add(c[i],b[i]);
return c;
}
vint polyMul(con(vint)a,con(vint)b){
int na=(int)a.size(),nb=(int)b.size(),nc=na+nb-1,l=BASIC_POLY::modelize(nc);
static int tmp1[N+10],tmp2[N+10];
for(int i=0;i<na;i++) tmp1[i]=a[i];fill(tmp1+na,tmp1+l,0);
for(int i=0;i<nb;i++) tmp2[i]=b[i];fill(tmp2+nb,tmp2+l,0);
BASIC_POLY::ntt(tmp1,l,1),BASIC_POLY::ntt(tmp2,l,1);
for(int i=0;i<l;i++) tmp1[i]=mul(tmp1[i],tmp2[i]);
BASIC_POLY::ntt(tmp1,l,-1);
vint c(nc);
for(int i=0;i<nc;i++) c[i]=tmp1[i];
return c;
}
vint polySubMul(con(vint)a,con(vint)b,con(int)nc){
vint c=a;
reverse(c.begin(),c.end());
c=polyMul(c,b);
c.resize(a.size());
reverse(c.begin(),c.end());
c.resize(nc);
return c;
}
vint seg[N+10];
#define idx(l,r) (((l)+(r))|((l)!=(r)))
void build(int *a,con(int)le,con(int)ri){
if(le==ri){
int u=idx(le,ri);
seg[u].clear();
seg[u].push_back(1),seg[u].push_back(sub(0,a[le]));
return;
}
int mi=(le+ri)>>1;
build(a,le,mi),build(a,mi+1,ri);
seg[idx(le,ri)]=polyMul(seg[idx(le,mi)],seg[idx(mi+1,ri)]);
}
void solve1(int *res,con(int)le,con(int)ri,con(vint)p){
if(le==ri){res[le]=p[0];return;}
int mi=(le+ri)>>1;
solve1(res,le,mi,polySubMul(p,seg[idx(mi+1,ri)],mi-le+1));
solve1(res,mi+1,ri,polySubMul(p,seg[idx(le,mi)],ri-mi));
}
void multiVal(int *f,con(int)n,int *pos,con(int)m,int *res){
build(pos,0,m-1);
static int tmp1[N+10],tmp2[N+10];
int rt=idx(0,m-1);
for(int i=0,ii=(int)seg[rt].size();i<n;i++)
tmp1[i]=i<ii?seg[rt][i]:0;
BASIC_POLY::polyInv(tmp1,tmp2,n);
vint p1(n),p2(n);
for(int i=0;i<n;i++) p1[i]=f[i],p2[i]=tmp2[i];
solve1(res,0,m-1,polySubMul(p1,p2,m));
}
void build2(int *a,con(int)le,con(int)ri){
if(le==ri){
seg[idx(le,ri)].clear();
seg[idx(le,ri)].push_back(sub(0,a[le])),seg[idx(le,ri)].push_back(1);
return;
}
int mi=(le+ri)>>1;
build2(a,le,mi),build2(a,mi+1,ri);
seg[idx(le,ri)]=polyMul(seg[idx(le,mi)],seg[idx(mi+1,ri)]);
}
vint solve2(int *y,con(int)le,con(int)ri){
if(le==ri) return vint(1,y[le]);
int mi=(le+ri)>>1;
return polyAdd(polyMul(solve2(y,le,mi),seg[idx(mi+1,ri)]),polyMul(seg[idx(le,mi)],solve2(y,mi+1,ri)));
}
vint multiPoly(int *pos,int *val,int n){
static int tmp1[N+10],tmp2[N+10];
build2(pos,0,n-1);
for(int i=0,rt=idx(0,n-1);i<=n;i++) tmp1[i]=seg[rt][i];
BASIC_POLY::polyDet(tmp1,tmp1,n+1);
multiVal(tmp1,n,pos,n,tmp2);
build2(pos,0,n-1);
//getinv
tmp1[0]=tmp2[0];for(int i=1;i<n;i++) tmp1[i]=mul(tmp1[i-1],tmp2[i]);
tmp1[n-1]=ina_pow(tmp1[n-1],MOD-2);
for(int i=n-2;~i;i--){
int tmp1i=mul(tmp1[i+1],tmp2[i+1]);
tmp1[i+1]=mul(tmp1[i+1],tmp1[i]);
tmp1[i]=tmp1i;
}
for(int i=0;i<n;i++) tmp1[i]=mul(tmp1[i],val[i]);
return solve2(tmp1,0,n-1);
}
}
THE END
Thank for reading!
若陪伴是我所有
那痕跡 藏匿了整個宇宙
——《聽風捕夢》By (×28)雙笙/封茗囧菌/司南
> Link 聽風捕夢-網易云
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/261654.html
標籤:其他
