多项式牛顿迭代及其应用

前置知识：泰勒展开、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;		
}