多项式牛顿迭代及其应用

前置知识：泰勒展开、NTT

多项式牛顿迭代

已知函数$ G(x)$，求满足函数方程：

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

的$ F(x)$。

依然考虑倍增法：

假设我们已经求出\(F_0(x)\)满足

\[G(F_0(x))\equiv 0\pmod {x^{\frac n2}} \]

在\(F_0(x)\)处泰勒展开：

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

不了解泰勒展开的可以看看这个问题下的第一篇回答，前半段可以直接跳过。。。

回归正题，因为显然有

\[G(F(x))\equiv 0\pmod {x^{\frac n2}}\\ \therefore F(x)-F_0(x)\equiv0\pmod {x^{\frac n2}} \]

所以如果\(i\ge2\),那么

\[(F(x)-F_0(x))^i\equiv0\pmod {x^n} \]

因此，之前泰勒展开的式子在\(\pmod {x^n}\)下只需要考虑前2项，等价于：

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

最终我们得到：

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

这就是牛顿迭代的式子了，有了它，我们就可以简易的解决许多问题，比如下面这些：

---

多项式指数函数（多项式 \(\exp\))

我们要求\(B(x)\)满足

\[B(x)\equiv e^{A(x)}\pmod {x^n} \]

也就是说

\[\ln B(x)\equiv A(x)\pmod {x^n} \]

相当于已知函数\(G(x)=\ln(x)+A(x)\)，

求满足函数方程：

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

的$ F(x)$。

于是直接套公式：

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

