zoukankan      html  css  js  c++  java
  • 多项式基础学习笔记(1)

    写这玩意儿的感想就是:一开始常数比你小了一点的,到后来随着各种函数调用来调用去,会变成亿点……

    关于界

    在大多数题目中,我们只对多项式前若干项感兴趣,此时为所有运算设定一个次数上界,即 (mod x^n) 。不难发现有性质:

    [(F(x)mod x^n+((x)mod x^n)mod x^n)mod x^n=(F(x)+G(x))mod x^n\\ (F(x)mod x^n)*(G(x)mod x^n)mod x^n=(F(x)*G(x))mod x^n ]

    这样能避免对无穷幂级数的处理,保证复杂度。

    多项式四则运算

    约定([x^n]F(x))(F[n]) 均表示 (F(x))(n) 次项系数。

    定义多项式的加减为:

    [F(x)+G(x)=sum_{i=0}(F[i]+G[i])x^i\\ F(x)-G(x)=sum_{i=0}(F[i]-G[i])x^i ]

    多项式乘法(加法卷积)为:

    [F(x)*G(x)=sum_{i=0}x^isum_{k=0}^iF[k]G[i-k] ]

    使用 NTT 在 (mathcal{O}(nlog n)) 内完成计算。FFT & NTT

    多项式乘法逆

    如果多项式 (F(x),G(x)) 满足 (F(x)*G(x)=1) ,那么称 (F(x),G(x)) 互为乘法逆。

    这里采用和卷积同样的倍增思想。假设现在已经求出了 (H(x)) 满足 (F(x)H(x)equiv 1(mod x^{lceil n/2 ceil})) ,现在需要通过 (H(x))(G(x)) .

    [ecause F(x)G(x)equiv 1(mod x^{lceil n/2 ceil}) herefore G(x)-H(x)equiv 0(mod x^{lceil n/2 ceil}) ]

    两边平方,得(注意,平方的时候模数也要同步)

    [G(x)^2-2G(x)H(x)+H(x)^2equiv 0(mod x^{n}) ]

    现在用 (F(x)) 消去平方项 (G(x))

    [G(x)-2H(x)+H(x)^2F(x)equiv 0(mod x^n) ]

    移项得

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

    这样就能直接用 NTT 递推计算了。模板题链接

    如果你发现人均 500ms ,就你 1s+ ,不要慌张,开个 O2 试试 /xyx ,我直接 1.22s ( o) 490ms .

    没想到老厌氧选手也有这一天 但是还是不喜 C++17 QAQ 一开准慢。

    注:这里只有当前段代码,也就是针对本题删去了无用部分和封装,但事实上完整代码是……大缝合怪!

    #define Clear(a,n) memset(a,0,sizeof(int)*(n))
    #define Copy(a,b,n) memcpy(a,b,sizeof(int)*(n))
    #define Rev(a,b) reverse(a,b)
    
    int power( int a,int b ) { int res=1; for (;b;b>>=1,a=1ll*a*a%Mod) if (b&1) res=1ll*a*res%Mod; return res; }
    inline void bmod( int &x ) { x-=Mod; x+=x>>31&Mod; }
    
    int rev[M],rt[M],lim,lg,lasn;	//反转数组,原根,上界的值和log
    void Poly_Init( int n ) 
    {
    	if ( n==lasn ) return; lasn=n;
    	for ( lg=0,lim=1; lim<=n; lg++,lim<<=1 );
    	for ( int i=0; i<lim; i++ ) rev[i]=(rev[i>>1]>>1) | ((i&1)<<(lg-1));
    	for ( int i=1,w; i<lim; i<<=1 )
    	{
    		w=power(3,(Mod-1)/(i<<1)); rt[i]=1;
    		for ( int j=1; j<i; j++ ) rt[i+j]=1ll*rt[i+j-1]*w%Mod;
    	}
    }
    void NTT( int *f,int opt )
    {
    	int i,invn,len;
    	for ( i=0; i<lim; i++ ) if ( i>rev[i] ) swap( f[i],f[rev[i]] );
    	for ( i=1; i<lim; i<<=1 )
    		for ( int j=0; j<lim; j+=i<<1 )
    			for ( int k=0,t1,t2; k<i; k++ )
    			{
    				t1=f[j+k]; t2=1ll*rt[i+k]*f[i+j+k]%Mod;
    				bmod(f[j+k]=t1+t2); bmod(f[i+j+k]=t1+Mod-t2);
    			}
    	if ( opt ) return; invn=power(lim,Mod-2);
    	for ( i=0; i<lim; i++ ) f[i]=1ll*f[i]*invn%Mod;
    }
    void Poly_Mul( int n,int m,int *f,int *g,int *res )
    {
    	Poly_Init(n+m); static int tf[M],tg[M];
    	Copy(tf,f,n); Clear(tf+n,lim-n); NTT(tf,1); 
    	Copy(tg,g,m); Clear(tg+m,lim-m); NTT(tg,1);
     	for ( int i=0; i<lim; i++ ) res[i]=1ll*tf[i]*tg[i]%Mod;
     	Rev(res+1,res+lim); NTT(res,0);
    }
    void Poly_Inv( int n,int *f,int *g )
    {	//G(x)equiv 2H(x)-H(x)^2F(x)(mod x^n)
    	static int tf[M];
    	if ( n==1 ) { g[0]=power(f[0],Mod-2); return; }
    	Poly_Inv((n+1)>>1,f,g); Poly_Init(n<<1);
    	Copy(tf,f,n); Clear(tf+n,lim-n); Clear(g+n,lim-n); NTT(tf,1); NTT(g,1);
    	for ( int i=0; i<lim; i++ ) g[i]=1ll*g[i]*(2+Mod-1ll*g[i]*tf[i]%Mod)%Mod;
    	Rev(g+1,g+lim); NTT(g,0); Clear(g+n,lim-n);
    }
    

    多项式除法

    给定 (F(x),G(x)) ,求 (Q(x)) 使得 (F(x)=G(x)Q(x)+R(x))(Q(x)) 次数为 (n-m)(R(x)) 次数小于 (m) .

    (F^R(x)) 表示 (F(x)) 系数翻转之后的结果。设 (F(x)) 最高次项为 (x^{n-1}) ,那么有

    [F^R(x)=x^{n-1}Fleft(dfrac 1x ight) ]

    于是开始推式子

    [Fleft(dfrac 1x ight)=Gleft(dfrac 1x ight)Qleft(dfrac 1x ight)+Rleft(dfrac 1x ight)\\ x^{n-1}Fleft(frac 1x ight)=x^{m-1}Gleft(frac 1x ight)x^{n-m}Qleft(frac 1x ight)+x^{n-m+1}x^{m-2}Rleft(frac 1x ight)\\ F^R(x)=G^R(x)Q^R(x)+x^{n-m+1}R^R(x) ]

    由于 (Q) 的次数是 (n-m) ,所以可以进行模意义下的转化

    [F^R(x)equiv G^R(x)Q^R(x)pmod {x^{n-m+1}} ]

    于是就能求出 (Q^R(x)=dfrac{F^R(x)}{G^R(x)}pmod{x^{n-m+1}}) ,这个用求逆就能完成,求 (R(x)) 就直接 (F-G*Q) 即可。

    void Poly_DivMod( int n,int m,int *f,int *g,int *q,int *r )
    {	//F^R(x)equiv G^R(x)Q^R(x)pmod {x^{n-m+1}}
    	static int tf[M],tg[M],nw[M];
    	Copy(tf,f,n); Copy(tg,g,m); Rev(tf,tf+n); Rev(tg,tg+m);
    	for ( int i=n-m+1; i<(n<<1); i++ ) tg[i]=0;
    	Poly_Inv(n-m+1,tg,nw); Poly_Mul(n,n-m+1,tf,nw,q);
    	for ( int i=n-m+1; i<(n<<1); i++ ) q[i]=0; 
    	Rev(q,q+n-m+1); Poly_Mul(n-m+1,m,q,g,r);
    }
    

    拉格朗日插值

    给定 (n) 个点值,确定多项式 (f(x)) 并求出 (f(k)mod 998244353) .

    拉格朗日他老人家留下的美妙遗产之一:

    [f(k)=sum_{i=0}^ny_iprod_{i eq j}dfrac{k-x_j}{x_i-x_j},给定点值为(x_i,y_i) ]

    详细推导参见 OI-Wiki . 这样直接做就达到了 (mathcal{O}(n^2)) 比消元少了一个 (n) 。在 (x) 连续的时候,可以把 (x_i) 替换成 (i) ,并令 (pre[i]=prodlimits_{j=0}^ik-j)(suf[i]=prodlimits_{j=i}^nk-j)(fac[i]=prodlimits_{j=1}^ij) ,那么有

    [f(k)=sumlimits_{i=0}^{n}y_ifrac{pre_{i-1}*suf_{i+1}}{fac_{i}*fac_{n-i}}*(-1)^{n-i} ]

    优化到了 (mathcal{O}(n)) .

    重心拉格朗日插值

    回到这个式子:

    [f(k)=sum_{i=0}^ny_iprod_{i eq j}dfrac{k-x_j}{x_i-x_j} ]

    (g=prodlimits_{i=0}^nk-x[i]) ,原式可以变为

    [f(k)=gsumlimits_{i=0}^{n}frac{y_i}{k-x_i}prodlimits_{i e j}frac{1}{x_i-x_j} ]

    (t[i]=dfrac{y[i]}{prodlimits_{i eq j}x_i-x_j}) ,有

    [f(k)=gsumlimits_{i=0}^{n}dfrac{t[i]}{k-x[i]} ]

    每次修改或插入直接修改 (t[i]) 即可。注意不能处理 (k=x[i]) 的情况。模板题链接

    暴力 (mathcal{O}(n^2)) 的插值 82ms 跑进第一页也是没想到啊……

    //Author: RingweEH
    const int N=2e3+10,Mod=998244353;
    int n,k,x[N],y[N];
    int power( int a,int b ) { int res=1; for (;b;b>>=1,a=1ll*a*a%Mod) if (b&1) res=1ll*a*res%Mod; return res; }
     
    int main()
    {
    	n=read(); k=read();
    	for ( int i=1; i<=n; i++ ) x[i]=read(),y[i]=read();
    
    	ll ans=0;
    	for ( int i=1; i<=n; i++ )
    	{
    		ll a=1,b=1;
    		for ( int j=1; j<=n; j++ )
    		{
    			if ( i==j ) continue;
    			a=a*(k-x[j]+Mod)%Mod; b=b*(x[i]-x[j]+Mod)%Mod;
    		}
    		ans=(ans+a*power(b,Mod-2)%Mod*y[i]%Mod)%Mod;
    	}
    
    	printf( "%lld
    ",ans );
    
    	return 0;
    }
    

    泰勒展开 & 牛顿迭代

    原理:用一个 (n) 次多项式高度拟合一个函数,新的函数某一点 (x_0)(n) 阶导数和原函数的 (n) 阶导数相同,并且都经过 ((x_0,f(x_0))) .

    具体可以参考《具体数学》5.3技巧2:高阶差分 关于牛顿级数和泰勒级数的阐述。

    函数 (F(x))(x_0) 位置泰勒展开,令 (F^n(x)) 表示 (F(x))(n) 阶导数,有

    [G(x_0+x)=dfrac{F(x_0)}{0!}x^0+dfrac{F^1(x_0)}{1!}x+dfrac{F^2(x_0)}{2!}x^2+cdots ]

    得到这个式子的推导可以类比牛顿级数的推导,先得到 (F^n(x_0)) 为当前项的系数且在更大的时候均为 (0) (根据求导易证),然后再根据求导得到的系数部分得到下面的分子。

    往往在精度满足的时候就停止,所以大部分都是有误差的。


    考虑多项式求逆的过程,假设已经有 (F(x)*H(x)equiv 1(mod x^{lceil n/2 ceil})) ,函数 (G(F(x))equiv 0(mod x^n)) .

    (这里显然有 (G(H(x))equiv 0pmod x^{lceil n/2 ceil})

    (G(F(x)))(H(x)) 泰勒展开

    [G(F(x))=dfrac{G^0(H(x))}{0!}(F(x)-H(x))^0+dfrac{G^1(H(x))}{1!}(F(x)-H(x))+dfrac{G''(H(x))}{2!}(F(x)-H(x))^2+cdots\ ]

    从第三项开始,由于 (F(x))(H(x))(leftlceildfrac{n}{2} ight ceil) 项相同,所以模意义下为 (0) .

    所以有

    [G(F(x)) equiv 0equiv G(H(x))+G^1(H(x))(F(x)-H(x))pmod{x^n}\ ]

    这样就能求得 (F(x))

    [F(x)=H(x)-dfrac{G(H(x))}{G'(H(x))}\ ]

    多项式 ln

    如果不太明白为什么多项式会有自然对数这种东西,可以看 这里 (但是貌似还是需要一些高级的东西……先感性理解吧)

    (B(x)=ln(A(x))) 两边求导

    [B'(x)=dfrac{A'(x)}{A(x)} ]

    然后对 (A) 求逆,求导,卷起来,然后积分,就得到了 (B) . 题目链接

    //Author: RingweEH
    void Itg( int n,int *f,int *g ) //Intergral
    {
    	for ( int i=1; i<=n; i++ ) g[i]=1ll*f[i-1]*Math::inv[i]%Mod; g[0]=0;
    }
    void Drv( int n,int *f,int *g ) //Derivative
    {
    	for ( int i=0; i<n-1; i++ ) g[i]=1ll*(i+1)*f[i+1]%Mod; g[n-1]=0;
    }
    void Poly_ln( int n,int *f,int *g )
    {
    	static int tf[M],tg[M];
    	Drv(n,f,tf); Clear(tg,n); Poly_Inv(n,f,tg); 
    	Poly_Mul(n,n,tf,tg,tf);  Itg(n,tf,g);
    }
    

    多项式 exp

    已知 (A(x)) ,求 (F(x)=e^{A(x)}) ,保证 (A(0)=0) .

    对两边求导得 (ln(F(x))=A(x)) . 令 (G(F(x))=ln(F(x))-A(x)) ,那么有 (G'(F(x))=dfrac{1}{F(x)}) .

    由于 (A(x)) 给出,在 (x) 固定时是固定的,可以看做常量。

    然后对其牛顿迭代

    [F(x)=H(x)-dfrac{G(H(x))}{G'(H(x))}=H(x)(1-ln(H(x))+A(x)) ]

    同样递归计算 (H(x)) 即可,边界为 (G(0)=1) . 这样的复杂度是 (mathcal{O}(nlog n)) 的。


    如果 (n) 很小且模数不一定是 NTT 模数,那么可以换一种做法。

    [F(x)=exp(A(x))=>ln(F(x))=A(x)=>dfrac{F'(x)}{F(x)}=A'(x)=>F'(x)=F(x)A'(x) ]

    所以有

    [F'[n]=sum_{i=0}^nA'[i]F[n-i] ]

    积分回去

    [(n+1)F[n+1]sum_{i=0}^n(i+1)A[i+1]F[n-i] ]

    模板题链接

    void Poly_Exp( int n,int *f,int *g )
    {	//G(x)=H(x)(1-ln(H(x))+A(x))
    	if ( n==1 ) { g[0]=1; return; } static int tf[M];
    	Poly_Exp( (n+1)>>1,f,g ); Clear(tf,n); Poly_ln(n,g,tf); tf[0]--;
    	for ( int i=0; i<n; i++ ) bmod( tf[i]=f[i]+Mod-tf[i] );
    	Poly_Mul(n,n,g,tf,g); Clear(g+n,lim-n);
    }
    

    多项式开根

    给定 (A(x)) ,求 (F^2(x)equiv A(x)pmod x^n) ,保证 (a_0=1) .

    (G(F(x))=F^2(x)-A(x)) ,那么 (G'(F(x))=2F(x)) . 对其牛顿迭代

    [egin{aligned} F(x) &=H(x)-dfrac{G(H(x))}{G'(H(x))}\\ &=H(x)-dfrac{H^2(x)-A(x)}{2H(x)}=dfrac{H(x)^2+A(x)}{2H(x)} end{aligned} ]

    同样递归求解即可。(G(0)=1) . 模板题链接

    void Poly_Sqrt( int powe,int *f,int *g )
    {//F(x)=dfrac{H(x)^2+A(x)}{2H(x)}
    	if ( powe==1 ) { g[0]=1; return; }
    	static int tf[M]; Poly_Sqrt((powe+1)>>1,f,g);
    	Clear(tf,powe); Poly_Inv(powe,g,tf); Poly_Mul(powe,powe,tf,f,tf);
    	int inv2=power(2,Mod-2); for ( int i=0; i<powe; i++ ) g[i]=1ll*(g[i]+tf[i])*inv2%Mod;
    }
    

    多项式快速幂

    经典用法:(ln+exp) 把幂次转化成四则运算。

    [B(x)equiv A^k(x)pmod{x^n}=>ln(B(x))=kcdot ln(A(x)) ]

    然后 (exp(ln(B(x)))) 即可。

    模板题链接

    注意:本题 (kleq 10^{10^5}) ,所以要用到:

    任意质数 (p) 和任意多项式 (F(x))(F^p(x)equiv 1pmod{x^n}) (还得要求 (a_0=1) ),具体可以看具体数学和 这个帖子

    总之就是读入的幂次要 (mod p) 啦。

    void Poly_Power( int n,int k,int *f,int *g )
    {
    	int tf[M]; Clear(tf,n); Poly_ln(n,f,tf);
    	for ( int i=0; i<n; i++ ) tf[i]=1ll*tf[i]*k%Mod;
    	Poly_Exp(n,tf,g);
    }
    

    一些技巧

    差卷积

    [C_k=sum_{i=k}^nA_kB_{i-k} ]

    (B_1[n-i]=B[i]) ,那么有

    [C[k]=sum_{i=0}^kA[i]B_1[k-i]=sum_{i=0}^kA[i]B[n+i-k] ]

    所以 (C[i]=(A*B_1)_{n+i}) .

    多项式平移

    (F(x)) 得到 (F(x+c)) 的系数。设 (F(x)=sum_{i=0}^na_ix^i) .

    [egin{aligned} F(x+c)&=sum_{i=0}^{n}a_i(x+c)^i=sum_{i=0}^{n}a_isum_{j=0}^{i}inom{i}{j}x^{j}c^{i-j}\\ &=sum_{j=0}^{n}x^jsum_{i=j}^{n}inom{i}{j}a_ic^{i-j}=sum_{j=0}^{n}x^jdfrac{1}{j!}sum_{i=j}^{n}i!a_idfrac{c^{i-j}}{(i-j)!} end{aligned} ]

    差卷积即可。

    Hints

    • 如果你想对 Poly_Mul 卡常的话,可以试试小范围暴力大范围卷积。
    • 不妨对手写取模加个 inline ,效果会很不错。
    • 这里 F3 “卡常后的倍增”会得到一个 Poly_Inv 的小优化。
    • 似乎 i+j+ki|j|k 快。
    • memset 应该比 fill 和循环快
    • 指针很强 虽然我不会 ,如果会指令集那我只能 Orz
    • 一个避免大量 memset 的好方法是每次在封装函数里面开 static ,而且能有效避免传指针修改外界变量的情况

    习题

    • 不要再忘记 Math::Init(n) 了!
    • 时刻记住你的卷积开闭区间情况。

    [E_i=sum_{j=1}^{i-1}dfrac{q_j}{(j-i)^2}-sum_{j=i+1}^{n}dfrac{q_j}{(i-j)^2} ]

    前面那个式子可以卷积,后面那个式子反一反同样卷积,就没了。不会吧不会吧你都看过多项式除法了不会想不到 reverse 吧

    礼物

    设当前调整值为 (c)(x_i=a_i-b_i) ,先来考虑一个简化的问题:已知 (b) 序列,如何确定 (c) 使 (sum(x_i+c)^2) 最小?

    [sum(x_i+c)^2=sum x_i^2+2csum x_i+n imes c^2 ]

    整个式子可以看做是一个关于 (c) 的二次函数,那么直接枚举 (c) 就能 (mathcal{O}(m)) 求解了。

    现在来考虑不固定的情况。注意到 (sum x_i) 是不会变的,那么式子依然可以拆分成两部分,一部分是最小化 (nc^2+2cX) (依然可以直接做),一部分是最大化 (sum a_ib_i) 。将 (a) 反向,构造卷积 (sum a_ib_{n-i}) ,由于 (b) 是环形,所以倍长之后取后 (n) 个中的最大值即可。复杂度是 (mathcal{O}(nlog n+m)) .

    差分与前缀和

    先考虑前缀和,发现其实就是和一个系数全 (1) 的多项式做一次卷积。那么 (k) 次前缀和,根据卷积的结合律,乘上去的多项式就是 (dfrac{1}{(1-x)^k})

    类似思考差分,一维差分要卷积的式子是 (1-x) ,那么 (k) 次就是 ((1-x)^k) ,直接多项式快速幂。

    这是比较暴力的做法……所以不开 O2 就全 TLE ,否则就会 AC

    开拓者的卓识

    再次说明了“ 算贡献 ”的思想是真的很重要。

    我们考虑对于一个原始序列的 (a_i) ,会对哪些 ([k,l,r]) 产生贡献。注意,这里不仅仅要选出最后的 (l,r) ,还要考虑之前的所有的 (k-1) 个区间,所以相当于求出满足 (1leq l_kleq l_{k-1}leq l_{k-2}leq cdots leq ileq cdots leq r_{k-1}leq r_k)({(l_1,r_1),(l_2,r_2),cdots ,(l_k,r_k)})集合个数,由于左右端点独立,那么可以直接插板得到

    [a_i imesinom{i+k-2}{k-1}inom{r-i+k-1}{k-1} ]

    注意到组合数的下指标都是 (k-1) ,所以可以一遍递推得到,令 (displaystyle b_i=inom{k-1+i}{k-1}) ,答案就是

    [sum a_ib_{i-1}b_{r-i} ]

    显然可以卷积。

    学习材料

    tommy0221多项式笔记(一)

    zhouzhendong多项式入门教程

  • 相关阅读:
    axios、ajax、fetch三者的区别
    React与Vue的相同与不同点
    react-redux
    redux【react】
    react高阶组件
    基于WebGL无插件虚拟场景漫游关键技术(完整版)ThingJS
    基于WebGL的三维交通监控可视化技术应用(实践版) ThingJS
    地下管线监控系统中互联网WebGL三维可视化构建技术 ThingJS
    基于WebGL实现智慧校园的全景漫游技术研究 三维可视化
    基于WebGL的3D可视化告警系统关键技术解析 ThingJS
  • 原文地址:https://www.cnblogs.com/UntitledCpp/p/Poly_AllInOne_1.html
Copyright © 2011-2022 走看看