zoukankan      html  css  js  c++  java
  • 多项式开根

    给定(n-1)次多项式(A(x)),求一个在(mod\, x^n)意义下的多项式(B(x)),使得(B(x)^2≡A(x) (mod\, x^n)),取零次项系数最小的作为答案

    系数(mod\, 998244353)

    根据题意:

    (B(x)^2-A(x)≡0 (mod\, x^n))

    (G(B(x))=B(x)^2-A(x)),套牛顿迭代

    (B(x)≡B_0(x)-frac{G(B_0(x))}{G'(B_0(x))})

    整理一下

    (B(x)=frac{A(x)+B_0(x)^2}{2B_0(x)})

    直接乘常数有点大,稍微变化一下

    (B(x)=frac{1}{2}(A(x)B_0(x)^{-1}+B_0(x)))

    递归求解

    #include<bits/stdc++.h>
    using namespace std;
    namespace red{
    #define int long long
    #define y1 miao
    #define eps (1e-8)
    	inline int read()
    	{
    		int x=0;char ch,f=1;
    		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
    		if(ch=='-') f=0,ch=getchar();
    		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    		return f?x:-x;
    	}
    	const int N=4e5+10,p=998244353;
    	int n,inv2;
    	int a[N],b[N],c[N],pos[N];
    	int invb[N];
    	inline int fast(int x,int k)
    	{
    		int ret=1;
    		while(k)
    		{
    			if(k&1) ret=ret*x%p;
    			x=x*x%p;
    			k>>=1;
    		}
    		return ret;
    	}
    	inline void ntt(int limit,int *a,int inv)
    	{
    		for(int i=0;i<limit;++i)
    			if(i<pos[i]) swap(a[i],a[pos[i]]);
    		for(int mid=1;mid<limit;mid<<=1)
    		{
    			int Wn=fast(3,(p-1)/(mid<<1));
    			for(int r=mid<<1,j=0;j<limit;j+=r)
    			{
    				int w=1;
    				for(int k=0;k<mid;++k,w=w*Wn%p)
    				{
    					int x=a[j+k],y=w*a[j+k+mid]%p;
    					a[j+k]=(x+y)%p;
    					a[j+k+mid]=(x-y+p)%p;
    				}
    			}
    		}
    		if(inv) return;
    		inv=fast(limit,p-2);reverse(a+1,a+limit);
    		for(int i=0;i<limit;++i) a[i]=a[i]*inv%p;
    	}
    	inline void poly_inv(int pw,int *a,int *b)
    	{
    		if(pw==1){b[0]=fast(a[0],p-2);return;}
    		poly_inv((pw+1)>>1,a,b);
    		int len=0,limit=1;
    		while(limit<(pw<<1)) limit<<=1,++len;
    		for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
    		for(int i=0;i<pw;++i) c[i]=a[i];
    		for(int i=pw;i<limit;++i) c[i]=0;
    		ntt(limit,c,1);ntt(limit,b,1);
    		for(int i=0;i<limit;++i) b[i]=((2-c[i]*b[i]%p)+p)%p*b[i]%p;
    		ntt(limit,b,0);
    		for(int i=pw;i<limit;++i) b[i]=0;
    	}
    	inline void poly_sqrt(int pw,int *a,int *b)
    	{
    		if(pw==1){b[0]=1;return;}
    		poly_sqrt((pw+1)>>1,a,b);
    		for(int i=0;i<(pw<<1);++i) invb[i]=0;//记得清空数组
    		poly_inv(pw,b,invb);
    		int len=0,limit=1;
    		while(limit<(pw<<1)) limit<<=1,++len;
    		for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
    		for(int i=0;i<pw;++i) c[i]=a[i];
    		for(int i=pw;i<limit;++i) c[i]=0;
    		ntt(limit,c,1);ntt(limit,invb,1);
    		for(int i=0;i<limit;++i) c[i]=c[i]*invb[i]%p;
    		ntt(limit,c,0);
    		for(int i=0;i<pw;++i) b[i]=(b[i]+c[i])%p*inv2%p;
    		for(int i=pw;i<limit;++i) b[i]=0;
    	}
    	inline void main()
    	{
    		n=read();inv2=fast(2,p-2);
    		for(int i=0;i<n;++i) a[i]=read();
    		poly_sqrt(n,a,b);
    		for(int i=0;i<n;++i) printf("%lld ",b[i]);
    	}
    }
    signed main()
    {
    	red::main();
    return 0;
    }
    
  • 相关阅读:
    (转)SVN命令详解
    CentOS 64位系统安装32位兼容库
    Eclipse4.4 Tomcat插件下载地址
    记:Android 安装apk的代码实现
    记: Android adb远程调试
    记:Ubuntu14.04 Android加载项目失败
    记:Ubuntu 14.04 安装32位库支持库失败
    记:Android 服务站 问题记录与解决方案
    记:Android 知识点整理 20140329
    微信小程序和微信公众号的id是一个吗
  • 原文地址:https://www.cnblogs.com/knife-rose/p/12131591.html
Copyright © 2011-2022 走看看