因为\(A(x)\)是已知量，可以看作常数，而\(\ln'(x)=\frac 1x\)

所以\(G'(x)=\frac 1x\)

于是

\[F(x)\equiv F_0(x)-\frac{\ln(F_0(x))-A(x)}{\frac{1}{F_0(x)}}\pmod {x^n}\\ \equiv F_0(x)-F_0(x)(\ln(F_0(x))-A(x))\pmod {x^n} \]

在此之前，我们已经知道了多项式\(\ln\)的求法，于是我们就可以递归下去完成此题了。

当\(A(x)\)只有一项时，因为保证了\(a_0=0\)，所以\(F(x)=1\)（如果不保证，那么\(e\)的其他次方都是无理数，无法取模）

根据主定理，多项式\(\exp\)的复杂度是\(\mathcal O(nlog(n))\)，但是它的常数较大，需要注意。

这里给出\(\exp\)部分的代码，其余部分请参考博主之前的博客，或者看博主的多项式模板集合

inline poly getexp(int n,poly f){
      if(n==1){poly g;g[0]=1;return g;}
      poly g=getexp(n+1>>1,f),B=getln(n,g);
      for(int i=0;i<n;++i) B[i]=dec((i==0),dec(B[i],f[i]));
      return mul(n,n,g,B);
}

---

多项式快速幂

这个和牛顿迭代没有什么关系，但和$ \exp$有关，就顺便讲了。

这一次我们需要求：

\[B(x)\equiv A^k(x)\pmod {x^n} \]

看到幂，直接两边取对数：

\[\ln B(x)\equiv k\ln A(x)\pmod {x^n} \]

于是：

\[B(x)=e^{k\ln A(x)} \]

所以对\(A\)求出\(\ln\)，乘\(k\)后再\(\exp\)即可：

inline poly Pow(int n,int k,poly f){
	poly g=getln(n,f);
	for(int i=0;i<n;++i) g[i]=1ll*g[i]*k%mod;
	return getexp(n,g);
}//求n次多项式f的k次方(f0=1)

这道题目中\(k\)很大，但是根据多项式定理，对于任意质数\(p\)和多项式\(F(x)\)，\(F^p(x)\equiv 1\pmod {x^n}\)

于是读入的时候直接对\(mod\)取模即可：

inline int readpow(){
	scanf("%s",s+1);
	long long x=0,f=1;int len=strlen(s+1);
	for(int i=1;i<=len;++i){char ch=s[i];x=((x<<3)+(x<<1)+(ch^48));if(x>=n) flag=1;x%=mod;}
	return f==-1?mod-x:x;
}

---

多项式幂函数 (加强版)

这道题目和上一道唯一的区别就在于是否保证\(a_0=1\)，这有什么影响呢？

考虑之前的做法，我们会先对\(A\)求出\(\ln\)，可是如果不保证\(a_0=1\)那么\(A\)不可能有对数函数，因为\(e\)的幂中除了\(e^0\)都是无法取模的。

考虑将\(a_0\)强行变成\(1\)，事实上因为这道题是求幂函数，所以我们有：

\[A^k(x)=(\frac{A(x)}{a_0})^ka_0^k \]

所以我们可以直接让\(A\)的每一项系数都除以\(a_0\)，计算完后再将\(a_0^k\)乘回来。

你以为这样就完了吗？\(a_0\)可能为\(0\)：

对于这个问题，我们需要找到\(A\)中最低的系数非零的项，记为\(a_tx^t\)，然后整个多项式除以\(a_tx^t\)，于是我们成功的让\(a_0\)变成了\(1\)，计算之后我们需要将答案右移\(t^k\)位，再乘\(a_t^k\)。

注意以下几点：

• \(t^k\)可能大于\(n\)，需要特判，同时\(k\)取模后可能会变得很小影响这一特判，可以再记一个标记表示读入\(k\)时是否出现过\(\ge n\)的时候。
• 之前我们讲过：多项式求\(k\)次幂时，\(k\)可以直接对\(p\)取模，但是在计算\(t^k\)时，这里的\(k\)我们应该对\(p-1\)取模。

代码如下：

inline poly Pow(int n,int k,poly f){
	poly g=getln(n,f);
	for(int i=0;i<n;++i) g[i]=1ll*g[i]*k%mod;
	return getexp(n,g);
}
namespace fastpow{
	int flag,kk,ne;
	poly ret;
	char s[N];
	inline int readpow1(){
		scanf("%s",s+1);
		long long x=0,f=1;int len=strlen(s+1);
		for(int i=1;i<=len;++i){char ch=s[i];x=((x<<3)+(x<<1)+(ch^48));if(x>=n) flag=1;while(x>=mod) x-=mod;}
		return f==-1?mod-x:x;
	}
	inline int readpow2(){
		long long x=0,f=1;int m=mod-1,len=strlen(s+1);
		for(int i=1;i<=len;++i){char ch=s[i];x=((x<<3)+(x<<1)+(ch^48));while(x>=m) x-=m;}
		return f==-1?m-x:x;
	}
	inline int init(int &n,poly& f){
		int now=0;
		while(f[now]==0) ++now;
		int iv=ksm(f[now],mod-2),ans=ksm(f[now],kk);
		for(int i=now;i<n;++i) f[i-now]=1ll*f[i]*iv%mod;
		long long x=1ll*now*k;
		if(x>=n||(flag&&now!=0)){
			for(int i=0;i<n;++i) ret[i]=0;
			return -1;
		} 
		for(int i=0;i<x;++i) ret[i]=0;
		ne=x;
		n-=now;
		return ans;
	}//init:将f的最低非零系数变为1 
	inline poly notone_pow(int n,int k,poly f){
		ne=n;kk=readpow2();int t=n;int p=init(n,f);
		if(p==-1) return ret;
		f=Poly::Pow(n,k,f);
		for(int i=ne;i<t;++i) ret[i]=1ll*f[i-ne]*p%mod;
		return ret;
	}
}

---

多项式开根

我们要求\(F(x)\)满足：

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

令

\[G(x)=x^2-A(x) \]

于是我们又转化为了：

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

当然也是直接上牛顿迭代：

\[F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\\ \equiv F_0(x)-\frac{F_0^2(x)-A(x)}{2F_0(x)}\pmod {x^n}\\ \equiv \frac 12(F_0(x)+\frac{A(x)}{F_0(x)})\pmod {x^n} \]

于是递归处理即可，注意递归的边界为：\(A(x)\)仅有一项时，因为题目保证了\(a_0=1\)，所以\(F(x)=1\)

inline poly Sqrt(int n,poly f){
	if(n==1){poly g;g[0]=1;return g;}
	poly g=Sqrt((n+1)>>1,f);
	poly p=getinv(n,g);
	poly h=mul(n,n,p,f);
	for(int i=0;i<n;++i) h[i]=1ll*((mod+1)/2)*add(h[i],g[i])%mod;
	return h;		
}