题目链接: https://www.luogu.org/problemnew/show/P4238
题意: 给定 n n n次多项式 A ( x ) A(x) A(x), 求 n n n次多项式 B ( x ) B(x) B(x)满足 B ( x ) A ( x ) ≡ 1 ( m o d    x n ) B(x)A(x)\equiv 1(\mod x^n) B(x)A(x)≡1(modxn)
题解:
DFT,每个数对
998244353
998244353
998244353求逆元。IDFT回来。
发现,错了。
为什么呢?
因为要对
x
n
x^n
xn取模。
例如,
1
−
x
1-x
1−x在模
x
2
x^2
x2意义下的逆元是
1
+
x
1+x
1+x, 但是在实际上逆元是
1
+
x
+
x
2
+
x
3
+
.
.
.
1+x+x^2+x^3+...
1+x+x2+x3+..., 是无穷和式。
所以此路不通。
考虑求解多项式问题的常用方法——分治法。
设已求
B
0
(
x
)
B_0(x)
B0(x)满足
B
0
(
x
)
A
(
x
)
≡
1
(
m
o
d
  
x
n
)
B_0(x)A(x)\equiv 1(\mod x^n)
B0(x)A(x)≡1(modxn), 现要求
B
(
x
)
B(x)
B(x)满足
B
(
x
)
A
(
x
)
≡
1
(
m
o
d
  
x
2
n
)
B(x)A(x)\equiv 1(\mod x^{2n})
B(x)A(x)≡1(modx2n)
显然有
B
(
x
)
−
B
0
(
x
)
≡
0
(
m
o
d
  
x
n
)
B(x)-B_0(x)\equiv 0(\mod x^n)
B(x)−B0(x)≡0(modxn)
为了凑出
x
2
n
x^{2n}
x2n两边平方得
B
2
(
x
)
−
2
B
0
(
x
)
B
(
x
)
+
B
0
2
(
x
)
≡
0
(
m
o
d
  
x
2
n
)
B^2(x)-2B_0(x)B(x)+B_0^2(x)\equiv 0(\mod x^{2n})
B2(x)−2B0(x)B(x)+B02(x)≡0(modx2n)
如何求
B
B
B呢?因为
A
(
x
)
B
(
x
)
≡
0
(
m
o
d
  
x
2
n
)
A(x)B(x)\equiv 0(\mod x^{2n})
A(x)B(x)≡0(modx2n), 我们将式子两边同乘
A
(
x
)
A(x)
A(x)
A
(
x
)
B
(
x
)
B
(
x
)
−
2
B
0
(
x
)
A
(
x
)
B
(
x
)
+
A
(
x
)
B
0
2
(
x
)
≡
0
(
m
o
d
  
x
2
n
)
A(x)B(x)B(x)-2B_0(x)A(x)B(x)+A(x)B_0^2(x)\equiv 0(\mod x^{2n})
A(x)B(x)B(x)−2B0(x)A(x)B(x)+A(x)B02(x)≡0(modx2n)
B
(
x
)
≡
2
B
0
(
x
)
−
A
(
x
)
B
0
2
(
x
)
(
m
o
d
  
x
2
n
)
B(x)\equiv 2B_0(x)-A(x)B_0^2(x) (\mod x^{2n})
B(x)≡2B0(x)−A(x)B02(x)(modx2n)
右边的式子FFT计算即可。
时间复杂度?
T
(
n
)
=
T
(
n
2
)
+
O
(
n
log
n
)
T(n)=T(\frac{n}{2})+O(n\log n)
T(n)=T(2n)+O(nlogn)
T
(
n
)
=
O
(
n
log
n
)
T(n)=O(n\log n)
T(n)=O(nlogn).
常数?首先隐藏在递归复杂度背后有一个
2
2
2倍常数。
然后我们把两个多项式相乘需要
3
3
3次FFT, 三个就要
6
6
6次。
因此常数为
12
12
12倍。
如何优化?
I
D
F
T
(
D
F
T
(
I
D
F
T
(
D
F
T
(
A
)
×
D
F
T
(
B
0
)
)
)
×
D
F
T
(
B
0
)
)
IDFT(DFT(IDFT(DFT(A)\times DFT(B_0)))\times DFT(B_0))
IDFT(DFT(IDFT(DFT(A)×DFT(B0)))×DFT(B0))
变成
I
D
F
T
(
D
F
T
(
A
)
×
D
F
T
2
(
B
0
)
IDFT(DFT(A)\times DFT^2(B_0)
IDFT(DFT(A)×DFT2(B0)
3
3
3次即可!
常数变为
6
6
6倍。
UPD: 关于这里的常数问题: 因为我递归里DFT的范围是
2
n
2n
2n,最终的复杂度是
T
(
2
n
)
=
6
(
2
n
)
log
(
2
n
)
T(2n) = 6(2n)\log (2n)
T(2n)=6(2n)log(2n), 因此我认为应为
12
12
12倍常数。
空间?空间复杂度
O
(
n
)
O(n)
O(n), 但是要开
4
d
4d
4d的数组,其中
d
d
d是
>
n
>n
>n的最小的
2
2
2的幂。
代码
因为FFT数组清零等原因代码(对我来说)很难写。
贴一下我刚刚写的版本吧,还算是比较简单。
(话说怎么CSDN突然傲娇了啊。。缩进变成1格??)
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define llong long long
#define ldouble long double
#define uint unsigned int
#define ullong unsigned long long
#define udouble unsigned double
#define uldouble unsigned long double
#define modinc(x) {if(x>=P) x-=P;}
#define pii pair<int,int>
#define piii pair<pair<int,int>,int>
#define piiii pair<pair<int,int>,pair<int,int> >
#define pli pair<llong,int>
#define pll pair<llong,llong>
#define Memset(a,x) {memset(a,x,sizeof(a));}
using namespace std;
const int N = 1<<19;
const int P = 998244353;
const int LGN = 19;
const int G = 3;
llong a[N+3];
llong b[N+3];
llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3],tmp5[N+3],tmp6[N+3];
int id[N+2];
int n;
void initid(int _len)
{
id[0] = 0;
for(int i=1; i<(1<<_len); i++) id[i] = (id[i>>1]>>1)|((i&1)<<(_len-1));
}
llong quickpow(llong x,llong y)
{
llong cur = x,ret = 1ll;
for(int i=0; y; i++)
{
if(y&(1ll<<i))
{
y-=(1ll<<i); ret = ret*cur%P;
}
cur = cur*cur%P;
}
return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}
void ntt(int dgr,int coe,llong poly[],llong ret[])
{
int len = 0; for(int i=0; i<=LGN; i++) if((1<<i)==dgr) {len = i; break;}
initid(len); for(int i=0; i<dgr; i++) ret[i] = 0ll;
for(int i=0; i<dgr; i++) ret[i] = poly[i];
for(int i=0; i<dgr; i++) if(i<id[i]) swap(ret[i],ret[id[i]]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
llong tmp = quickpow(G,(P-1)/(i<<1));
if(coe==-1) tmp = mulinv(tmp);
for(int j=0; j<dgr; j+=(i<<1))
{
llong expn = 1ll;
for(int k=0; k<i; k++)
{
llong x = ret[j+k],y = (expn*ret[j+i+k])%P;
ret[j+k] = x+y; modinc(ret[j+k]);
ret[j+i+k] = x-y+P; modinc(ret[j+i+k]);
expn = (expn*tmp)%P;
}
}
}
}
void polyinv(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr; i++) ret[i] = 0ll;
ret[0] = mulinv(poly[0]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
for(int j=0; j<(i<<2); j++) tmp1[j] = j<i ? ret[j] : 0ll;
for(int j=0; j<(i<<2); j++) tmp2[j] = j<(i<<1) ? poly[j] : 0ll;
ntt((i<<2),1,tmp1,tmp3); ntt((i<<2),1,tmp2,tmp4);
for(int j=0; j<(i<<2); j++) tmp5[j] = tmp3[j]*tmp3[j]%P*tmp4[j]%P;
ntt((i<<2),-1,tmp5,tmp6); llong tmp = mulinv(i<<2);
for(int j=0; j<(i<<2); j++) tmp6[j] = tmp6[j]*tmp%P;
for(int j=0; j<(i<<1); j++) ret[j] = (tmp1[j]+tmp1[j]-tmp6[j]+P)%P;
}
}
int main()
{
scanf("%d",&n); int dgr = 1; while(dgr<=n) dgr<<=1;
for(int i=0; i<n; i++) scanf("%lld",&a[i]);
polyinv(dgr,a,b);
for(int i=0; i<n; i++) printf("%lld ",b[i]);
return 0;
}