<更新提示>
<第一次更新> 在教练的要求下开始学习多项式算法了，不过因为不太会积分和求导先把多项式牛顿迭代，多项式指数函数，多项式幂函数，多项式快速幂等内容咕掉了，于是这一篇博客就是其他基础多项式内容的总结。
<正文>
Fast Fourier Transform
(FFT)，快速傅里叶变换，可以在(O(nlog_2n))的时间内计算多项式乘法。
先回忆一下(FFT)的思路，首先是多项式的系数表达法，如果直接计算的话，时间复杂度是(O(n^2))的，可以用分治算法优化到(O(n^{log_23}))，通常没有其他更优的算法了。
然后就是考虑多项式的另一个表达方式，点值表达法，这样就可以实现(O(n))乘法。
但是使用普通的求值方法和拉格朗日插值法实现系数表达法和点值表达法的转换是(O(n^2))的，毫无意义。
后来一个分治算法出现了，我们可以在求值的时候代入(n)个单位根，根据单位根的性质，我们可以只计算一半的值得到另一半，从而递归下去分治解决，实现时间复杂度(O(nlog_2n))。
说到单位根，就必须提一下有关复数的基础知识，可以看这篇博客。
然后就是单位根的性质了：
(1.) (omega_{n}^k=cos kfrac{2pi}{n}+isin kfrac{2pi}{n})
(2.) (omega_{2n}^{2k}=omega_n^k)
(3.) (omega_{n}^{k+frac{n}{2}}=-omega_{n}^{k})
(4.) (k ot=0)时，(sum_{i=0}^{n-1}(omega_{n}^{k})^i=0)，(k=0)时，(sum_{i=0}^{n-1}(omega_{n}^{k})^i=n)。
然后把多项式下标按奇偶性分类，就能分治了。
通过构造可知，我们可以用同样的方式也可以把点值表达法转换回系数表达法，需要除以一个序列长度的大小。
最后的最后就是把递归的(FFT)用迭代实现，就是一点二进制的小(Trick)，叫蝴蝶定理，找规律的事。
可以看这篇博客。
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = (1<<22) + 20;
namespace Mathcalc
{
const double PI = acos(-1.0) , eps = 5e-3;
struct Complex
{
double x,y;
Complex ( double _x = 0 , double _y = 0 ) { x = _x , y = _y; }
friend Complex operator + (Complex a,Complex b) { return Complex( a.x + b.x , a.y + b.y ); }
friend Complex operator - (Complex a,Complex b) { return Complex( a.x - b.x , a.y - b.y ); }
friend Complex operator * (Complex a,Complex b) { return Complex( a.x*b.x - a.y*b.y , a.x*b.y + a.y*b.x ); }
};
}
using namespace Mathcalc;
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
int n,m,lenth,cnt,rev[N];
Complex a[N],b[N],c[N];
inline void input(void)
{
n = read() , m = read();
for (int i=0;i<=n;i++) a[i].x = read();
for (int i=0;i<=m;i++) b[i].x = read();
}
inline void calcrev(void)
{
lenth = 1;
while ( lenth <= n + m ) lenth <<= 1 , cnt++;
for (int i=0;i<lenth;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) );
}
inline void FFT(Complex *p,int op)
{
for (int i=0;i<lenth;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<lenth;i<<=1)
{
Complex w ( cos(PI/i) , op * sin(PI/i) );
for (int j=0;j<lenth;j+=i<<1)
{
Complex omega ( 1 , 0 );
for (int k=0;k<i;k++)
{
Complex x = p[j+k] , y = omega * p[i+j+k];
p[j+k] = x + y , p[i+j+k] = x - y;
omega = omega * w;
}
}
}
if ( op == -1 )
for (int i=0;i<lenth;i++)
p[i].x /= lenth;
}
int main(void)
{
input();
calcrev();
FFT( a , 1 ) , FFT( b , 1 );
for (int i=0;i<lenth;i++)
c[i] = a[i] * b[i];
FFT( c , -1 );
for (int i=0;i<=n+m;i++)
printf("%d ",int(c[i].x + eps));
return 0;
}
Number Theoretic Transform
(NTT)，快速数论变换，可以在(O(nlog_2n))的时间内计算模意义下的多项式乘法。
思路和(FFT)是一模一样的，关键在于我们在模意义下找到了单位根的替代品：原根，满足单位根的所有性质。
然后按照(FFT)的方式来就可以了。
但是(NTT)的模数是有限定的，可以看这个。
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = (1<<22)+20;
namespace Mathcalc
{
const int Mod = 998244353 , G = 3 , INVG = 332748118;
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Mul(int &a,int b) { a = mul( a , b ); }
inline int quickpow(int a,int b) { int res = 1; for (;b;Mul(a,a),b>>=1) if ( b & 1 ) Mul( res , a ); return res; }
inline int inv(int a) { return quickpow( a , Mod-2 ); }
};
using namespace Mathcalc;
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
int n,m,lenth,cnt,rev[N];
int a[N],b[N],c[N];
inline void input(void)
{
n = read() - 1 , m = read() - 1;
for (int i=0;i<=n;i++) a[i] = read();
for (int i=0;i<=m;i++) b[i] = read();
}
inline void calcrev(void)
{
lenth = 1;
while ( lenth <= n + m ) lenth <<= 1 , cnt++;
for (int i=0;i<lenth;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) );
}
inline void NTT(int *p,int op)
{
for (int i=0;i<lenth;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<lenth;i<<=1)
{
int w = quickpow( op == 1 ? G : INVG , (Mod-1) / (i<<1) );
for (int j=0;j<lenth;j+=i<<1)
{
int omega = 1;
for (int k=0;k<i;k++)
{
int x = p[j+k] , y = mul( omega , p[i+j+k] );
p[j+k] = add( x , y ) , p[i+j+k] = sub( x , y );
Mul( omega , w );
}
}
}
if ( op == -1 )
{
int val = inv( lenth );
for (int i=0;i<lenth;i++)
Mul( p[i] , val );
}
}
int main(void)
{
input();
calcrev();
NTT( a , 1 ) , NTT( b , 1 );
for (int i=0;i<lenth;i++)
c[i] = mul( a[i] , b[i] );
NTT( c , -1 );
for (int i=0;i<=n+m;i++)
printf("%d ",c[i]);
return 0;
}
分治FFT
我们知道，(FFT/NTT)是可以求形容(c_k=sum_{i+j=k}a_i imes b_j)的卷积的，但是如果是长这样的卷积呢？
[f_k=sum_{i+j=k}f_i imes g_j
]
仿佛套一个(cdq)分治算贡献就可以了。
没错，事实上也就是这样，中间的卷积还是要(NTT)，时间复杂度(O(nlog^2n))。
(Cocd:)
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5+20 , SIZE = (1<<18)+20;
const int Mod = 998244353 , G = 3 , INVG = 332748118;
int n,cnt,lenth,rev[SIZE];
long long f[N],g[N],a[SIZE],b[SIZE],c[SIZE];
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
inline void input(void)
{
n = read() - 1 , f[0] = 1;
for (int i=1;i<=n;i++)
g[i] = read();
}
inline long long quickpow(long long a,long long b)
{
long long res = 1;
while ( b )
{
if ( 1 & b ) res = res * a % Mod;
b >>= 1 , a = a * a % Mod;
}
return res;
}
inline void NTT(long long *p,int op)
{
for (int i=0;i<lenth;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<lenth;i<<=1)
{
long long w = quickpow( op == 1 ? G : INVG , (Mod-1) / (i<<1) );
for (int j=0;j<lenth;j+=i<<1)
{
long long omega = 1;
for (int k=0;k<i;k++)
{
long long x = p[j+k] , y = omega * p[i+j+k] % Mod;
p[j+k] = ( x + y ) % Mod;
p[i+j+k] = ( ( x - y ) % Mod + Mod ) % Mod;
omega = omega * w % Mod;
}
}
}
}
inline void cdq(int l,int r)
{
if ( l == r ) return void();
int mid = l+r >> 1 , len = (r-l-1) + (mid-l);
cdq( l , mid );
lenth = 1 , cnt = 0;
while ( lenth <= len ) lenth <<= 1 , cnt++;
for (int i=0;i<lenth;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) ),
a[i] = b[i] = c[i] = 0;
for (int i=l;i<=mid;i++) a[i-l] = f[i];
for (int i=1;i<=r-l;i++) b[i-1] = g[i];
NTT( a , 1 ) , NTT( b , 1 );
for (int i=0;i<lenth;i++) c[i] = a[i] * b[i] % Mod;
NTT( c , -1 );
long long inv = quickpow( lenth , Mod-2 );
for (int i=0;i<=len;i++)
c[i] = c[i] * inv % Mod;
for (int i=mid+1;i<=r;i++)
f[i] = ( f[i] + c[i-l-1] ) % Mod;
cdq( mid + 1 , r );
}
int main(void)
{
input();
cdq( 0 , n );
for (int i=0;i<=n;i++)
printf("%lld ",f[i]);
puts("");
return 0;
}
Fast Walsh-Hadamard Transform
(FWT)，快速沃尔什变换，用于计算位运算卷积。
形式是这样的(c_k=sum_{i op j=k}a_i imes b_j)，(op)是(and,or,xor)中的一个位运算。
思路和(FWT)是一样的，考虑先对系数表达法做一点变换，然后点乘(O(n))解决问题，再逆变换回来就可以了。
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = 20;
namespace Mathcalc
{
const int Mod = 998244353 , INV2 = 499122177;
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Add(int &a,int b) { a = add( a , b ); }
inline void Sub(int &a,int b) { a = sub( a , b ); }
inline void Mul(int &a,int b) { a = mul( a , b ); }
};
using namespace Mathcalc;
int n; int a[1<<N],b[1<<N],c[1<<N];
inline void input(void)
{
scanf("%d",&n);
for (int i=0;i<1<<n;i++)
scanf("%d",&a[i]);
for (int i=0;i<1<<n;i++)
scanf("%d",&b[i]);
}
inline void FWTor(int *p,int op)
{
for (int i=1;i<1<<n;i<<=1)
for (int j=0;j<1<<n;j+=i<<1)
for (int k=0;k<i;k++)
if ( op == 1 ) Add( p[i+j+k] , p[j+k] );
else Sub( p[i+j+k] , p[j+k] );
}
inline void FWTand(int *p,int op)
{
for (int i=1;i<1<<n;i<<=1)
for (int j=0;j<1<<n;j+=i<<1)
for (int k=0;k<i;k++)
if ( op == 1 ) Add( p[j+k] , p[i+j+k] );
else Sub( p[j+k] , p[i+j+k] );
}
inline void FWTxor(int *p,int op)
{
for (int i=1;i<1<<n;i<<=1)
for (int j=0;j<1<<n;j+=i<<1)
for (int k=0;k<i;k++)
{
int x = p[j+k] , y = p[i+j+k];
p[j+k] = mul( op , add( x , y ) );
p[i+j+k] = mul( op , sub( x , y ) );
}
}
inline void solve(void)
{
FWTor( a , 1 ) , FWTor( b , 1 );
for (int i=0;i<1<<n;i++) c[i] = mul( a[i] , b[i] );
FWTor( c , -1 );
for (int i=0;i<1<<n;i++) printf("%d ",c[i]);
puts("") , FWTor( a , -1 ) , FWTor( b , -1 );
FWTand( a , 1 ) , FWTand( b , 1 );
for (int i=0;i<1<<n;i++) c[i] = mul( a[i] , b[i] );
FWTand( c , -1 );
for (int i=0;i<1<<n;i++) printf("%d ",c[i]);
puts("") , FWTand( a , -1 ) , FWTand( b , -1 );
FWTxor( a , 1 ) , FWTxor( b , 1 );
for (int i=0;i<1<<n;i++) c[i] = mul( a[i] , b[i] );
FWTxor( c , INV2 );
for (int i=0;i<1<<n;i++) printf("%d ",c[i]);
puts("") , FWTxor( a , INV2 ) , FWTxor( b , INV2 );
}
int main(void)
{
input();
solve();
return 0;
}
多项式求逆
给定多项式(F(x))，求多项式(G(x))，满足(F(x)G(x)equiv1pmod {x^n})。
当多项式只有一个常数项的时候，求逆元即可。
反之，考虑一个倍增算法：
假设(F imes G'equiv 1pmod{x^{frac{n}{2}}})，我们要求(F imes Gequiv1pmod{x^n})。
[F imes(G-G')equiv 0pmod {x^{frac{n}{2}}}\ \G-G'equiv 0pmod {x^{frac{n}{2}}}\ \(G-G')^2equiv 0pmod {x^n}\ \G^2-2GG'+G'^2equiv 0pmod{x^n}\ \F imes(G^2-2GG'+G'^2)equiv 0pmod{x^n}\ \G-2G'+FG'^2equiv 0pmod{x^n}\ \Gequiv2G'+FG'^2pmod{x^n}
]
结合(NTT)，就可以倍增计算了，时间复杂度(O(nlog_2n))。
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = (1<<20)+20;
namespace Mathcalc
{
const int Mod = 998244353 , G = 3 , INVG = 332748118;
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Add(int &a,int b) { a = add( a , b ); }
inline void Sub(int &a,int b) { a = sub( a , b ); }
inline void Mul(int &a,int b) { a = mul( a , b ); }
inline int quickpow(int a,int b) { int res = 1; for (;b;Mul(a,a),b>>=1) if ( b & 1 ) Mul( res , a ); return res; }
inline int inv(int a) { return quickpow( a , Mod-2 ); }
};
using namespace Mathcalc;
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
int temp[N],rev[N];
inline void calcrev(int len,int cnt)
{
for (int i=0;i<len;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) );
}
inline void NTT(int *p,int op,int len)
{
for (int i=0;i<len;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<len;i<<=1)
{
int w = quickpow( op == 1 ? G : INVG , (Mod-1) / (i<<1) );
for (int j=0;j<len;j+=i<<1)
{
int omega = 1;
for (int k=0;k<i;k++)
{
int x = p[j+k] , y = mul( omega , p[i+j+k] );
p[j+k] = add( x , y ) , p[i+j+k] = sub( x , y );
Mul( omega , w );
}
}
}
if ( op == -1 )
{
int val = inv( len );
for (int i=0;i<len;i++)
Mul( p[i] , val );
}
}
inline void Polyinv(int *a,int *b,int n)
{
int p = 1; while ( p < n ) p <<= 1;
for (int i=0;i<p;i++) b[i] = 0;
b[0] = inv( a[0] );
for (int len=1,cnt=0;len<=p;len<<=1,cnt++)
{
calcrev( len<<1 , cnt+1 );
for (int i=0;i<len;i++) temp[i] = a[i];
for (int i=len;i<len<<1;i++) temp[i] = 0;
NTT( b , 1 , len<<1 ) , NTT( temp , 1 , len<<1 );
for (int i=0;i<len<<1;i++)
Mul( b[i] , sub( 2 , mul( temp[i] , b[i] ) ) );
NTT( b , -1 , len<<1 );
for (int i=len;i<len<<1;i++) b[i] = 0;
}
for (int i=n;i<p;i++) b[i] = 0;
}
int n,m,a[N],b[N];
inline void input(void)
{
n = read() , m = read();
for (int i=0;i<n;i++) a[i] = read();
}
int main(void)
{
input();
Polyinv( a , b , m );
for (int i=0;i<m;i++)
printf("%d ",b[i]);
return 0;
}
多项式除法
给定(n)次多项式(F(x))和(m)次多项式(G(x))，求多项式(Q(x),R(x))，满足：
(1.) (Q(x))为(n-m)次多项式，(R(x))的次数小于(m)。
(2.) (F(x)=G(x) imes Q(x)+R(x))。
考虑一个操作(rev)，对于多项式(F(x))，(F^{rev}(x)=x^nF(frac{1}{x}))，就是系数翻转。
然后考虑推一波式子：
[F(x)=Q(x) imes G(x)+R(x)
\ \
Fleft (frac{1}{x}
ight )=Gleft (frac{1}{x}
ight ) imes Qleft (frac{1}{x}
ight )+Rleft (frac{1}{x}
ight )
\ \
x^nFleft (frac{1}{x}
ight )=x^mGleft (frac{1}{x}
ight ) imes x^{n-m} Qleft (frac{1}{x}
ight )+x^{n-m+1}x^{m-1}Rleft (frac{1}{x}
ight )
\ \
F^{rev}(x)=Q^{rev}(x) imes G^{rev}(x)+x^{n-m+1}R^{rev}(x)
\ \
F^{rev}(x)=Q^{rev}(x) imes G^{rev}(x)pmod{x^{n-m+1}}
\ \
Q^{rev}(x)=F^{rev}(x) imes (G^{rev})^{-1}(x)pmod{x^{n-m+1}}
]
考虑对(G^{rev})做一个多项式求逆，然后再用(NTT)把它和(F^{rev})乘起来就是(Q_{rev})了。
然后(R)直接算就可以了：
[R(x)=F(x)-Q(x) imes G(x)
]
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = (1<<20)+20;
namespace Mathcalc
{
const int Mod = 998244353 , G = 3 , INVG = 332748118;
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Add(int &a,int b) { a = add( a , b ); }
inline void Sub(int &a,int b) { a = sub( a , b ); }
inline void Mul(int &a,int b) { a = mul( a , b ); }
inline int quickpow(int a,int b) { int res = 1; for (;b;Mul(a,a),b>>=1) if ( b & 1 ) Mul( res , a ); return res; }
inline int inv(int a) { return quickpow( a , Mod-2 ); }
};
using namespace Mathcalc;
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
int temp[N],t1[N],t2[N],rev[N];
inline void calcrev(int len,int cnt)
{
for (int i=0;i<len;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) );
}
inline void NTT(int *p,int op,int len)
{
for (int i=0;i<len;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<len;i<<=1)
{
int w = quickpow( op == 1 ? G : INVG , (Mod-1) / (i<<1) );
for (int j=0;j<len;j+=i<<1)
{
int omega = 1;
for (int k=0;k<i;k++)
{
int x = p[j+k] , y = mul( omega , p[i+j+k] );
p[j+k] = add( x , y ) , p[i+j+k] = sub( x , y );
Mul( omega , w );
}
}
}
if ( op == -1 )
{
int val = inv( len );
for (int i=0;i<len;i++)
Mul( p[i] , val );
}
}
inline void Polyinv(int *a,int *b,int n)
{
int p = 1; while ( p < n ) p <<= 1;
for (int i=0;i<p;i++) b[i] = 0;
b[0] = inv( a[0] );
for (int len=1,cnt=0;len<=p;len<<=1,cnt++)
{
calcrev( len<<1 , cnt+1 );
for (int i=0;i<len;i++) temp[i] = a[i];
for (int i=len;i<len<<1;i++) temp[i] = 0;
NTT( b , 1 , len<<1 ) , NTT( temp , 1 , len<<1 );
for (int i=0;i<len<<1;i++)
Mul( b[i] , sub( 2 , mul( temp[i] , b[i] ) ) );
NTT( b , -1 , len<<1 );
for (int i=len;i<len<<1;i++) b[i] = 0;
}
for (int i=n;i<p;i++) b[i] = 0;
}
inline void Polymul(int *a,int *b,int *c,int len)
{
NTT( a , 1 , len ) , NTT( b , 1 , len );
for (int i=0;i<len;i++)
c[i] = mul( a[i] , b[i] );
NTT( c , -1 , len );
}
inline void Polydiv(int *f,int *g,int *q,int *r,int n,int m)
{
int len = 1 , cnt = 0;
while ( len < 2*(n-m+1) ) len <<= 1 , cnt++;
for (int i=0;i<n-m+1;i++) t2[i] = g[m-i-1];
for (int i=n-m+1;i<len;i++) t2[i] = 0;
Polyinv( t2 , t1 , n-m+1 );
for (int i=0;i<n-m+1;i++) t2[i] = f[n-i-1];
for (int i=n-m+1;i<len;i++) t2[i] = 0;
calcrev( len , cnt ) , Polymul( t1 , t2 , q , len );
reverse( q , q+n-m+1 );
len = 1 , cnt = 0;
while ( len < n ) len <<= 1 , cnt++;
for (int i=0;i<n-m+1;i++) t1[i] = q[i];
for (int i=n-m+1;i<len;i++) t1[i] = 0;
for (int i=0;i<m;i++) t2[i] = g[i];
for (int i=m;i<len;i++) t2[i] = 0;
calcrev( len , cnt ) , Polymul( t1 , t2 , r , len );
for (int i=0;i<len;i++)
r[i] = sub( f[i] , r[i] );
}
int n,m,f[N],g[N],q[N],r[N];
int main(void)
{
n = read() , m = read();
for (int i=0;i<n;i++) f[i] = read();
for (int i=0;i<m;i++) g[i] = read();
Polydiv( f , g , q , r , n , m );
for (int i=0;i<n-m+1;i++)
printf("%d ",q[i]);
puts("");
for (int i=0;i<m-1;i++)
printf("%d ",r[i]);
puts("");
return 0;
}
多项式开根
给定一个(n-1)次多项式(A(x))，求一个在(mod x^n)意义下的多项式(B(x))，使得(B^2(x) equiv A(x) (mod x^n))。
当(A)只有一个常数项时，(B)就是(A_0)的二次剩余，可以用(cipolla)求，当然这道题不用。
然后考虑倍增：
假设(Aequiv B'^2pmod{x^{frac{n}{2}}})，我们要求(Aequiv B^2pmod{x^n})。
[B^2-B'^2equiv 0pmod{x^{frac{n}{2}}}\ \B^4+B'^4-2B^2B'^2equiv0pmod{x^n}\ \ B^2+B'^2equiv2BB'pmod{x^n}\ \ A+B'^2equiv2BB'pmod{x^n}\ \ Bequiv frac{A+B'^2}{2B'}equiv frac{1}{2}left (frac{A}{B'}+B'^2
ight)pmod{x^n}
]
多项式求一下逆，然后倍增即可，时间复杂度(O(nlog_2n))。
(Code:)
#include <bits/stdc++.h>
using namespace std;
const int N = (1<<20)+20;
namespace Mathcalc
{
const int Mod = 998244353 , G = 3 , INVG = 332748118 , INV2 = 499122177;
inline int add(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int sub(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Add(int &a,int b) { a = add( a , b ); }
inline void Sub(int &a,int b) { a = sub( a , b ); }
inline void Mul(int &a,int b) { a = mul( a , b ); }
inline int quickpow(int a,int b) { int res = 1; for (;b;Mul(a,a),b>>=1) if ( b & 1 ) Mul( res , a ); return res; }
inline int inv(int a) { return quickpow( a , Mod-2 ); }
};
using namespace Mathcalc;
inline int read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch=='-' , ch = getchar();
while ( isdigit(ch) ) x = x*10 + ch-48 , ch = getchar();
return w ? -x : x;
}
int temp[N],t1[N],t2[N],rev[N];
inline void calcrev(int len,int cnt)
{
for (int i=0;i<len;i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << (cnt-1) );
}
inline void NTT(int *p,int op,int len)
{
for (int i=0;i<len;i++)
if ( i < rev[i] ) swap( p[i] , p[rev[i]] );
for (int i=1;i<len;i<<=1)
{
int w = quickpow( op == 1 ? G : INVG , (Mod-1) / (i<<1) );
for (int j=0;j<len;j+=i<<1)
{
int omega = 1;
for (int k=0;k<i;k++)
{
int x = p[j+k] , y = mul( omega , p[i+j+k] );
p[j+k] = add( x , y ) , p[i+j+k] = sub( x , y );
Mul( omega , w );
}
}
}
if ( op == -1 )
{
int val = inv( len );
for (int i=0;i<len;i++)
Mul( p[i] , val );
}
}
inline void Polyinv(int *a,int *b,int n)
{
int p = 1; while ( p < n ) p <<= 1;
for (int i=0;i<p;i++) b[i] = 0;
b[0] = inv( a[0] );
for (int len=1,cnt=0;len<=p;len<<=1,cnt++)
{
calcrev( len<<1 , cnt+1 );
for (int i=0;i<len;i++) temp[i] = a[i];
for (int i=len;i<len<<1;i++) temp[i] = 0;
NTT( b , 1 , len<<1 ) , NTT( temp , 1 , len<<1 );
for (int i=0;i<len<<1;i++)
Mul( b[i] , sub( 2 , mul( temp[i] , b[i] ) ) );
NTT( b , -1 , len<<1 );
for (int i=len;i<len<<1;i++) b[i] = 0;
}
for (int i=n;i<p;i++) b[i] = 0;
}
inline void Polymul(int *a,int *b,int *c,int len)
{
NTT( a , 1 , len ) , NTT( b , 1 , len );
for (int i=0;i<len;i++)
c[i] = mul( a[i] , b[i] );
NTT( c , -1 , len );
}
inline void Polysqrt(int *a,int *b,int n)
{
int p = 1; while ( p < n ) p <<= 1;
for (int i=0;i<p;i++) b[i] = 0;
b[0] = 1;
for (int len=1,cnt=0;len<=p;len<<=1,cnt++)
{
calcrev( len<<1 , cnt+1 );
Polyinv( b , t1 , len );
for (int i=0;i<len;i++) t2[i] = a[i];
for (int i=len;i<len<<1;i++) t2[i] = 0;
Polymul( t1 , t2 , t1 , len<<1 );
for (int i=0;i<len;i++)
b[i] = mul( add( b[i] , t1[i] ) , INV2 );
for (int i=len;i<len<<1;i++) b[i] = 0;
}
for (int i=n;i<p;i++) b[i] = 0;
}
int n,a[N],b[N];
int main(void)
{
n = read();
for (int i=0;i<n;i++) a[i] = read();
Polysqrt( a , b , n );
for (int i=0;i<n;i++)
printf("%d ",b[i]);
return 0;
}
<后记>