zoukankan html css js c++ java

Note_4.1

2019/4/1 奇奇怪怪的笔记

多项式除法

问题描述

给定(n)次多项式(A(x))和(m)次多项式(B(x))

求：

[A(x)=B(x)*C(x)+R(x) ]

我们要求(C(x))的次数必须是(n-m)，(R(x))的次数小于(m) ，所以不能简单地用求逆解决

方法

首先我们考虑：

[x^nF(frac{1}{x})=sum_{i=0}^{n}a_ix^{n-i+1} ]

根据(A=B*C+R)，我们可以得到：

[egin{equation} A(frac{1}{x})=B(frac{1}{x})*C(frac{1}{x})+R(frac{1}{x})\ x^nA(frac{1}{x})=x^mB(frac{1}{x})*x^{n-m}C(frac{1}{x})+x^nR(frac{1}{x})\ x^nA(frac{1}{x})=x^mB(frac{1}{x})*x^{n-m}C(frac{1}{x}) mod x^{n-m} end{equation} ]

只要求逆算出(x^{n-m}C(frac{1}{x}))，翻转就可以得到(C(x))

代码

#include <bits/stdc++.h>
using namespace std;
#define reg register
#define ll long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}
const int N=1<<18|5,P=998244353,g=3,invg=332748118;
int pos[N];
#define Mul(x,y) (1ll*(x)*(y)%P)
inline int fpow(int x,int m){int r=1;for(;m;m>>=1,x=Mul(x,x))if(m&1)r=Mul(r,x);return r;}
void init(int len)
{
	reg int i;
	for(i=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(len>>1));
}
void NTT(int *a,int n,int type)
{
	register int i,j,k,w,wn,X,Y;
	for(i=0;i<n;++i) if(i<pos[i]) swap(a[i],a[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		wn=fpow(type>0?g:invg,(P-1)/(i<<1));
		for(j=0;j<n;j+=i<<1)for(w=1,k=0;k<i;++k,w=Mul(w,wn))
		{
			X=a[j+k],Y=Mul(w,a[j+k+i]);
			a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
		}
	}
	reg int invn=fpow(n,P-2);
	if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn);
}
inline void Dec(int *a,int *b,int *c,int n){for(reg int i=0;i<n;++i)c[i]=(a[i]-b[i]%P+P)%P;}
inline void Pro(int *a,int *b,int *c,int n)
{
	static int len,A[N],B[N];
	for(len=1;len<(n<<1);len<<=1);
	memcpy(A,a,sizeof(int[n]));memcpy(B,b,sizeof(int[n]));
	memset(A+n,0,sizeof(int[len-n]));memset(B+n,0,sizeof(int[len-n]));
	reg int i;init(len);
	NTT(A,len,1);NTT(B,len,1);
	for(i=0;i<len;++i) c[i]=Mul(A[i],B[i]);
	NTT(c,len,-1);
	memset(c+n,0,sizeof(int[len-n]));
}
void _Inv(int *A,int *b,int n)
{
	if(n==1){b[0]=fpow(A[0],P-2);return;}
	int t=(n+1)>>1;_Inv(A,b,t);
    static int len,a[N];
    for(len=1;len<(n<<1);len<<=1);
    memcpy(a,A,sizeof(int[n]));
    memset(a+n,0,sizeof(int[len-n]));memset(b+t,0,sizeof(int[len-t]));
    reg int i;init(len);
    NTT(a,len,1),NTT(b,len,1);
    for(i=0;i<len;++i) b[i]=(Mul(2ll,b[i])-Mul(a[i],Mul(b[i],b[i]))+P)%P;
    NTT(b,len,-1);memset(b+n,0,sizeof(int[len-n]));
}
inline void Inv(const int *A, int *b, const int n)
{
    static int a[N];
    memcpy(a,A,sizeof(int[n]));
    _Inv(a,b,n);
}
inline void Div(int *A,int *B,int *c,int *r,int n,int m)
{
	static int a[N],b[N],inv_b[N];
    memcpy(a,A,sizeof(int[n]));memcpy(b,B,sizeof(int[m]));
    reverse(a,a+n),reverse(b,b+m);
    int t=n-m+1;
    if(t>=m) memset(b+m,0,sizeof(int[t-m]));
    Inv(b,inv_b,t);Pro(a,inv_b,c,t);
    reverse(c,c+t);reverse(a,a+n),reverse(b,b+m);
    Pro(b,c,b,n);Dec(a,b,r,n);
}
int a[N],b[N],c[N],r[N];
int main()
{
	reg int i,n=read()+1,m=read()+1;
	for(i=0;i<n;++i) a[i]=read();
	for(i=0;i<m;++i) b[i]=read();
	Div(a,b,c,r,n,m);
	for(i=0;i<=n-m;++i) printf("%d ",c[i]);puts("");
	for(i=0;i<m-1;++i) printf("%d ",r[i]);
	return 0;
}

