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！