【Learning】 多项式全家桶（加减乘除逆元ln exp）

欢迎订购KFC多项式全家桶

　　　

约定的记号

　　
　　 对于一个多项式(A(x))，若其最高次系数不为零的项是(x^k)，则该多项式的次数为(k).
　　
​　　 记为(deg(A)=k).
　　
​　　 对于(xin(k,+ infty))，称(x)都为(A(x))的次数界. 但一般地，我们都使用(k+1)作为(A(x))的次数界。
　　
　　
　　

取模意义下的多项式运算

　　
​　　 这里要消除一下疑惑，两个次数界为(n)的多项式在取模(x^n)意义下相乘，虽然NTT完了之后我们会得到一个次数界为(2n-1)的多项式，但是(x^n....x^{2n-1})的系数是没有意义的。答案多项式就是(x^0...x^{n-1})的系数描述的多项式。
　　
​　　 求分式同理，(frac {F(x)}{G(x)}pmod {x^n})的答案多项式，就是(F(x)*G^{-1}(x))的前(n)个系数描述的多项式。
　　
　　
　　

多项式加法减法乘法

　　
​　　 略去，加减为(mathcal O(n))系数运算，乘法用NTT在(mathcal O(nlg n))内计算。
　　
　　
　　

多项式求逆

　　
​　　 已知次数界为(n)的多项式(A(x))，求其在模(x^n)意义下的逆多项式(B(x))，使得

[egin{equation} label{eqn1} A(x)B(x)equiv1 pmod {x^n} end{equation} ]

