多项式牛顿迭代 学习笔记

零. 前置芝士

在正式进入牛顿迭代的学习前，我们需要预先了解 (Taylor) Swift 展开，以便处理我们推导过程中遇到的问题。

所谓泰勒展开，是一个近似模拟原函数的一个公式，和拉格朗日插值大概做的是同一件事情(？)，不过拉格朗日插值是通过多项式点值表达求出的，而泰勒展开是基于原函数表达式已知，利用另一种方式来近似它以便计算&研究的。

想要复制这段曲线，首先得找一个切入点，可以是这条曲线最左端的点，也可以是最右端的点，anyway，可以是这条线上任何一点。他选了最左边的点。由于这段曲线过 ((0,1)) 这个点，仿造的第一步，就是让仿造的曲线也过这个点，完成了仿造的第一步，很粗糙，甚至完全看不出来这俩有什么相似的地方，那就继续细节化。

开始考虑曲线的变化趋势，即导数，保证在此处的导数相等。经历了第二步，现在起始点相同了，整体变化趋势相近了，可能看起来有那么点意思了。想进一步精确化，应该考虑凹凸性。

高中学过：表征图像的凹凸性的参数为“导数的导数”。所以，下一步就让二者的导数的导数相等。起始点相同，增减性相同，凹凸性相同后，仿造的函数更像了。如果再继续细化下去，应该会无限接近。所以泰勒认为“仿造一段曲线，要先保证起点相同，再保证在此处导数相同，继续保证在此处的导数的导数相同……”

(以上内容摘自知乎)

基于以上文字的描述，我们来考虑严谨的计算证明：假设我们求(n)次导来近似这个函数，那么我们需要构造一个多项式函数 (g(x))，使其与原函数 (f(x)) 在某一点的初始值相等，且求 (i（i in [1,n]）) 阶导相等
如果这个函数 (g(x)) 可求 (n) 阶导，那么它必然是一个 (n)次多项式，不妨设为 (sumlimits_{i=0}^n{a_ix^i})（貌似很像普通生成函数）
现在假设我们选定的点是 ((0,f(0)=a_0)) ，那么有：(f^{(n)}(0)=g^{(n)}(0))，由于求 (n) 阶导时只有最后一项非零，为(n!a_n)，那么应有：(a_n=frac{f^{(n)}(0)}{n!})

同理有：$$a_i=frac{f^{(i)}(0)}{i!}(i in [1,n])$$

这就是泰勒展开的精髓，所以原函数的近似函数(g(x))为：(g(x)=g(0)+frac{f^{(1)}(0)}{1!}x^1+frac{f^{(2)}(0)}{2!}x^2+...+frac{f^{(n)}(0)}{n!}x^n)
如果 (n) 趋于无穷大，那么这两条曲线无限近似，同时我们将 (0) 换为 (x_0) ，有

