「NOTE」進階多項式小札

時間  2021-02-22
標籤 優化 spa 翻譯 code get 博客 ast class 简体版
原文   原文鏈接

好像又好久沒寫博客了
總結成一個小札仍是香啊 awa
把 @tly 的博客學了一遍而後翻譯成本身看得懂的備忘優化

# 目錄

目錄

# 基礎多項式版子

你可能會在下面看到 BASIC_POLY 這麼一個命名空間，代碼以下。ui

點擊展開/摺疊代碼 
#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})\)。spa

- 新的解法

以前的作法要用到多項式取模這種大常數計算，並且代碼還很複雜……翻譯

先定義一種所謂的「差卷積」，記做code

\[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)\)。get

點擊展開/摺疊證實

把 $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$。it

這樣咱們能夠把一個點的求值問題轉化成卷積問題。那怎麼擴展到多點求值呢？ast

關於差卷積，咱們還須要一個性質——class

\[(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\) 項，概括證實以下：

  1. 概括邊界：葉子處只須要保留常數項；
  2. 考慮 \([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 聽風捕夢-網易雲

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程