多项式基础学习笔记(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多项式入门教程

原文地址:https://www.cnblogs.com/UntitledCpp/p/Poly_AllInOne_1.html