[f(x)=g(x)=g(x_0)+frac{f^{'}(x_0)}{1!}(x-x_0)^1+frac{f^{''}(x_0)}{2!}(x-x_0)^2+...+frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+... ]

这就是泰勒展开公式

一.(Newton's) (Method)

理解完上面的泰勒展开，我们进入正题：

牛顿迭代法解决如下问题：已知多项式(g(x))表达式，且有(g(f(x))=0(mod x^n))（不会打同余符号，以下均用等号表示），求(f(x))
考虑倍增（以下叙述用类似数学归纳法的叙述方式展开）
在(n=1)时，单独计算([x^0]g(f(x))=0)的解单独求出
在倍增的前提下，若我们现在已经得到模(x^{lceil frac{n}{2} ceil})意义下的解 (f_0(x))
考虑将 (g(f(x)))在(f_0(x))处泰勒展开，有：(sumlimits_{i=1}^{+ infty}{frac{g^{(i)}f_0(x)}{i!}(f(x)-f_0(x))^i}=0 (mod x^n))
容易得到 (f(x)-f_0(x))的最低非零项次数为 (lceil frac{n}{2} ceil)，有：对于任意(2 leq i)，((f(x)-f_0(x))^i=0(mod x^n))，那么重新写柿子，移项，有：

[f(x)=f_0(x)-frac{g(f_0(x))}{g'(f_0(x))}(mod x^n) ]

这样就得到了 (f(x)) 的表达式

二.牛顿迭代应用1：多项式求逆

• 给定多项式为 (h(x))，且(f(x) imes h(x) = 1 (mod x^n))，欲求 (f(x))
  有：(h(x) = h(x) (mod m))，即(f^{-1}=h(x) (mod m))，移项有(f^{-1}(x)-h(x) (mod m))
  由牛顿迭代，有(f(x)=f_0(x)- frac{f_0^{-1}(x)-h(x)}{- frac{1}{f_0^2(x)}} (mod m))，整理有：

[f(x)=2f_0(x)-f_0^2(x)h(x) (mod m) ]

由此得到(f(x))递推式，利用(ntt)可将复杂度优化至(O(n log n))
代码如下：

#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
	typedef long long ll;
	typedef double db;
	#define pf printf
	#define pc putchar
	#define fr(i,x,y) for(register ll i=(x);i<=(y);++i)
	#define pfr(i,x,y) for(register ll i=(x);i>=(y);--i)
	#define go(u) for(ll i=head[u];i;i=e[i].nxt)
	#define enter pc('
')
	#define space pc(' ')
	#define fir first
	#define sec second
	#define MP make_pair
	const ll inf=0x3f3f3f3f;
	const ll inff=1e15;
	inline ll read()
	{
		ll sum=0,f=1;
		char ch=0;
		while(!isdigit(ch))
		{
			if(ch=='-') f=-1;
			ch=getchar();
		}
		while(isdigit(ch))
		{
			sum=sum*10+(ch^48);
			ch=getchar();
		}
		return sum*f;
	}
	inline void write(ll x)
	{
		if(x<0)
		{
			x=-x;
			pc('-');
		}
		if(x>9) write(x/10);
		pc(x%10+'0');
	}
	inline void writeln(ll x)
	{
		write(x);
		enter;
	}
	inline void writesp(ll x)
	{
		write(x);
		space;
	}
}
using namespace my_std;
const ll N=1e7+50,G=3,mod=998244353;
ll n,limit=1,a[N],b[N],tmp[N],l,r[N];
inline ll ksmod(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=(ans*a)%mod;
		a=(a*a)%mod;
		b>>=1;
	}
	return ans%mod;
}
inline void ntt(ll *a,ll limit,ll pd)
{
	fr(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
	for(ll i=1;i<limit;i<<=1)
	{
		ll gn=ksmod(G,(mod-1)/(i<<1));
		for(ll j=0;j<limit;j+=(i<<1))
		{
			ll g=1;
			for(ll k=0;k<i;++k,g=(g*gn)%mod)
			{
			 	ll x=a[j+k],y=g*a[j+k+i]%mod;
				a[j+k]=(x+y)%mod,a[j+k+i]=(x-y+mod)%mod;
			}
		}
	}
	if(pd==1) return ;
	ll inv=ksmod(limit,mod-2);
	reverse(a+1,a+limit);
	fr(i,0,limit-1) a[i]=a[i]*inv%mod;
}
void solve(ll n,ll *a,ll *b)
{
	if(n==1)
	{
		b[0]=ksmod(a[0],mod-2);
		return ;
	}
	solve((n+1)>>1,a,b);
	static ll l=0,limit=1;
	while(limit<(n<<1)) limit<<=1,++l;
	fr(i,1,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	fr(i,0,n-1) tmp[i]=a[i];
	fr(i,n,limit-1) tmp[i]=0;
	ntt(b,limit,1),ntt(tmp,limit,1);
	fr(i,0,limit-1) b[i]=(2*b[i]%mod-(b[i]*b[i]%mod*tmp[i])%mod+mod)%mod;
	ntt(b,limit,-1);
	fr(i,n,limit-1) b[i]=0; 
}
int main(void)
{
	n=read();
	fr(i,0,n-1) a[i]=read();
	solve(n,a,b);
	fr(i,0,n-1) writesp(b[i]);
	return 0;
}

三.牛顿迭代应用2：多项式 (ln)

• 给定多项式 (f(x))，求 (g(x)=ln f(x)(mod x^n))
  其实这玩意用不着牛顿迭代，但是过程中需要用到多项式求逆所以只是来凑数的
  对这个式子两边分别求导，有： (g'(x)=frac{f'(x)}{f(x)} (mod x^n))，
  同时积分，有：

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

于是对 (f(x)) 求逆、求导、(ntt)相乘，积分就可以得到 (g(x))

• 多项式求导：

inline void polyderi(ll n,ll *a)
{
	fr(i,1,n) a[i-1]=(a[i]*i)%mod;
	a[n]=0;
}

• 多项式积分

inline void polyinte(ll n,ll *a)
{
	pfr(i,n,1) a[i]=(a[i-1]*ksmod(i,mod-2))%mod;
	a[0]=0;
}

代码：

#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
	typedef long long ll;
	typedef double db;
	#define pf printf
	#define pc putchar
	#define fr(i,x,y) for(register ll i=(x);i<=(y);++i)
	#define pfr(i,x,y) for(register ll i=(x);i>=(y);--i)
	#define go(u) for(ll i=head[u];i;i=e[i].nxt)
	#define enter pc('
')
	#define space pc(' ')
	#define fir first
	#define sec second
	#define MP make_pair
	const ll inf=0x3f3f3f3f;
	const ll inff=1e15;
	inline ll read()
	{
		ll sum=0,f=1;
		char ch=0;
		while(!isdigit(ch))
		{
			if(ch=='-') f=-1;
			ch=getchar();
		}
		while(isdigit(ch))
		{
			sum=sum*10+(ch^48);
			ch=getchar();
		}
		return sum*f;
	}
	inline void write(ll x)
	{
		if(x<0)
		{
			x=-x;
			pc('-');
		}
		if(x>9) write(x/10);
		pc(x%10+'0');
	}
	inline void writeln(ll x)
	{
		write(x);
		enter;
	}
	inline void writesp(ll x)
	{
		write(x);
		space;
	}
}
using namespace my_std;
namespace polynomial
{
	const ll N=1e7+50,G=3,mod=998244353;
	ll r[N],tmp[N];
	inline ll ksmod(ll a,ll b)
	{
		ll ans=1;
		while(b)
		{
			if(b&1) ans=(ans*a)%mod;
			a=(a*a)%mod;
			b>>=1;
		}
		return ans%mod;
	}
	inline void ntt(ll *a,ll limit,ll pd)
	{
		fr(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
		for(ll i=1;i<limit;i<<=1)
		{
			ll gn=ksmod(G,(mod-1)/(i<<1));
			for(ll j=0;j<limit;j+=(i<<1))
			{
				ll g=1;
				for(ll k=0;k<i;++k,g=(g*gn)%mod)
				{
				 	ll x=a[j+k],y=g*a[j+k+i]%mod;
					a[j+k]=(x+y)%mod,a[j+k+i]=(x-y+mod)%mod;
				}
			}
		}
		if(pd==1) return ;
		ll inv=ksmod(limit,mod-2);
		reverse(a+1,a+limit);
		fr(i,0,limit-1) a[i]=a[i]*inv%mod;
	}
	inline void polyinv(ll n,ll *a,ll *b)
	{
		if(n==1)
		{
			b[0]=ksmod(a[0],mod-2);
			return ;
		}
		polyinv((n+1)>>1,a,b);
		static ll l=0,limit=1;
		while(limit<(n<<1)) limit<<=1,++l;
		fr(i,1,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
		fr(i,0,n-1) tmp[i]=a[i];
		fr(i,n,limit-1) tmp[i]=0;
		ntt(b,limit,1),ntt(tmp,limit,1);
		fr(i,0,limit-1) b[i]=(2*b[i]%mod-(b[i]*b[i]%mod*tmp[i])%mod+mod)%mod;
		ntt(b,limit,-1);
		fr(i,n,limit-1) b[i]=0; 
	}
	inline void polyderi(ll n,ll *a)
	{
		fr(i,1,n) a[i-1]=(a[i]*i)%mod;
		a[n]=0;
	}
	inline void polyinte(ll n,ll *a)
	{
		pfr(i,n,1) a[i]=(a[i-1]*ksmod(i,mod-2))%mod;
		a[0]=0;
	}
	ll b[N],ans[N];
	inline void polyln(ll n,ll *a)
	{
		memcpy(b,a,sizeof(b));
		polyderi(n,a),polyinv(n,b,ans);
		ll limit=1,l=0;
		while(limit<=(n<<1)) limit<<=1,++l;
		fr(i,0,limit-1) b[i]=ans[i];
		fr(i,1,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
		ntt(a,limit,1),ntt(b,limit,1);
		fr(i,0,limit) a[i]=(a[i]*b[i])%mod;
		ntt(a,limit,-1);
		polyinte(n,a);
		fr(i,n+1,limit) a[i]=0;
	}
}
using namespace polynomial;
ll n,l,limit=1,f[N],g[N];
int main(void)
{
	n=read();
	fr(i,0,n-1) f[i]=read();
	while(limit<=(n<<1)) limit<<=1,++l;
	polyln(n,f);
	fr(i,0,n-1) writesp(f[i]);
	return 0;
}

(To) (Be) (Continued...)