zoukankan      html  css  js  c++  java
  • 多项式:从门都没入到刚迈过门槛

    多项式:从门都没入到刚迈过门槛

    还记得我上次讲多项式的文章么?经过CSP的洗礼和最近一段时间的学习后,我对多项式的理解又加深了吧。

    初三党文化课压力是大,尽管上次月考年级第一,但是文化课还是不能掉以轻心,毕竟要中考,也不能没学上吧。

    零、快速数论变换(NTT)

    多项式乘法是多项式其他操作的一个基础。上次我们研究了一下FFT。但是不少题目都是要求取模的。显然,FFT中的复数是不能取模的。而且,很多时候题目中多项式的系数也都是整数。这个时候,NTT就派上用场了!

    前置知识:原根

    (其实我本来是打算写一篇博客讲这玩意的,但是被我咕了)

    阶:当 ((a,m)=1)时,令最小的(x)使得(a^xequiv 1pmod{m})(apmod{m})的阶。

    有一个很重要的性质:(x|varphi(m))

    (a)的取值为(g)时,如果它的阶就是(varphi(m)),那么我们称(g)(m)的原根。

    原根怎么求呢?考虑到大部分情况下原根都比较小,我们只需从2开始枚举原根,然后从(1)(varphi(m))都试一下就行了。

    常见的模数(998244353,1004535809,469762049)的原根都是(3)

    下面是关于原根的几个结论。我这里就不证明了。(毕竟重点在NNT不在原根上)

    1.(g^0,g^1,g^2,dots,g^{varphi(m)})两两不同,构成了(m)的完全剩余系。

    2.只有(2,4,p^k,2p^k)存在原根。

    3.若(g)(p)的原根,(g)也是(p^k)原根。

    我们回顾一下FFT的过程。FFT结合了单位根的美妙性质,完成了蝴蝶操作,从而使复杂度降低为(O(n log n))

    而在模意义下,原根也具有单位根的性质!

    我们只需要把FFT中的(pi)换成模数,把单位根换成原根,就是NTT了!

    核心代码:

    inline void NTT(int *a,int len,int f){
    	for(int i=0;i<len;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    	for(int mid=1;mid<len;mid<<=1){//枚举合并区间的一半长度
    		int wn=ksm((f==1?g:gi),(MOD-1)/(mid<<1));//求当前的原根
    		for(int j=0;j<len;j+=(mid<<1)){//进行蝴蝶操作
    			int w=1;
    			for(int k=j;k<j+mid;k++){
    				int u=a[k];
    				int t=(1ll*w*a[k+mid])%MOD;
    				a[k]=(u+t)%MOD;
    				a[k+mid]=(u-t+MOD)%MOD;
    				w=(1ll*w*wn)%MOD;
    			}
    		}
    	} 
    }
    

    那么NTT和FFT的效率对比?

    1

    1

    好吧,至少在这道题里,以及我写的代码里,NTT是完爆FFT的。

    在我们开始多项式的其他操作之前,我们先给大家讲个概念:卷积

    设有两个数列(a)(b),那么这两个数列的卷积(c)定义为:

    [c_k=sum_{i=0}^ka_ib_{k-i} ]

    我们常常记作:

    [c=a*b ]

    我们很容易的得到,多项式的乘法与多项式的卷积是一致的:

    (F(x)=sum a_i x^i)(G(x)=sum b_ix^i)(A(x)=sum c_ix^i),且(c=a*b),于是我们可以得到:

    [F=A*G ]

    下面就让我们看看多项式的一些其他操作吧:

    一、多项式求逆

    已知一个多项式(G(x)),求一个多项式(F(x)),使得:

    [F(x)G(x)equiv 1 pmod{x^n} ]

    (n=1)时,(G(x))是一个常数(c),那么,(F(x)=c^{-1})

    假设我们已经知道了一个多项式(F_0(x)),使得:

    [F_0(x)G(x)equiv 1 pmod{x^{lceil n/2 ceil}} ]

    那么显然,我们有:

    [F(x)G(x)equiv 1 pmod{x^{lceil n/2 ceil}} ]

    将上述两式相减,我们得到:

    [G(x)(F(x)-F_0(x))equiv 0 pmod{ x^{lceil x/2 ceil}} ]

    由于(G(x))有常数项,所以我们有:

    [F(x)-F_0(x)equiv 0 pmod{x^{lceil n/2 ceil}} ]

    将上式平方得:

    [(F(x)-F_0(x))^2equiv 0 pmod{x^n} ]

    展开:

    [F(x)^2 -2F(x)F_0(x)+F_0(x)^2 equiv 0 pmod{x^n} ]

    因为(F(x))现在是二次的,我们想方法把他变成一次。根据最开始的式子,我们有:

    [F(x)-2F_0(x)+G(x)F_0(x)^2equiv0 pmod{x^n} ]

    所以:

    [F(x)equiv 2F_0(x)-G(x)F_0(x)^2 pmod{x^n} ]

    每次多项式乘法是(O(n log n ))的,所以:

    [T(n)=T(n/2)+O(n log n) = O(n log n) ]

    核心代码:

    inline void P_inv(int* a,int* ans,int n){
    	static int tmp[maxn<<1];
    	if(n==1){
    		ans[0]=ksm(a[0],MOD-2);
    		return;
    	}
    	P_inv(a,ans,(n+1)>>1);
    	int len=1,cnt=0;
    	while(len<(n<<1)) len<<=1,cnt++;
    	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
    	for(int i=0;i<n;i++) tmp[i]=a[i];
    	for(int i=n;i<len;i++) tmp[i]=0;
    	NTT(tmp,len,1); NTT(ans,len,1);
    	for(int i=0;i<len;i++) ans[i]=1ll*(2ll-1ll*tmp[i]*ans[i]%MOD+MOD)%MOD*ans[i]%MOD;
    	NTT(ans,len,-1);
    	for(int i=n;i<len;i++) ans[i]=0;
    }
    

    二、多项式(ln)

    已知一个多项式(G(x)),求一个多项式(F(x)),使得:

    [F(x)equiv ln G(x) pmod{x^n} ]

    等一下,多项式有ln这种东西?

    没错,它有,定义如下:

    [ln(1-F(x))=-sum_{ige 0}frac{F^i(x)}{i} ]

    事实上,你可以把它看作多项式与麦克劳林级数的复合。

    所以根据定义,它的常数项一定为1,否则不存在ln了。

    见到(ln),我们一般都是要求导的:

    [F'(x)equiv frac{G'(x)}{G(x)} pmod{x^n} ]

    然后我们再积分就行了:

    [F(x)equiv int frac{G'(x)}{G(x)}\,dx pmod{x^n} ]

    时间复杂度(O(n log n))

    核心代码:

    inline void P_der(int *a,int *ans,int n){
    	for(int i=1;i<n;i++) ans[i-1]=1ll*i*a[i]%MOD;
    	ans[n-1]=0;
    } 
    inline void P_int(int *a,int *ans,int n){
    	for(int i=1;i<n;i++) ans[i]=1ll*ksm(i,MOD-2)*a[i-1]%MOD;
    	ans[0]=0;
    }
    inline void P_ln(int *a,int *ans,int n){
    	static int inva[maxn<<1],dera[maxn<<1];
    	int len=1,cnt=0;
    	while(len<(n<<1)) len<<=1,cnt++;
    	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
    	for(int i=0;i<len;i++) inva[i]=dera[i]=0;
    	P_inv(a,inva,n);
    	P_der(a,dera,n);
    	NTT(inva,len,1); NTT(dera,len,1);
    	for(int i=0;i<len;i++) inva[i]=(1ll*inva[i]*dera[i])%MOD;
    	NTT(inva,len,-1);
    	P_int(inva,ans,n);
    	for(int i=n;i<len;i++) ans[i]=0;
    }
    

    三、多项式牛顿迭代

    已知一个函数(G(x)),求一个多项式(F(x)),使得:

    [G(F(x))equiv 0 pmod{x^n} ]

    (n=1)时,(G(F(x))equiv 0 pmod{x}),这个单独求。

    假设我们已经知道了一个多项式(F_0(x)),使得:

    [G(F_0(x))equiv 0 pmod{x^{lceil n/2 ceil}} ]

    考虑将(G(F(x)))(F_0(x))处泰勒展开:

    [G(F(x))equiv G(F_0(x))+frac{G'(F_0(x))}{1!}(F(x)-F_0(x))+frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+... pmod{x^{lceil n/2 ceil}} ]

    由于(F(x)-F_0(x)))最低的非(0)次项次数最低为(lceil n/2 ceil),所以我们有:

    [G(F(x))equiv G(F_0(x))+G(F_0(x))(F(x)-F_0(x)) pmod{x^n} ]

    移项,合并同类项,我们有:

    [F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))} pmod{x^n} ]

    四、多项式exp

    已知一个多项式(G(x)),求一个多项式 (F(x)),使得:

    [F(x)=e^{G(x)}pmod{x^n} ]

    多项式有exp?

    和多项式对数函数一样,我们也可以把它看作是麦克劳林级数:

    [e^{F(x)} = sum _{ige 0} frac{F^i(x)}{i!} ]

    这里,多项式的常数项必须为0,否则常数项不收敛。

    两边取自然对数得:

    [ln F(x) equiv G(x) pmod{x^n} ]

    移项得:

    [ln F(x) -G(x) equiv 0 pmod{x^n} ]

    令函数(f(F(x))=ln F(x) -G(x)),对其进行牛顿迭代,我们就有了递推式:

    [F(x)equiv F_0(x) - frac{f(F_0(x))}{f'(F_0(x))} pmod{x^n} ]

    通过计算得到:

    [f'(F(x))=frac{1}{F(x)} ]

    代入,经过计算,我们便得到:

    [F(x)equiv F_0(1+G(x)-ln F_0(x)) pmod{x^n} ]

    时间复杂度

    [T(n)=T(n/2)+O(n log n) = O(n log n) ]

    核心代码:

    inline void P_exp(int *a,int *ans,int n){
    	static int lna[maxn<<1],tmp[maxn<<1];
    	if(n==1){
    		ans[0]=1;
    		return;
    	}
    	P_exp(a,ans,(n+1)>>1);
    	P_ln(ans,lna,n);
    	int len=1,cnt=0;
    	while(len<(n<<1)) len<<=1,cnt++;
    	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
    	for(int i=0;i<n;i++) tmp[i]=a[i];
    	for(int i=n;i<len;i++) tmp[i]=0;
    	NTT(tmp,len,1); NTT(lna,len,1);
    	NTT(ans,len,1);
    	for(int i=0;i<len;i++) ans[i]=1ll*ans[i]*((1+tmp[i]-lna[i]+MOD)%MOD)%MOD;
    	NTT(ans,len,-1);
    	for(int i=n;i<len;i++) ans[i]=0;
    }
    

    (这看似短短的代码调了我两天。。。主要是前面写的ln和inv中有很多数组清空的不彻底。。。)

    五、多项式开根

    已知一个多项式(G(x)),求一个多项式(F(x)),使得:

    [F^2(x)equiv G(x) pmod{x^n} ]

    移项得:

    [F^2(x)-G(x)equiv 0pmod{x^n} ]

    令函数(f(F(x))=F^2(x)-G(x)),对其进行牛顿迭代,我们就有了递推式:

    [F(x)equiv F_0(x) - frac{f(F_0(x))}{f'(F_0(x))} pmod{x^n} ]

    通过计算得到:

    [f'(F(x))=2F(x) ]

    代入,经过计算,我们便得到:

    [F(x)equiv frac{F_0^2(x)+G(x)}{2F_0(x)} pmod{x^n}时间复杂度$O(n log n)$ ]

    将式子变形一下:

    [F(x) equiv frac{F_0(x)}{2} + frac{G(x)}{2F_0(x)} pmod{x^n} ]

    时间复杂度(O(n log n))

    刚刚还看到了一种解法,和大家分享一下:

    两边同取自然对数:

    [ln (F^2(x)) equiv ln (G(x)) pmod{x^n} ]

    所以

    [ln (F(x)) equiv ln(frac{G(x)}{2}) pmod{x^n} ]

    所以

    [F(x)equiv e^{ln(frac{G(x)}{2})} pmod{x^n} ]

    牛顿迭代版:

    inline void P_sqrt(int *a,int *ans,int n){
    	static int tmp[maxn<<1],inva[maxn<<1];
    	if(n==1){
    		ans[0]=1;
    		return;
    	}
    	P_sqrt(a,ans,(n+1)>>1); 
    	int len=1,cnt=0;
    	while(len<(n<<1)) len<<=1,cnt++;
    	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
    	for(int i=0;i<len;i++) inva[i]=0;
    	P_inv(ans,inva,n);
    	for(int i=0;i<n;i++) tmp[i]=a[i];
    	for(int i=n;i<len;i++) tmp[i]=0;
    	NTT(tmp,len,1); NTT(inva,len,1);
    	for(int i=0;i<len;i++) tmp[i]=(1ll*tmp[i]*inva[i])%MOD;
    	NTT(tmp,len,-1);
    	for(int i=0;i<len;i++) ans[i]=1ll*(tmp[i]+ans[i])%MOD*inv2%MOD;
    	for(int i=n;i<len;i++) ans[i]=0;
    }
    

    exp版:

    inline void P_sqrt(int *a,int *ans,int n){
    	static int lna[maxn<<1];
    	P_ln(a,lna,n);
    	for(int i=0;i<n;i++) lna[i]=(1ll*lna[i]*inv2)%MOD;
    	P_exp(lna,ans,n);
    }
    

    然而用exp写的T了3个点。。。还是老老实实写牛迭吧。

    大家有没有发现,在上面牛顿迭代的代码中,我们默认了常数项是1?那么如果常数项不是1呢?

    这时候,我们就要计算二次剩余了。

    又是一个被我咕了的东西

    六、多项式快速幂

    已知一个多项式(G(x)),求一个多项式(F(x)),使得

    [F(x)equiv G^k(x) pmod{x^n} ]

    其中(kin Q)

    因为

    [G(x)=e^{ln G(x)} ]

    所以

    [F(x)equiv e^{k ln G(x)} pmod{x^n} ]

    时间复杂度(O(n log n))

    核心代码:

    inline void P_power(int* a,int *ans,int n){
    	static int lna[maxn];
    	for(int i=0;i<n;i++) lna[i]=0;
    	P_ln(a,lna,n);
    	for(int i=0;i<n;i++) lna[i]=(1ll*lna[i]*k)%MOD;
    	P_exp(lna,ans,n);
    }
    

    在洛谷的模板题里,指数达到了(10^{10^5})级别,那么我们该怎么办呢?其实很简单,我们只需将指数对(MOD)取模就行了。(这个还是比较显然的吧)

    注意,因为ln的缘故,这里常数项也必须是1。

    好了,先将这么多操作,以后我会补充的。

    那我们来一道例题吧:

    P4091 [HEOI2016/TJOI2016]求和

    求:

    [f(n)=sum_{i=0}^nsum_{j=0}^iS(i,j) imes2^j imes(j!) ]

    其中(S(i,j))为第二类斯特林数

    首先我我们呢得知道第二类斯特林数是啥。

    (S(i,j))是用来解决LUB问题的(就是把(i)个不同的球放入(j)个相同的盒子里,且盒子不能为空。有兴趣的可以去了解TwelveFoldway)。

    我们先考虑盒子不同时怎么做。考虑枚举空盒的个数,那么我们就可以大力容斥了,最后因为盒子相同,除以(j!)即可。于是我们得到:

    [S(i,j)=frac{sum_{k=0}^j(-1)^kC_j^k(j-k)^i}{j!} ]

    现在我们开始推式子。

    改变枚举顺序,将(2^j)提取出来,我们得到:

    [f(n)=sum_{j=0}^n2^jsum_{i=0}^nS(i,j) imes(j!) ]

    (注意当(i<j)时,(S(i,j)=0))

    代入第二类斯特林数的计算公式,我们得到:

    [f(n)=sum_{j=0}^n2^jsum_{i=0}^nsum_{k=0}^j(-1)^kC_j^k(j-k)^i ]

    (k)的枚举提前,我们得到:

    [f(n)=sum_{j=0}^n2^jsum_{k=0}^j(-1)^kC_j^ksum_{i=0}^n(j-k)^i ]

    我们可以看到,后面两个sigma相乘,其实可以把他们写成多项式卷积。

    [A[k]=(-1)^kC_j^k \ B[k]=sum_{i=0}^n k^i = frac{k^{n+1}-1}{k-1} ]

    所以我们有:

    [f(n)=sum^n_{j=0}2^j(A*B) ]

    于是我们直接NTT就好了。

    我们再来一道题:P3321 [SDOI2015]序列统计

    给定整数 (x),求所有可以生成出的,且满足数列中所有数的乘积 (mod m)的值等于 (x) 的不同的数列的有多少个。

    既然是计数,那么我们考虑(dp)

    (dp[i][j])表示选了(i)个数,数列的积模(m)的值为(j),那么我们有以下的状态转移方程:

    [dp[2*i][c] = sum_{a*bequiv cpmod{m}} dp[i][a]*dp[i][b] ]

    我们观察这个式子,发现如果是(a+b=c)的话,我们就完全可以卷积嘛!

    那么我们这里就有了一个奇技淫巧:取对数

    我们在高一学必修一的时候就学过:(虽然我还没有高一

    [log _a n*m=log_a n + log _a m ]

    注意我们是模意义下的对数,联想到原根的一些性质(我保证这坑我以后会补),我们只需以原根为底取对数即可

    (A=log_g a \%m,B=log _g b\%m,C=log_g c\%m),于是式子就被我们改写成:

    [dp[2*i][C]=sum_{A+Bequiv Cpmod{varphi(m)}} dp[i][A]*dp[i][B] ]

    好了,我今天就分享到这里了。

    如有不足敬请指正,谢谢!

    参考资料:

    1.pmt巨佬的讲义(%%%pmt巨佬)

    2.OI wiki

    3.Miskcoo's space

  • 相关阅读:
    BZOJ 2527 Meteors 整体二分
    BZOJ 1176: [Balkan2007]Mokia
    DP杂题2
    点分治
    一些图论模板
    一些字符串的题
    斐波那契+线段树
    BZOJ 2957楼房重建
    POJ
    BZOJ 2002 弹飞绵羊
  • 原文地址:https://www.cnblogs.com/ybwowen/p/12198379.html
Copyright © 2011-2022 走看看