多项式牛顿迭代 学习笔记

零. 前置芝士

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

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

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

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

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

(以上内容摘自知乎)

基于以上文字的描述，我们来考虑严谨的计算证明：假设我们求\(n\)次导来近似这个函数，那么我们需要构造一个多项式函数 \(g(x)\)，使其与原函数 \(f(x)\) 在某一点的初始值相等，且求 \(i（i \in [1,n]）\) 阶导相等
如果这个函数 \(g(x)\) 可求 \(n\) 阶导，那么它必然是一个 \(n\)次多项式，不妨设为 \(\sum\limits_{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} \rceil}\)意义下的解 \(f_0(x)\)
考虑将 \(g(f(x))\)在\(f_0(x)\)处泰勒展开，有：\(\sum\limits_{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} \rceil\)，有：对于任意\(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) \times 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('\n')
	#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('\n')
	#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...\)

posted @ 2020-07-23 20:45  L_G_J  阅读(326)  评论(1)  编辑  收藏  举报