　　​ 其中，(B(x))的次数小于等于(A(x))的次数
　　
　　
　　
​　　 采用倍增思路向上倍增求解
　　
​　　 当(n=1)时，(A(x))与(B(x))仅有一个常数项，且(B(x))的常数项为(A(x))常数项的逆元。多项式有无逆元也取决于这个常数是否有逆元
　　
​　　 假设已经求得(A(x))在模(x^{lceilfrac n 2 ceil})意义下的逆元(B'(x))

[egin{equation} A(x)B'(x)=1 pmod{x^{lceilfrac n 2 ceil}} end{equation} ]

　　​ 而((1))放在模(x^{lceilfrac n 2 ceil})意义下同样成立，有

[egin{equation} label{eqn2} A(x)B(x)equiv1 pmod {x^{lceilfrac n 2 ceil}} end{equation} ]

　　​ 将((3)-(2))得到

[egin{equation} B(x)-B'(x)equiv 0 pmod {x^{lceilfrac n 2 ceil}} end{equation} ]

​　　 平方得到

[egin{equation} B^2(x)-2B(x)B'(x)+B'^2(x)equiv 0 pmod {x^n} end{equation} ]

　　​ 模数同时平方的原因是：原本多项式模(x^{lceilfrac n 2 ceil})为0，说明(0)~(lceilfrac n 2 ceil-1)项的系数为0，平方后由于系数相乘，这些0系数会导致0~n-1项系数为0，也即模(x^n)为0
　　
　　​ 两边同乘(A(x))，消去(B(x))并移项，得到

[B(x) equiv 2B'(x) - A(x)B'^2(x) pmod {x^n} ]

　　​ 至此可以递归倍增求解，伪代码如下，忽略清零、取模操作：
　　

void polyInv(int *a,int *b,int n){ // a是要求逆元的多项式，b是在模x^n意义下的a的逆元
	if(n==1) 令b为一个次数为1,常数项为a[0]的逆元的多项式，返回;
  	polyInv(a,b,(n+1)>>1);
  	ntt_init(dega+degb+degb); //ntt长度
  	static int A[]=a,B[]=b; //临时数组，防止影响到传入的指针
  	ntt(A);
  	ntt(B);
  	for(int i=0;i<nttlen;i++)
    	A[i]=2*B[i]-A[i]*B[i]*B[i]; //点值计算
  	ntt(A,-1);
  	b=A;
}

　　
　　完整代码如下：
　　

void polyInv(int *a,int *b,int n){
	if(n==1){
		b[0]=fmi(a[0],MOD-2);
		return;
	}
	int m=(n+1)>>1;
	polyInv(a,b,m);
	NTT::init(2*m+n);
	static int ta[S],tb[S];
	for(int i=0;i<n;i++) ta[i]=a[i],tb[i]=b[i];
	for(int i=n;i<NTT::n;i++) ta[i]=tb[i]=0;//NTT前对后半也要清空
	NTT::ntt(ta,0);
	NTT::ntt(tb,0);
	for(int i=0;i<NTT::n;i++) 
		ta[i]=(2LL*tb[i]-1LL*ta[i]*tb[i]%MOD*tb[i]%MOD)%MOD;
	NTT::ntt(ta,1);
	for(int i=0;i<n;i++) b[i]=ta[i];
}

　　
​　　 Tips：调用整个递归之前需要清空(b)数组，否则(b)在赋值给(B)准备NTT时可能有处于([(n+1)>>1,n))中的杂乱数字，将影响NTT。或者每次递归完成后都清空一次也可以。总之记得清空无用数据。
　　
​　　 时间复杂度(O(n log n))，不过我不知道怎么证。
　　
　　
　　
　　
　　

多项式除法及取模

　　
　　​ 给定多项式(A(x))和多项式(B(x))，求两个多项式(D(x))与(R(x))，使得

[egin{equation} label{1} A(x)=D(x)B(x)+R(x) end{equation} ]

　　​ 其中，(deg(A)>=deg(B)) , (deg(D)leq deg(A)-deg(B))，(deg(R)<deg(B))
　　
　　
　　
​　　 除法用奇妙变化解决.
　　
​　　 引入一个操作：翻转操作. 对于一个次数为(n)的多项式(A(x))，定义

[A^R(x)=x^nA(frac 1 x) ]

​　　 会发现(A^R(x))的系数相对于(A(x))完全(reverse)了一下
　　
​　　 将((1))的(x)换成(frac 1 x)，并两边同乘(x^n)，记(n=deg(A)), (m=deg(B)), (deg(D)=n-m), (deg(R)=m-1).

[egin{aligned} x^nA(frac1x)&=x^{n-m}D(frac1x)x^mB(frac1x)+x^{n-m+1}x^{m-1}R(frac1x)\ A^R(x)&=D^R(x)B^R(x)+x^{n-m+1}R^R(x) end{aligned} ]

​　　 我们发现(D^R(x))的次数仍然等于(n-m)，如果把上式放在模(x^{n-m+1})意义下，我们会发现(D^R(x))不会受到任何影响，而(x^{n-m+1}R^R(x))被完全模掉了！
　　
​　　 于是

[egin{aligned} A^R(x)&equiv D^R(x)B^R(x)pmod{x^{n-m+1}}\ D^R(x)&equiv A^R(x){B^R}^{-1}(x)pmod{x^{n-m+1}} end{aligned} ]

​　　 直接对(A)系数翻转，(B)系数翻转，求(B)的逆，(A)和(B)一乘，把结果的系数再次翻转，就是(D)了！
　　

void polyDiv(int *a,int *b,int *d){// a和b对应上文的A和B，d是上文的D，是答案
  	if(deg(a)<deg(b)) 令d为一个0多项式，返回;
	reverse(a);
  	reverse(b);
  	static int invb[];
  	polyInv(b,invb,deg(b)+1); //b次数界是deg(b)+1
  	d=a*invb;
  	reverse(d);
}

　　
　　​ 如果你还要求(R(x))，带回最初的式子直接算就好，多项式乘法用*直接代替了
　　

void polyMod(int *a,int *b,int *r){
  	if(deg(a)<deg(b)) 令r为a，返回;
  	static int d[];
  	polyDiv(a,b,d);
  	r=a-d*b;
}

　　
　　
　　
　　
　　

多项式求(ln)

　　
​　　 已知次数界为(n)的多项式(f(x))，并知关系(f(x)=e^{g(x)})。求未知多项式(g(x)=ln f(x))
　　
​　　 两边求导，得到

[g'(x)=frac {f'(x)}{f(x)} ]

　　​ 所以

[g(x)=intfrac{f'(x)}{f(x)} ]

​　　 先对(f)求导（(mathcal O(n))），再与(f)的逆多项式相乘（都是(mathcal O(n lg n))），对结果再积分（(mathcal O(n))），得到(g(x))
　　

void polyDeri(int *a,int *b,int n){//计算一个次数界为n的多项式a的导数，放入b
	for(int i=0;i<n-1;i++) b[i]=1LL*(i+1)*a[i+1]%MOD;
	b[n-1]=0;
}
void polyInte(int *a,int *b,int n){//计算一个次数界为n的多项式a的积分，放入b
	for(int i=n-1;i>0;i--) 
		b[i]=1LL*a[i-1]*inv[i]%MOD;
	b[0]=0;	
}
void polyLn(int *f,int *g,int n){//整体计算于模x^n意义下进行
	static int t1[S],t2[S];
	polyDeri(f,t1,n);//求导，放入t1
	memset(t2,0,sizeof t2);
	polyInv(f,t2,n);//求逆，放入t2
	NTT::init(n*2-1);//接下来是t1与t2相乘
	for(int i=n;i<NTT::n;i++) t1[i]=t2[i]=0;
	NTT::ntt(t1,0);
	NTT::ntt(t2,0);
	for(int i=0;i<NTT::n;i++) t1[i]=1LL*t1[i]*t2[i]%MOD;
	NTT::ntt(t1,1);//此时t1只有0~n-1项的系数有意义，即为上述分式
	polyInte(t1,g,n);//积分，放入g
}

　　
　　
　　
　　
　　

多项式牛顿迭代法

　　
​　　 已知多项式(g(x))，求(f(x))使得(g(f(x))equiv0pmod{x^n})
　　
​ 做法：倍增模数(x^n)。
　　
​　　 当(n=1)时，(f(x)=1)（令常数项为1）。
　　
​　　 假设已知(f_0(x))使得

[g(f_0(x))equiv0pmod{x^{lceil frac n2 ceil}} ]

​　　 则构造(f(x))

[f(x)equiv f_0(x)-frac{g(f_0(x))}{g'(f_0(x))}pmod{x^n} ]

​　　 可满足(g(f(x))equiv0pmod{x^n})
　　
　　
　　
　　
　　

多项式求exp

　　
​　　 已知次数界为(n)的多项式(f(x))，并知关系(g(x)equiv e^{f(x)}pmod{x^n})。求未知多项式(g(x))
　　
​　　 前置技能：多项式牛顿迭代法。
　　
​　　 先来转化一下式子：

[egin{aligned} g(x)-e^{f(x)}&equiv 0pmod{x^n}\ ln g(x)-f(x)&equiv 0 pmod{x^n} end{aligned} ]

​　　 其实就是要求解多项式函数

[H(g(x))=ln g(x)-f(x)pmod{x^n} ]

​　　 的零点多项式。其中(g(x))是变量，(f(x))是常量多项式。
　　
​　　 套用多项式牛顿迭代，假设已知多项式(g_0(x))满足

[ln g_0(x)-f(x)equiv 0pmod{x^{lceil frac n2 ceil}} ]

　　​ 则构造(g(x))

[egin{aligned} g(x)&equiv g_0(x)-frac{ln g_0(x)-f(x)}{frac 1 {g_0(x)}}\ &equiv g_0(x)*(1-ln g_0(x)+f(x))pmod{x^n} end{aligned} ]

　　​ 可满足(ln g(x)-f(x)equiv 0pmod{x^{n}})。
　　
　　​ 向上倍增即可。
　　

void polyExp(int *f,int *g,int n){
	if(n==1){
		g[0]=1;
		return;
	}
	static int t[S];
	polyExp(f,g,(n+1)>>1);//求得g0，放在g中
	polyLn(g,t,n);//求g0的ln，放在t中
	for(int i=0;i<n;i++) t[i]=(f[i]-t[i])%MOD;//构造后面一个括号的多项式，还是放在t中
	t[0]++;
	NTT::init(n*2-1);//g0和t相乘
	for(int i=n;i<NTT::n;i++) t[i]=g[i]=0;
	NTT::ntt(t,0);
	NTT::ntt(g,0);
	for(int i=0;i<NTT::n;i++) g[i]=1LL*t[i]*g[i]%MOD;
	NTT::ntt(g,1);
	for(int i=n;i<NTT::n;i++) g[i]=0;//清空不必要的东西
}

　　
​　　 同求逆，最好在整体调用之前清空(g)数组，以免不必要的错误。
　　　　
​