• 你不可不读的多项式教程



    写在前面

    看懂这篇博客,需要掌握

    • 一般程度的代数知识
    • 入门程度的复数知识

    多项式乘法与卷积

    对于给定的两个多项式

    [F(x),G(x) ]

    我们称

    [S(x)=sum_ksum_{i+j==k}F(i)G(j) ]

    (F,G)的卷积。不难证明,这就是多项式乘法,即

    [F(x)G(x)=sum_ksum_{i+j==k}F(i)G(j) ]

    (简单地说,(F,G)的乘积中次数为(k)的项系数为

    [sum_{i+j==k}F(i)G(j) ]

    计算(F,G)的卷积显然有(O(n^2))的朴素思路,但是对于较高的数据范围这样的时间复杂度是无法接受的。


    多项式的点值表达

    除了常见的系数表达以外,多项式还可以用点值表达来表示。

    对于已知次数为(n)的多项式(F(x)),我们可以用

    [X={(x_1,F(x_1)),(x_2,F(x_2),...,(x_{n+1},F(x_{n+1})))} ]

    来唯一确定其系数表达式。证明可以由代数基本定理推出。

    多项式的点值表达相对于系数表达有一个优点:在计算卷积时时间复杂度缩减为(O(n))。更准确地说,有

    [S(x)={(x_1x_2,F(x_1)G(x_2))...}=F(x)G(x) ]

    但是将系数表达转换为点值表达的过程仍然是(O(n^2))的,于是我们想到对这个过程进行优化。


    单位根

    快速傅里叶变换能够在(O(nlogn))的时间内解决系数表达到点值表达的过程,是因为引入了单位根的概念。称方程

    [z^n=1 ]

    的复数根为(n)次单位根。这样的单位根一共有(n)个,记为

    [omega^k_n=e^{frac{2pi ki}{n}}(k=0,1,...,n-1) ]

    两条基本性质如下:

    [|omega_n^k|=1 ag{1} ]

    [omega_n^xomega_n^y=omega_n^{x+y} ag{2} ]

    ((1))的证明是显然的:复数根在复平面上的表示即为单位圆的半径。

    ((2))可以通过

    [omega_n^k=cosfrac{2kpi}{n}+sinfrac{2kpi}{n}i ag{*} ]

    计算得出。本条性质还可以推出另外两条常用推论:

    [omega_n^k=(omega_n^1)^k ]

    [(omega_n^k)^l=omega_n^{kl} ]

    此外,单位根还有两个重要的性质

    [omega_n^k=omega_{nl}^{kl} ag{3} ]

    [omega_n^{k+n/2}=-omega_n^k ag{4} ]

    ((3))((4))的证明都可以通过直接套用((*))完成。


    快速傅里叶变换(FFT)

    快速傅里叶变换最核心的环节(之一)就是多项式的系数表达与点值表达的转换,可以采用分治的方式解完成。(此处先介绍递归写法)

    对于多项式(为什么为最高次为(n-1)次下面解释)

    [F(x)=sum_{i=0}^{n-1}a_ix^i (n=2^k) ]

    [F(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})\+(a_1x+a_3x^3+...+a_{n-1}x^{n-1}) ]

    不妨设

    [L(x)=a_0+a_2x+...+a_nx^{frac{n-2}{2}}\ R(x)=a_1+a_3x+...+a_n-1x^{frac{n-2}{2}} ]

    那么有

    [F(x)=L(x^2)+xR(x^2) ]

    (此时可以看出为什么此前将(F)最高次设为(n-1):美观)

    现在我们要求(F)的点值表达,于是代入单位根,也就是

    [F(omega_n^k)=L((omega_n^k)^2)+omega_n^kR((omega_n^k)^2)\ =L(omega_{n/2}^k)+omega_n^kR(omega_{n/2}^k) ]

    同理,可以得到

    [F(omega_n^{k+n/2})=L((omega_n^{k+n/2})^2)+omega_n^{k+n/2}R((omega_n^{k+n/2})^2)\ =L(omega_n^{2k+n})-omega_n^kR(omega_n^{2k+n})\ =L(omega_n^{2k})-omega_n^kR(omega_n^{2k})\ =L(omega_{n/2}^k)-omega_n^kR(omega_{n/2}^k) ]

    也就是说,如果我们求出了(L,R)(omega_{n/2}^0,...,omega_{n/2}^{n/2-1})的点值表达,我们就可以得到(F)(omega_{n}^0,...,omega_{n}^n-1)的点值表达。

    由于每一次长度减半,所以可以在(O(nlogn))的时间内求出多项式(F)的点值表示。

    对于多项式的点值表达到系数表达的转换,可以推出类似的公式。

    (FFT)的递归实现如下:

    const double pi=acos(-1);
    typedef complex<double> cx;
    	
    void fft (vector<cx> x,int n,int type) {
    	if (n==1) return;
    	vector<cx> l,r;
    	for (register int i=0; i<n; i+=2) 
    		l.push_back(x[i]),r.push_back(x[i+1]);
    	fft(l,n>>1,type),fft(r,n>>1,type);
    	cx wn(cos(2*pi/n),type*sin(2*pi/n)),w(1,0);
    	for (register int i=0; i<(n>>1); ++i,w*=wn) {
    		x[i]=l[i]+w*r[i],x[i+(n>>1)]=l[i]-w*r[i];
    	}
    }
    

    在得到了递归实现以后,我们发现其常数较大,因此考虑是否存在非递归实现。观察递归过程中的数组下标,其变换过程如下

    [0~~~1~~~2~~~3~~~4~~~5~~~6~~~7\0~~~2~~~4~~~6~|~1~~~3~~~5~~~7\0~~~4~|~2~~~6~|~1~~~5~|~3~~~7\0~|~4~|~2~|~6~|~1~|~5~|~3~|~7 ]

    观察其二进制表示

    [000~~~001~~~010~~~011~~~100~~~101~~~110~~~111\000~~~010~~~100~~~110~|~001~~~011~~~101~~~111\000~~~100~|~010~~~110~|~001~~~101~|~011~~~111\000~|~100~|~010~|~110~|~001~|~101~|~011~|~111 ]

    不难发现每一项递归之后的最终位置为其二进制表示翻转的结果。

    那么我们可以通过预处理一个(R)数组来直接得到最终位置。不难发现,每次递归实质上是按照对应位进行排序,于是我们可以得到实现:

    for (register int i=0; i<n; ++i) {
    	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    }
    

    那么根据(R),我们可以写出(FFT)的非递归实现:

    void fft (vector<cx> x,int n,int type) {
    	for (register int i=0; i<n; ++i)  // 求出递归后的数组
    		if (i<R[i]) swap(x[i],x[R[i]]); // 如果不判断就会交换两次,等于没交换
    	for (register int i=1; i<n; i<<=1) { // 枚举区间长度的一半(当然也可以枚举区间长度)
    		cx wn(cos(pi/i),type*sin(pi/i)); // 此处消掉了公式里的2
    		for (register int j=0; j<n; j+=(i<<1)) { // 枚举区间,对每一个区间做修改
    			cx w(1,0);
    			for (register int k=0; k<i; ++k,w*=wn) { // 枚举区间左半侧,同时修改右半侧
    				cx t=x[j+k],s=x[j+k+i]*w;
    				x[j+k]=t+s,x[j+k+i]=t-s;
    			}
    		}
    	}
    }
    

    快速数论变换(NTT)

    (NTT)(FFT)的唯一区别在于将单位根改为了原根,从而有效避免了复数取模运算中的精度丢失。原根满足此前提到的单位根的性质,其细节不做赘述。实现如下:

    void ntt (vector<int> x,int n,int type) {
    	for (register int i=0; i<n; ++i) 
    		if (i<R[i]) swap(x[i],x[R[i]]);
    	for (register int i=1; i<n; i<<=1) {
    		int wn=ksm(type==1?G:invG,(MOD-1)/(i<<1));
    		for (register int j=0; j<n; j+=(i<<1)) {
    			int w=1;
    			for (register int k=0; k<i; ++k,w*=wn) {
    				int t=x[j+k],s=(x[j+k+i]*w)%MOD;
    				x[j+k]=(t+s)%MOD,x[j+k+i]=(t-s)%MOD;
    			}
    		}
    	}
    }
    

    (G)为模数的原根,(invG)为该原根在模数意义下的逆元。对于很大一部分题目(准确地说,对于满足(MOD=a*2^b+1)的题目),我们能够直接由暴力求出其原根,对于另一类任意模数题目,则需要利用中国剩余定理进行求解。暴力实现如下:

    void getG () {
    	int tmp=MOD;
    	for (register int i=2; i*i<=MOD; ++i) {
    		if (tmp%i==0) {
    			fac[++cnt]=i;
    			while (tmp%i==0) tmp/=i;
    		}
    	}
    	if (tmp!=1) fac[++cnt]=tmp;
    	for (register int i=2; i<MOD; ++i) {
    		int flag=1;
    		for (register int j=1; j<=cnt; ++j) {
    			if (ksm(i,(MOD-1)/fac[j])==1) {
    				flag=0;
    				break;
    			}
    		}
    		if (flag) {
    			write(i);
    			break;
    		}
    	}
    }
    

    任意模数NTT暂时先坑着,回头写完多项式求逆和开根再填。


    多项式求逆

    对于(n-1)次多项式(F),其逆元为满足

    [F(x)G(x)equiv1(mod~x^n) ]

    的多项式(G)

    假设我们已经求出满足

    [F(x)G^{'}(x)equiv1(mod~x^{frac{n}{2}}) ]

    (G^{'}),由于显然有

    [F(x)G(x)equiv1(mod~x^n) ]

    那么有

    [F(x)(G(x)-G^{'}(x))equiv0(mod~x^{frac{n}{2}}) ]

    如果(F(x)equiv0(mod~x^{frac{n}{2}})),结果是显然的。那么考虑

    [G(x)-G^{'}(x)equiv0(mod~x^{frac{n}{2}}) ]

    等式两侧同时平方,得到

    [G^2(x)+G^{'2}(x)-2G(x)G^{'}(x)equiv0(mod~x^n) ]

    两侧同乘(A(x)),有

    [G(x)+F(x)G^{'2}(x)-2F(x)G^{'}(x)equiv0(mod~x^n) ]

    所以

    [G(x)equiv G^{'}(x)(2-F(x)G^{'}(x))(mod~x^n) ]

    于是说明了我们可以用模(x^{frac{n}{2}})意义下逆元推出模(x^n)的意义下的逆元。

    非递归的多项式求逆实现如下:

    void inv (vector<int> x,vector<int> y,int n) {
    	b[0]=ksm(a[0],MOD-2);
    	int len;
    	for (len=1; len<(n<<1); len<<=1) {
    		for (register int i=0; i<len; ++i) A[i]=x[i],B[i]=y[i];
    		for (register int i=0; i<len<<1; ++i) R[i]=(R[i>>1]>>1)|((i&1)?len:0);
    		ntt(A,len<<1,1),ntt(B,len<<1,1);
    		for (register int i=0; i<len<<1; ++i) y[i]=((2-A[i]*B[i]%MOD)*B[i]%MOD+MOD)%MOD;
    		ntt(y,len<<1,1);
    		for (register int i=len; i<len<<1; ++i) y[i]=0;
    	}
    	for (register int i=0; i<len; ++i) A[i]=B[i]=0;
    	for (register int i=len; i<len<<1; ++i) y[i]=0;
    }
    

    多项式开根

    对于(n-1)次多项式(F),其开根结果为满足

    [F(x)equiv G^2(x)(mod~x^n) ]

    的多项式(G)

    假设我们已经求出满足

    [F(x)equiv G^{'2}(x)(mod~x^{frac{n}{2}}) ]

    (G^{'}),由于显然有

    [F(x)equiv G^2(x)(mod~x^n) ]

    那么有

    [G^2(x)-G^{'2}(x)equiv0(mod~x^{frac{n}{2}}) ]

    可以化简得到

    [(G(x)+G^{'}(x))(G(x)-G^{'}(x))equiv0(mod~x^{frac{n}{2}}) ]

    所以

    [G(x)-G^{'}(x)equiv0(mod~x^frac{n}{2}) ]

    两侧同时平方,得到

    [G^2(x)+G^{'2}(x)-2G(x)G^{'}(x)equiv0(mod~x^n) ]

    化简得

    [F(x)+G^{'2}(x)-2G(x)G^{'}(x)equiv0(mod~x^n) ]

    所以

    [G(x)equiv frac{F(x)+G^{'2}(x)}{2G^{'}(x)}(mod~x^n) ]

    于是说明了我们可以用模(x^{frac{n}{2}})意义下开根结果推出模(x^n)的意义下的开根结果。

  • 相关阅读:
    SpringCloud(三):服务消费以及负载均衡(RestTemplate+Ribbon)
    SpringCloud(二):服务的注册与发现(Eureka)
    SpringCloud(一):了解SpringCloud
    SpringBoot(十二):SpringBoot整合Mybatis-Plus
    SpringBoot(十一):SpringBoot整合Redis
    idea使用maven中的tomcat插件开启服务出现java.net.BindException: Address already in use: JVM_Bind :8080错误原因
    SSM框架之SpringMVC(6)异常处理及拦截器
    SSM框架之SpringMVC(5)文件上传
    SSM框架之SpringMVC(4)返回值类型及响应数据类型
    SSM框架之SpringMVC(3)常用注解
  • 原文地址:https://www.cnblogs.com/ilverene/p/11644405.html
Copyright © 2020-2023  润新知