zoukankan      html  css  js  c++  java
  • 【学习笔记】牛顿迭代

    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(mod x^n) ]

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

    考虑已经解决了:

    [G(F_0(x))equiv 0(mod x^{leftlceil frac{n}{2} ight ceil} ) ]

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

    在这里泰勒展开一下:

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

    注意到当(ige 2)((F(x)-F_0(x))^i)的最低非(0)次项的次数是严格大于(2leftlceilfrac{n}{2} ight ceil,)所以:

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

    注意到由题设得:

    [G(F(x))equiv 0(mod x^n) ]

    所以:

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

    整理可以得到:

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

    例题:

    1.多项式 exp

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

    (F(x)equiv e^{A(x)}(mod x^n),)两边取对数:

    [ln F(x)equiv A(x)(mod x^n) ]

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

    (A(x))看成常数,设:

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

    (G(F(x))equiv 0mod (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)}}(mod x^n) ]

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

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

    递归求解(O(nlog 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)(mod x^n).)

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

    (G(B(x))equiv 0(mod x^n).)

    (F(x)=B(x).)套用牛顿迭代公式:

    [F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))}(mod 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)}(mod x^n) ]

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

    牛顿迭代即可。复杂度(O(nlog 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;
    }
    
  • 相关阅读:
    redis相关
    Ubuntu安装之python开发
    Shell编程实战
    saltstack高效运维
    docker网络
    docker入门
    python学习博客地址集合。。。
    vue+uwsgi+nginx部署路飞学城
    部署你的CRM程序
    Java Netty教程(目录)
  • 原文地址:https://www.cnblogs.com/h-lka/p/14489483.html
Copyright © 2011-2022 走看看