【学习笔记】牛顿迭代

Taylor 展开

对于一个函数\(f(x),\)如果我们知道它在\(x_0\)处的各阶导数，那么：

\[f(x)=\sum_{i=0}^n \frac{f^{(i)}(x_0)(x-x0)^i}{i!} \]

即 我们在\(x_0\)处逼近了\(f(x).\)

牛顿迭代

考虑求：

\[G(F(x))\equiv 0(\bmod x^n) \]

对于\(n=1\)特殊求出来

考虑已经解决了：

\[G(F_0(x))\equiv 0(\bmod x^{\left\lceil \frac{n}{2}\right \rceil} ) \]

考虑如何拓展到\(x^n.\)

在这里泰勒展开一下：

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

注意到当\(i\ge 2\)时 \((F(x)-F_0(x))^i\)的最低非\(0\)次项的次数是严格大于\(2\left\lceil\frac{n}{2}\right\rceil,\)所以：

\[G(F(x))\equiv G(F_0(x))+(F(x)-F_0(x))G'(F_0(x))(\bmod x^n) \]

注意到由题设得：

\[G(F(x))\equiv 0(\bmod x^n) \]

所以:

\[G(F_0(x))+(F(x)-F_0(x))G'(F_0(x))\equiv 0(\bmod x^n) \]

整理可以得到：

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

例题：

1.多项式 exp

给定多项式\(A(x),\)求\(e^{A(x)}(\bmod x^n).\)

设\(F(x)\equiv e^{A(x)}(\bmod x^n),\)两边取对数：

\[\ln F(x)\equiv A(x)(\bmod x^n) \]

\[\ln F(x)-A(x)\equiv 0(\bmod x^n) \]

将\(A(x)\)看成常数，设：

\[G(F(x))=\ln F(x)-A(x) \]

则\(G(F(x))\equiv 0\bmod (x^n)\)

\(G'(F(x))=\frac{1}{F(x)}\)

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

\[F(x)\equiv F_0(x)-F_0(x)(\ln F_0(x)-A(x))(\bmod x^n) \]

\[F(x)\equiv F_0(x)(1-\ln F_0(x)+A(x))(\bmod x^n) \]

递归求解\(O(n\log n).\)

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=310000;
const int mod=998244353;
int rev[N],a[N],b[N],c[N],d[N],e[N],f[N],g[N];
int lnb[N],G[N],n,k[N];
inline int add(int x,int y){return (x+y)%mod;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);b>>=1;
	}
	return res;
}
void NTT(int *A,int lim,int tp){
	for(int i=0;i<lim;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int i=1;i<lim;i<<=1){
		int gn=qpow(3,(mod-1)/(i<<1));
		if(tp!=1)gn=qpow(gn,mod-2);
		for(int j=0;j<lim;j+=(i<<1)){
			int G=1;
			for(int k=0;k<i;++k,G=mul(G,gn)){
				int x=A[j+k],y=mul(G,A[i+j+k]);
				A[j+k]=add(x,y);A[i+j+k]=add(x,mod-y);
			}
		}
	}
	if(tp==1)return;
	int inv=qpow(lim,mod-2);
	for(int i=0;i<lim;++i)A[i]=mul(A[i],inv);
}
void Inv(int d,int *a,int *b){
	if(d==1){b[0]=qpow(a[0],mod-2);return;}
	Inv((d+1)>>1,a,b);
	int lim=1,len=0;
	while(lim<(d<<1))lim<<=1,len++;
	for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	for(int i=0;i<d;++i)c[i]=a[i];
	for(int i=d;i<lim;++i)c[i]=0;
	NTT(c,lim,1);NTT(b,lim,1);
	for(int i=0;i<lim;++i)b[i]=1ll*(2-1ll*b[i]*c[i]%mod+mod)%mod*b[i]%mod;
	NTT(b,lim,-1);for(int i=d;i<lim;++i)b[i]=0;
}
inline void Dx(int *a,int *b,int L){for(int i=1;i<L;++i)b[i-1]=mul(i,a[i]);b[L-1]=0;}
inline void Int(int *a,int *b,int L){for(int i=1;i<L;++i)b[i]=mul(a[i-1],qpow(i,mod-2));b[0]=0;}
void Ln(int L,int *a,int *R){
	Dx(a,e,L);Inv(L,a,f);
	int lim=1,len=0;
	while(lim<(L<<1))lim<<=1,len++;
	for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	NTT(e,lim,1);NTT(f,lim,1);
	for(int i=0;i<lim;++i)e[i]=mul(e[i],f[i]);
	NTT(e,lim,-1);Int(e,R,lim);
	for(int i=0;i<lim;++i)e[i]=f[i]=0;
}
void Exp(int d,int *a,int *b){
	if(d==1){b[0]=1;return;}
	Exp((d+1)>>1,a,b);
	Ln(d,b,lnb);
	int lim=1,len=0;
	while(lim<(d<<1))lim<<=1,len++;
	for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	for(int i=0;i<d;++i)lnb[i]=a[i]>=lnb[i]?a[i]-lnb[i]:a[i]-lnb[i]+mod;
	for(int i=d;i<lim;++i)b[i]=lnb[i]=0;
	lnb[0]++;lnb[0]%=mod;
	NTT(lnb,lim,1);NTT(b,lim,1);
	for(int i=0;i<lim;++i)b[i]=mul(b[i],lnb[i]);
	NTT(b,lim,-1);
	for(int i=d;i<lim;++i)b[i]=lnb[i]=0;
}
signed main(){
	scanf("%lld",&n);
	for(int i=0;i<n;++i)scanf("%lld",&a[i]);
	int len=1;
	while(len<=n)len<<=1;
	Exp(len,a,k);
	//Ln(len,a,k);
	for(int i=0;i<n;++i)printf("%lld ",k[i]);
	puts("");
	return 0;
}

2.多项式开根

给定多项式\(A(x),\)求\(B^2(x)\equiv A(x)(\bmod x^n).\)

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

求\(G(B(x))\equiv 0(\bmod x^n).\)

令\(F(x)=B(x).\)套用牛顿迭代公式：

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

注意到，\(A(x)\)是常数，\(G(x)\)的导数\(G'(x)=2B(x)\)

所以原式：

\[F(x)\equiv F_0(x)-\frac{F_0^2(x)-A(x)}{2F_0(x)}(\bmod x^n) \]

\[F(x)\equiv \frac{F_0^2(x)+A(x)}{2F_0(x)} \]

牛顿迭代即可。复杂度\(O(n\log n).\)

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=310000;
const int mod=998244353;
int rev[N],a[N],b[N],c[N],F[N],G[N],ans[N];
int n,inv2,C[N];
inline int add(int x,int y){return (x+y)%mod;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);b>>=1;
	}
	return res;
}
void NTT(int *A,int lim,int tp){
	for(int i=0;i<lim;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int i=1;i<lim;i<<=1){
		int gn=qpow(3,(mod-1)/(i<<1));
		if(tp!=1)gn=qpow(gn,mod-2);
		for(int j=0;j<lim;j+=(i<<1)){
			int G=1;
			for(int k=0;k<i;++k,G=mul(G,gn)){
				int x=A[j+k],y=mul(G,A[i+j+k]);
				A[j+k]=add(x,y);A[i+j+k]=add(x,mod-y);
			}
		}
	}
	if(tp==1)return;
	int inv=qpow(lim,mod-2);
	for(int i=0;i<lim;++i)A[i]=mul(A[i],inv);
}
void Inv(int d,int *a,int *b){
	if(d==1){b[0]=qpow(a[0],mod-2);return;}
	Inv((d+1)>>1,a,b);
	int lim=1,len=0;
	while(lim<(d<<1))lim<<=1,len++;
	for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	for(int i=0;i<d;++i)c[i]=a[i];
	for(int i=d;i<lim;++i)c[i]=0;
	NTT(c,lim,1);NTT(b,lim,1);
	for(int i=0;i<lim;++i)b[i]=1ll*(2-1ll*b[i]*c[i]%mod+mod)%mod*b[i]%mod;
	NTT(b,lim,-1);for(int i=d;i<lim;++i)b[i]=0;
}
void Sqrt(int d,int *a,int *b){
	if(d==1){b[0]=1;return;}
	Sqrt((d+1)>>1,a,b);
	int lim=1,len=0;
	while(lim<(d<<1))lim<<=1,len++;
	for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	Inv(d,b,G);
	for(int i=0;i<d;++i)C[i]=a[i];
	for(int i=d;i<lim;++i)C[i]=0;
	NTT(C,lim,1);NTT(G,lim,1);
	for(int i=0;i<lim;++i)C[i]=mul(G[i],C[i]);
	NTT(C,lim,-1);//G=A/F_0
	for(int i=0;i<d;++i)b[i]=1ll*(b[i]+C[i])%mod*inv2%mod;
	for(int i=d;i<lim;++i)b[i]=G[i]=C[i]=0;
	for(int i=0;i<d;++i)C[i]=G[i]=0;
}
signed main(){
	scanf("%lld",&n);
	for(int i=0;i<n;++i)scanf("%lld",&a[i]);
	inv2=qpow(2,mod-2);
	int len=1;
	while(len<=n)len<<=1;
	Sqrt(len,a,ans);
	for(int i=0;i<n;++i)printf("%lld ",ans[i]);
	puts("");
	return 0;
}