由于复数域的神奇性质,我们在FFT的时候可以将计算
C
(
x
)
=
A
(
x
)
B
(
x
)
C(x)=A(x)B(x)
C(x)=A(x)B(x)这个原本需要三次DFT的操作变成只需要两次。
设
F
(
x
)
=
A
(
x
)
+
i
B
(
x
)
,
G
(
x
)
=
A
(
x
)
−
i
B
(
x
)
F(x)=A(x)+iB(x),G(x)=A(x)-iB(x)
F(x)=A(x)+iB(x),G(x)=A(x)−iB(x),
F
(
x
)
和
G
(
x
)
F(x)和G(x)
F(x)和G(x)是共轭的,那么
D
F
T
F
(
x
)
和
D
F
T
G
(
x
)
DFT_F(x)和DFT_G(x)
DFTF(x)和DFTG(x)也是共轭的。
即
D
F
T
G
(
x
)
=
c
o
n
j
(
D
F
T
F
(
n
−
x
)
)
DFT_G(x)=conj(DFT_F(n-x))
DFTG(x)=conj(DFTF(n−x)),这里记
c
o
n
j
(
x
)
conj(x)
conj(x)表示
x
x
x的共轭复数(
(
a
+
b
i
)
(a+bi)
(a+bi)的共轭复数是
(
a
−
b
i
)
(a-bi)
(a−bi))。
有了这个我们只需要算出
D
F
T
F
(
x
)
DFT_F(x)
DFTF(x)就能求出
D
F
T
G
(
x
)
DFT_G(x)
DFTG(x)了,那么就可以解出
D
F
T
A
(
x
)
和
D
F
T
B
(
x
)
DFT_A(x)和DFT_B(x)
DFTA(x)和DFTB(x)了,省去了一次DFT。
应用到拆系数FFT
A
(
x
)
=
a
W
+
b
,
B
(
x
)
=
c
W
+
d
A(x)=aW+b,B(x)=cW+d
A(x)=aW+b,B(x)=cW+d
A
(
x
)
B
(
x
)
=
a
c
W
2
+
(
a
d
+
b
c
)
W
+
b
d
A(x)B(x)=acW^2+(ad+bc)W+bd
A(x)B(x)=acW2+(ad+bc)W+bd
我们可以这样求
a
c
,
a
d
,
b
c
,
b
d
ac,ad,bc,bd
ac,ad,bc,bd:
P
(
x
)
=
(
a
+
i
b
)
(
c
+
i
d
)
=
(
a
c
−
b
d
)
+
(
a
d
+
b
c
)
i
P(x)=(a+ib)(c+id)=(ac-bd)+(ad+bc)i
P(x)=(a+ib)(c+id)=(ac−bd)+(ad+bc)i
Q
(
x
)
=
(
a
−
i
b
)
(
c
+
i
d
)
=
(
a
c
+
b
d
)
+
(
a
d
−
b
c
)
i
Q(x)=(a-ib)(c+id)=(ac+bd)+(ad-bc)i
Q(x)=(a−ib)(c+id)=(ac+bd)+(ad−bc)i
(
a
+
i
b
)
和
(
a
−
i
b
)
(a+ib)和(a-ib)
(a+ib)和(a−ib)的点值可以一次求出,
(
c
+
i
d
)
(c+id)
(c+id)也可以一次求出,再IDFT求出
P
(
x
)
和
Q
(
x
)
P(x)和Q(x)
P(x)和Q(x),最后就可以解出每一个位置
a
c
,
b
d
,
a
d
,
b
c
ac,bd,ad,bc
ac,bd,ad,bc的值了。
PS:为了保证精度要用long double,我的方法常数还是比较大的。
#include<cstdio>#include<cmath>#include<cstring>#include<algorithm>#define maxn 500005#define db long double #define ll long longusingnamespace std;const db Pi=acos(-1.0);int n,m,i,j,k,a[maxn],b[maxn],lim,bt[maxn],mo;
ll H[maxn],W[maxn];intI(db x){return(x<0)?-1:1;}struct Z{
db x,y;Z(db _x=0,db _y=0){x=_x,y=_y;}} A[maxn],B[maxn],C[maxn],F[maxn],G[maxn],w[maxn];
Z operator+(Z a,Z b){returnZ(a.x+b.x,a.y+b.y);}
Z operator*(Z a,Z b){returnZ(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
Z operator-(Z a,Z b){returnZ(a.x-b.x,a.y-b.y);}voiddft(Z *a,int sig){for(int i=0;i<lim;i++)if(i<bt[i])swap(a[i],a[bt[i]]);for(int mid=1;mid<lim;mid<<=1){
Z Wi(cos(Pi/mid),sig*sin(Pi/mid));for(int j=0;j<lim;j+=mid<<1){
Z w(1,0);for(int k=0;k<mid;k++,w=w*Wi){
Z x=a[j+k],y=a[j+k+mid]*w;
a[j+k]=x+y,a[j+k+mid]=x-y;}}}}intmain(){scanf("%d%d%d",&n,&m,&mo);for(i=0;i<=n;i++)scanf("%d",&a[i]),a[i]%=mo;for(i=0;i<=m;i++)scanf("%d",&b[i]),b[i]%=mo;for(lim=1;lim<=n+m;lim<<=1);for(i=1;i<lim;i++) bt[i]=(bt[i>>1]>>1)|((i&1)*(lim>>1));for(i=0;i<=n;i++) A[i].x=a[i]>>15,A[i].y=a[i]&((1<<15)-1);for(i=0;i<=m;i++) B[i].x=b[i]>>15,B[i].y=b[i]&((1<<15)-1);dft(A,1),dft(B,1);for(i=1;i<lim;i++) C[i]=A[lim-i],C[i].y=-C[i].y;
C[0]=A[0],C[0].y=-C[0].y;for(i=0;i<lim;i++) F[i]=A[i]*B[i],G[i]=C[i]*B[i];dft(F,-1),dft(G,-1);
ll K=(1<<15)%mo,K2=(1<<30)%mo;for(i=0;i<lim;i++){
ll x=(abs(F[i].x/lim)+0.5)*I(F[i].x),xx=(abs(F[i].y/lim)+0.5)*I(F[i].y);
ll y=(abs(G[i].x/lim)+0.5)*I(G[i].x),yy=(abs(G[i].y/lim)+0.5)*I(G[i].y);
H[i]=(((x+y)/2%mo*K2+xx%mo*K+(y-x)/2)%mo+mo)%mo;}for(i=0;i<=n+m;i++)printf("%lld ",H[i]);}