多项式对数函数

问题描述

给出一个(n-1)次的多项式(A(x))，求(B(x)≡ln A(x) mod 998244353)

方法

考虑对方程两边求导：

[B'(x)≡frac{A'(x)}{A(x)} mod 998244353 ]

所以我们只需要：

对(A(x))求逆，得到(A^{-1}(x))
对(A(x))求导，得到(A'(x))
多项式乘法，(B'(x)=A'(x)*A^{-1}(x))
对(B'(x))积分，得到(B(x))

怎么求导和积分？ 参照人教版高中选修(2-1)

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define reg register
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
const int N=1<<18|5,P=998244353,g=3,invg=332748118;
#define Mul(x,y) (1ll*(x)*(y)%P)
#define ri reg int i
int _[N],pos[N];
int fpow(int x,int m){int r=1;for(;m;m>>=1,x=Mul(x,x))if(m&1)r=Mul(r,x);return r;}
void der(int *a,int *b,int n){for(ri=0;i<n-1;++i) b[i]=Mul(i+1,a[i+1]);b[n-1]=0;}
void inte(int *a,int *b,int n)
{
	ri;for(_[1]=i=1;i<n;++i,_[i]=Mul(_[P%i],P-P/i));
	for(i=n-1;i;--i)b[i]=Mul(a[i-1],_[i]);b[0]=0;
}
void init(int len)
{
	reg int i;
	for(i=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(len>>1));
}
void NTT(int *a,int n,int type)
{
	register int i,j,k,w,wn,X,Y;
	for(i=0;i<n;++i) if(i<pos[i]) swap(a[i],a[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		wn=fpow(type>0?g:invg,(P-1)/(i<<1));
		for(j=0;j<n;j+=i<<1)for(w=1,k=0;k<i;++k,w=Mul(w,wn))
		{
			X=a[j+k],Y=Mul(w,a[j+k+i]);
			a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
		}
	}
	reg int invn=fpow(n,P-2);
	if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn);
}
void Comb(int *a,int *b,int *c,int n)
{
	static int A[N],B[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(A,a,sizeof(int[n]));memcpy(B,b,sizeof(int[n]));
	memset(A+n,0,sizeof(int[len-n]));memset(B+n,0,sizeof(int[len-n]));
	reg int i;init(len);NTT(A,len,1);NTT(B,len,1);
	for(i=0;i<len;++i) c[i]=Mul(A[i],B[i]);NTT(c,len,-1);
	memset(c+n,0,sizeof(int[len-n]));
}
void _I(int *A,int *b,int n)
{
	if(n==1) return(void)(b[0]=fpow(A[0],P-2));
	int t=(n+1)>>1;_I(A,b,t);
	static int a[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(a,A,sizeof(int[n]));memset(a+n,0,sizeof(int[n]));
	memset(b+t,0,sizeof(int[len-t]));init(len);
	reg int i;NTT(a,len,1);NTT(b,len,1);
	for(i=0;i<len;++i)b[i]=(Mul(2ll,b[i])-Mul(a[i],Mul(b[i],b[i]))+P)%P;
	NTT(b,len,-1);memset(b+n,0,sizeof(int[len-n]));
}
void Inv(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_I(a,b,n);
}
void ln(int *A,int *B,int n)
{
	static int a[N],da[N];
	memcpy(a,A,sizeof(int[n]));memcpy(da,A,sizeof(int[n]));
	Inv(a,a,n);der(da,da,n);Comb(da,a,B,n);inte(B,B,n);
}
int A[N],B[N],n;
int main()
{
	n=read();
	reg int i;
	for(i=0;i<n;++i) A[i]=read();
	ln(A,B,n);
	for(i=0;i<n;++i) printf("%d ",B[i]);
	return 0;
}

牛顿迭代法

迭代关系式

目标：求方程(f(x)=0)的根

假设当前，我们得到的近似值是(x_n)，要求得它的(n+1)次近似值：

过点((x_n,f(x_n)))做曲线的切线(L:y=f(x_n)+f'(x_n)(x-x0))

(L)于(x)轴的交点是((x_n-frac{f(x_n)}{f'(x_n)},0))

我们定义牛顿迭代公式就是

[x_{n+1}=x_n-frac{f(x_n)}{f'(x_n)} ]

神仙证明：“如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。并且，如果不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍”

多项式

牛顿迭代在多项式中同样可行。

假设我们要求(G(F(x))≡0mod x^n)

考虑已经求出(G(F_0(x)) ≡ 0 mod x^{⌈frac{n}{2}⌉})把(G(F(x)))在(F_0(z))处进行泰勒展开

[egin{equation} G(F(x))=G(F_0(x))+frac{G'(F_0(x))}{1!}(F(x)-F_0(x))+frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+...+frac{G^{(n)}(F_0(x))}{n!}(F(x)-F_0(x))^n+... end{equation} ]

可知，(F(x))和(F_0(x))模(x^{⌈frac{n}{2}⌉})同余，所以((F(x)-F_0(x))^2≡0mod x^n)

又因为

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

所以

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

多项式指数函数

问题描述

给出一个(n-1)次的多项式(A(x))，求(B(x)≡e^{A(x)} mod x^{n})，(a_0=0)，然后(mod 998244353)

方法

等同于：

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

对(F)求导

[F'(B(x))=frac{1}{B(x)} ]

套用牛顿迭代

[egin{equation} egin{split} B(x)&leftarrow B_0(x)-frac{F(B_0(x))}{F'(B_0(x))}\ &=B_0(x)-(ln(B_0(x))-A(x))B_0(x)\ &=B_0(x)(-ln(B_0(x))+A(x)+1) end{split} end{equation} ]

考虑边界情况，因为本题限定了(a_0=0)，而(e^0=1)，所以此时(b_0=1)

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define reg register
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return x*f;
}
const int N=1<<18|5,P=998244353,g=3,invg=332748118;
#define Mul(x,y) (1ll*(x)*(y)%P)
#define ri reg int i
int _[N],pos[N];
int fpow(int x,int m){int r=1;for(;m;m>>=1,x=Mul(x,x))if(m&1)r=Mul(r,x);return r;}
void der(int *a,int *b,int n){for(ri=0;i<n-1;++i) b[i]=Mul(i+1,a[i+1]);b[n-1]=0;}
void inte(int *a,int *b,int n)
{
	ri;for(_[1]=i=1;i<n;++i,_[i]=Mul(_[P%i],P-P/i));
	for(i=n-1;i;--i)b[i]=Mul(a[i-1],_[i]);b[0]=0;
}
void init(int len)
{
	reg int i;
	for(i=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(len>>1));
}
void NTT(int *a,int n,int type)
{
	register int i,j,k,w,wn,X,Y;
	for(i=0;i<n;++i) if(i<pos[i]) swap(a[i],a[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		wn=fpow(type>0?g:invg,(P-1)/(i<<1));
		for(j=0;j<n;j+=i<<1)for(w=1,k=0;k<i;++k,w=Mul(w,wn))
		{
			X=a[j+k],Y=Mul(w,a[j+k+i]);
			a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
		}
	}
	reg int invn=fpow(n,P-2);
	if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn);
}
void Comb(int *a,int *b,int *c,int n)
{
	static int A[N],B[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(A,a,sizeof(int[n]));memcpy(B,b,sizeof(int[n]));
	memset(A+n,0,sizeof(int[len-n]));memset(B+n,0,sizeof(int[len-n]));
	reg int i;init(len);NTT(A,len,1);NTT(B,len,1);
	for(i=0;i<len;++i) c[i]=Mul(A[i],B[i]);NTT(c,len,-1);
	memset(c+n,0,sizeof(int[len-n]));
}

void _I(int *A,int *b,int n)
{
	if(n==1) return(void)(b[0]=fpow(A[0],P-2));
	int t=(n+1)>>1;_I(A,b,t);
	static int a[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(a,A,sizeof(int[n]));memset(a+n,0,sizeof(int[len-n]));
	memset(b+t,0,sizeof(int[len-t]));init(len);
	reg int i;NTT(a,len,1);NTT(b,len,1);
	for(i=0;i<len;++i)b[i]=(Mul(2ll,b[i])-Mul(a[i],Mul(b[i],b[i]))+P)%P;
	NTT(b,len,-1);memset(b+n,0,sizeof(int[len-n]));
}
void Inv(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_I(a,b,n);
}
void ln(int *A,int *B,int n)
{
	static int a[N],da[N];
	Inv(A,a,n);der(A,da,n);Comb(da,a,B,n);inte(B,B,n);
}
void _e(int *a,int *b,int n)
{
	if(n==1)return(void)(b[0]=1);
	int t=(n+1)>>1;_e(a,b,t);
	static int xxx[N];
	ln(b,xxx,n);reg int i;
	for(i=0;i<n;++i)xxx[i]=(-xxx[i]+a[i]+P)%P;
	xxx[0]=(xxx[0]+1)%P;
	Comb(b,xxx,b,n);
	
}
void exp(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_e(a,b,n);
}
int A[N],B[N],n;
int main()
{
	n=read();
	reg int i;
	for(i=0;i<n;++i) A[i]=read();
	exp(A,B,n);
	for(i=0;i<n;++i) printf("%d ",B[i]);
	return 0;
}

多项式快速幂

问题描述

给出一个(n-1)次的多项式(A(x))，求(B(x)≡A^k(x) mod x^{n})，(a_0=1)，然后(mod 998244353)

方法

等同于：

[B(x)=e^{kln A(x)} ]

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define reg register
const int N=1<<18|5,P=998244353,g=3,invg=332748118;
inline int read()
{
	ll x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x*10ll+ch-'0')%P;ch=getchar();}
	return x*f;
}
#define Mul(x,y) (1ll*(x)*(y)%P)
#define ri reg int i
int _[N],pos[N];
int fpow(int x,int m){int r=1;for(;m;m>>=1,x=Mul(x,x))if(m&1)r=Mul(r,x);return r;}
void der(int *a,int *b,int n){for(ri=0;i<n-1;++i) b[i]=Mul(i+1,a[i+1]);b[n-1]=0;}
void inte(int *a,int *b,int n)
{
	ri;for(_[1]=i=1;i<n;++i,_[i]=Mul(_[P%i],P-P/i));
	for(i=n-1;i;--i)b[i]=Mul(a[i-1],_[i]);b[0]=0;
}
void init(int len)
{
	reg int i;
	for(i=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(len>>1));
}
void NTT(int *a,int n,int type)
{
	register int i,j,k,w,wn,X,Y;
	for(i=0;i<n;++i) if(i<pos[i]) swap(a[i],a[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		wn=fpow(type>0?g:invg,(P-1)/(i<<1));
		for(j=0;j<n;j+=i<<1)for(w=1,k=0;k<i;++k,w=Mul(w,wn))
		{
			X=a[j+k],Y=Mul(w,a[j+k+i]);
			a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
		}
	}
	reg int invn=fpow(n,P-2);
	if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn);
}
void Comb(int *a,int *b,int *c,int n)
{
	static int A[N],B[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(A,a,sizeof(int[n]));memcpy(B,b,sizeof(int[n]));
	memset(A+n,0,sizeof(int[len-n]));memset(B+n,0,sizeof(int[len-n]));
	reg int i;init(len);NTT(A,len,1);NTT(B,len,1);
	for(i=0;i<len;++i) c[i]=Mul(A[i],B[i]);NTT(c,len,-1);
	memset(c+n,0,sizeof(int[len-n]));
}
void _I(int *A,int *b,int n)
{
	if(n==1) return(void)(b[0]=fpow(A[0],P-2));
	int t=(n+1)>>1;_I(A,b,t);
	static int a[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(a,A,sizeof(int[n]));memset(a+n,0,sizeof(int[len-n]));
	memset(b+t,0,sizeof(int[len-t]));init(len);
	reg int i;NTT(a,len,1);NTT(b,len,1);
	for(i=0;i<len;++i)b[i]=(Mul(2ll,b[i])-Mul(a[i],Mul(b[i],b[i]))+P)%P;
	NTT(b,len,-1);memset(b+n,0,sizeof(int[len-n]));
}
void Inv(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_I(a,b,n);
}
void ln(int *A,int *B,int n)
{
	static int a[N],da[N];
	memcpy(a,A,sizeof(int[n]));memcpy(da,A,sizeof(int[n]));
	Inv(a,a,n);der(da,da,n);Comb(da,a,B,n);inte(B,B,n);
}
void _e(int *a,int *b,int n)
{
	if(n==1)return(void)(b[0]=1);
	int t=(n+1)>>1;_e(a,b,t);
	static int xxx[N];
	ln(b,xxx,n);reg int i;
	for(i=0;i<n;++i)xxx[i]=(-xxx[i]+a[i]+P)%P;
	xxx[0]=(xxx[0]+1)%P;
	Comb(b,xxx,b,n);
}
void exp(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_e(a,b,n);
}
int A[N],xxx[N],B[N],n,k;
int main()
{
	n=read();k=read();
	reg int i;
	for(i=0;i<n;++i) A[i]=read();
	ln(A,xxx,n);
	for(i=0;i<n;++i) xxx[i]=Mul(xxx[i],k);
	exp(xxx,B,n);
	for(i=0;i<n;++i) printf("%d ",B[i]);
	return 0;
}

多项式开根

问题描述

给出一个(n-1)次的多项式(A(x))，求(B^2(x)≡A(x) mod x^{n})，(a_0=1)，然后(mod 998244353)

方法

首先

[F(B(x))=B^2(x)-A(x)=0 ]

考虑牛顿迭代

[egin{equation} egin{split} B(x)&leftarrow B_0(x)-frac{F(B_0(x))}{F'(B_0(x))}\ &=B_0(x)-frac{{B^2_0}(x)-A(x)}{2B(x)}\ &=2^{-1}(B_0(x)+A(x)B^{-1}(x)) end{split} end{equation} ]

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define reg register
const int N=1<<18|5,P=998244353,g=3,invg=332748118,inv2=499122177;
inline int read()
{
	ll x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x*10ll+ch-'0')%P;ch=getchar();}
	return x*f;
}
#define Mul(x,y) (1ll*(x)*(y)%P)
#define ri reg int i
int _[N],pos[N];
int fpow(int x,int m){int r=1;for(;m;m>>=1,x=Mul(x,x))if(m&1)r=Mul(r,x);return r;}
void der(int *a,int *b,int n){for(ri=0;i<n-1;++i) b[i]=Mul(i+1,a[i+1]);b[n-1]=0;}
void inte(int *a,int *b,int n)
{
	ri;for(_[1]=i=1;i<n;++i,_[i]=Mul(_[P%i],P-P/i));
	for(i=n-1;i;--i)b[i]=Mul(a[i-1],_[i]);b[0]=0;
}
void init(int len){for(ri=0;i<len;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(len>>1));}
void NTT(int *a,int n,int type)
{
	register int i,j,k,w,wn,X,Y;
	for(i=0;i<n;++i) if(i<pos[i]) swap(a[i],a[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		wn=fpow(type>0?g:invg,(P-1)/(i<<1));
		for(j=0;j<n;j+=i<<1)for(w=1,k=0;k<i;++k,w=Mul(w,wn))
		{
			X=a[j+k],Y=Mul(w,a[j+k+i]);
			a[j+k]=(X+Y)%P;a[j+k+i]=(X-Y+P)%P;
		}
	}
	reg int invn=fpow(n,P-2);
	if(type==-1)for(i=0;i<n;++i)a[i]=Mul(a[i],invn);
}
void Comb(int *a,int *b,int *c,int n)
{
	static int A[N],B[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(A,a,sizeof(int[n]));memcpy(B,b,sizeof(int[n]));
	memset(A+n,0,sizeof(int[len-n]));memset(B+n,0,sizeof(int[len-n]));
	ri;init(len);NTT(A,len,1);NTT(B,len,1);
	for(i=0;i<len;++i) c[i]=Mul(A[i],B[i]);NTT(c,len,-1);
	memset(c+n,0,sizeof(int[len-n]));
}
void _I(int *A,int *b,int n)
{
	if(n==1) return(void)(b[0]=fpow(A[0],P-2));
	int t=(n+1)>>1;_I(A,b,t);
	static int a[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(a,A,sizeof(int[n]));memset(a+n,0,sizeof(int[len-n]));
	memset(b+t,0,sizeof(int[len-t]));init(len);
	ri;NTT(a,len,1);NTT(b,len,1);
	for(i=0;i<len;++i)b[i]=(Mul(2ll,b[i])-Mul(a[i],Mul(b[i],b[i]))+P)%P;
	NTT(b,len,-1);memset(b+n,0,sizeof(int[len-n]));
}
void Inv(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	_I(a,b,n);
}
void _S(int *A,int *b,int n)
{
	if(n==1) return(void)(b[0]=1);
	int t=(n+1)>>1;_S(A,b,t);
	static int inv_b[N],a[N],len;
	for(len=1;len<(n<<1);len<<=1);
	memcpy(a,A,sizeof(int[n]));memset(a+n,0,sizeof(int[len-n]));
	Inv(b,inv_b,n);Comb(inv_b,a,a,n);
	for(ri=0;i<n;++i) b[i]=Mul((b[i]+a[i])%P,inv2);
    memset(b+n,0,sizeof(int[len-n]));
}
void Sqrt(int *A,int *b,int n)
{
	static int a[N];
	memcpy(a,A,sizeof(int[n]));
	memset(b,0,sizeof(int[n]));
	_S(a,b,n);
}
int A[N],n,k;
int main()
{
	n=read();
	reg int i;
	for(i=0;i<n;++i) A[i]=read();
	Sqrt(A,A,n);for(i=0;i<n;++i) printf("%d ",A[i]);
	return 0;
}

Blog来自PaperCloud，未经允许，请勿转载，TKS！

查看全文

相关阅读:
从零开始学android开发-通过WebService获取今日天气情况
 android常见错误-E/AndroidRuntime(13678): java.lang.NoClassDefFoundError:
java 使用相对路径读取文件
 冒泡排序
 快速排序
 为什么使用抽象类？有什么好处？
为什么用抽象类，接口
 String.valueOf()
Python 资源
 文本相似度-BM25算法

原文地址：https://www.cnblogs.com/PaperCloud/p/10646165.html