多项式1.5·拆系数FFT

上一篇多项式:多项式1·普通的FFT


$\rm{0x01\quad Preface}$

今天是$Feb.19^{th}$,开学前一天,现在是北京时间下午$10:05$,我还剩除了英语物理的所有作业没做 … 耶,真开心。

由于特殊原因嘛,我暂时不会写的特别详细——毕竟还有作业没有做完,所以就先整理地仓促一点。

首先,为什么要拆系数?这是一个问题——直接$FFT$之后判个精度再取模不就得了?很喜闻乐见的是,这个题中的$value_{max}$到达了$1e9\cdot1e9 \cdot 1e5 = 1e23$的级别,不可以直接long long爆艹。

$\rm{0x02~Normal-Coefficient-Spliting~FFT}$

首先是最朴素的拆系数法。其原理简单得很,就是选择一个$M$,把$A(x)$和$B(x)$的各项系数表示成
$$
A_i = a_iM+b_i~(0\leq b_i <a_i) \\
B_i = c_iM+d_i~(0\leq d_i <c_i)
$$
的形式。然后我们做乘法的话,会有
$$
A_i\cdot B_i = a_ic_iM^2+(a_id_i+c_ib_i)M+b_id_i
$$
那么我们考虑,对于第一部分$a_ic_iM^2$我们可以通过一次$DFT$、一次$IDFT$算出来,;对于第二部分$(a_id_i+c_ib_i)M$ 我们可以通过两次$DFT$分别算出$a_id_i$、$c_ib_i$然后合并,之后一次$IDFT$求出来;对于最后一部分则只需要一次$DFT$、一次$IDFT$求出。所以算法流程一共$7$次$FFT$。

那我们考虑估计一下范围,此处不甚严谨地使用$\Theta$作为同阶的符号$^{[1]}$

$$b_i\cdot d_i \approx \Theta(M^2)$$

$$a_i\cdot c_i =\lfloor \frac{P}{M}\rfloor \cdot \lfloor \frac{P}{M}\rfloor = \Theta(\lfloor \frac{P^2}{M^2}\rfloor)$$

$$(a_id_i + c_ib_i)\cdot M =\Theta( M \cdot \lfloor \frac{P}{M}\rfloor) $$

那么我们取$M = \Theta(\sqrt P)$级别的,可以保证三个值的阶为$\Theta(P)$,大概是$1e9 \cdot 1e5 = 1e14$级别的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <cmath>
#include <cstdio>
#include <iostream>

#define MAXN 423333
#define rr register
#define rep(i, a, b) for(i = a ; i <= b ; ++ i)

using namespace std ;
const long double Pi = acos(-1) ;
int R[MAXN], L, x ; long long Ans[MAXN] ;
struct node{
long double x, y ;
node (long double xx = 0, long double yy = 0){
x = xx, y = yy ;
}
}A[MAXN], B[MAXN], C[MAXN], D[MAXN], F[MAXN], G[MAXN] ;
node H1[MAXN], H2[MAXN], H3[MAXN] ; int L1, L2, P, N, M ;
node operator + (node J, node Q){ return node(J.x + Q.x , J.y + Q.y); }
node operator - (node J, node Q){ return node(J.x - Q.x , J.y - Q.y); }
node operator * (node J, node Q){ return node(J.x * Q.x - J.y * Q.y , J.x * Q.y + J.y * Q.x) ; }

void FFT(node *J, int flag){
rr int i, j, k, l ;
for(i = 0; i < N; i ++)
if(i < R[i]) swap(J[i], J[R[i]]) ;
for(j = 1; j < N; j <<= 1){
node T(cos(Pi / j), flag * sin(Pi / j)) ;
for(k = 0; k < N; k += (j << 1) ){
node t(1, 0) ;
for(l = 0 ; l < j; l ++, t = t * T){
node Nx = J[k + l], Ny = t * J[k + j + l] ;
J[k + l] = Nx + Ny ;
J[k + j + l] = Nx - Ny ;
}
}
}
if (flag < 0) for (i = 0 ; i <= N ; ++ i) J[i].x = J[i].x / N + 0.5 ;
}
inline int qr(){
rr int k = 0 ; char c = getchar() ;
while (!isdigit(c)) c = getchar() ;
while (isdigit(c)) k = (k << 1) + (k << 3) + c - 48, c = getchar() ;
return k ;
}
signed main(){
rr int i, t ; cin >> L1 >> L2 >> P, M = 32767 ; N = 1, t = L1 + L2 ;
for (i = 0 ; i <= L1 ; ++ i) x = qr(), A[i].x = x / M, B[i].x = x % M ;
for (i = 0 ; i <= L2 ; ++ i) x = qr(), C[i].x = x / M, D[i].x = x % M ;
while(N <= t) N <<= 1, ++ L ; rep(i, 0, N) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;
FFT(A, 1), FFT(B, 1), FFT(C, 1), FFT(D, 1) ;
rep(i, 0, N) H1[i] = A[i] * C[i], H2[i] = B[i] * D[i], H3[i] = A[i] * D[i]+ C[i] * B[i] ;
FFT(H1, -1), FFT(H2, -1), FFT(H3, -1) ;
for (i = 0 ; i <= N ; ++ i)
Ans[i] = ((long long)H1[i].x * M % P * M % P + (long long)H2[i].x % P + (long long)H3[i].x * M % P) % P ;
for (i = 0 ; i <= L1 + L2 ; ++ i) printf("%lld ", Ans[i]) ; return 0 ;
}

所以最终的复杂度是$\boldsymbol{O(7\cdot P(n) \cdot n \log n)}$,其中$P(n)$是$FFT$自带的、大到不可忽略的常量因子。

但是……好像这个常数有点大诶,算上常数的话已经是$O(n \log^2n)$的级别了,于是——

$\rm{0x03~Conjugate~Optimization}$

源自毛啸的《再探快速傅立叶变换》。

我们思考这样两个多项式$\boldsymbol{P, Q} \in \mathbb{C}$
$$
\rm{P(x) = A(x) + iB(x)} \\
\rm{Q(x) = A(x) - iB(x)}
$$
我们不妨令$P’[k]$和$Q’[k]$为其$\text{DFT}$之后的序列,即$P’[k] =P(\omega_n^k),Q’[k] =Q(\omega_n^k) $。

同时,令$\text{conj(x)}$表示对$x$取共轭。

那么会有$^{[5]}$:
$$
\begin{align}
P’[k] &= A(\omega_{n}^{k}) + i B(\omega_{n}^{k}) \\
& = \sum_{j=0}^{n-1} A_{j} \omega_{n}^{jk} + i B_{j} \omega_{n}^{jk} \\
& = \sum_{j=0}^{n-1} (A_{j} + i B_{j}) \left(\cos \left(\frac{2 \pi jk}{n}\right) + i \sin \left(\frac{2 \pi jk}{n}\right)\right) \\
& = \sum_{j=0}^{n-1} (A_{j} + i B_{j}) \omega_{n}^{kj}
\\
Q’[k] &= A(\omega_{n}^{k}) - i B(\omega_{n}^{k}) \\
& = \sum_{j=0}^{n-1} A_{j} \omega_{n}^{jk} - i B_{j} \omega_{n}^{jk} \\
& = \sum_{j=0}^{n-1} (A_{j} - i B_{j}) \left(\cos \left(\frac{2 \pi jk}{n}\right) + i \sin \left(\frac{2 \pi jk}{n}\right)\right) \\
& = \sum_{j=0}^{n-1} \left(A_{j} \cos \left(\frac{2 \pi jk}{n}\right) + B_{j} \sin \left(\frac{2 \pi jk}{n}\right)\right) + i \left(A_{j} \sin \left(\frac{2 \pi jk}{n}\right) - B_{j} \cos \left(\frac{2 \pi jk}{n}\right)\right) \\
& = \text{conj} \left( \sum_{j=0}^{n-1} \left(A_{j} \cos \left(\frac{2 \pi jk}{n}\right) + B_{j} \sin \left(\frac{2 \pi jk}{n}\right)\right) - i \left(A_{j} \sin \left(\frac{2 \pi jk}{n}\right) - B_{j} \cos \left(\frac{2 \pi jk}{n}\right)\right) \right) \\
& = \text{conj} \left( \sum_{j=0}^{n-1} \left(A_{j} \cos \left(\frac{-2 \pi jk}{n}\right) - B_{j} \sin \left(\frac{-2 \pi jk}{n}\right)\right) + i \left(A_{j} \sin \left(\frac{-2 \pi jk}{n}\right) + B_{j} \cos \left(\frac{-2 \pi jk}{n}\right)\right) \right) \\
& = \text{conj} \left( \sum_{j=0}^{n-1} (A_{j} + i B_{j}) \left(\cos \left(\frac{-2 \pi jk}{n}\right) + i \sin \left(\frac{-2 \pi jk}{n}\right)\right)\right) \\
& = \text{conj} \left( \sum_{j=0}^{n-1} (A_{j} + i B_{j}) \omega_{n}^{-jk} \right) \\
& = \text{conj} \left( \sum_{j=0}^{n-1} (A_{j} + i B_{j}) \omega_{n}^{(n-k)j} \right) \\
& = \text{conj} \left( P’[n-k] \right)
\end{align}
$$
好吧我承认这段推导过程的代码甚是壮观,于是并不是我自己写的qwq。

那么我们发现其中$A(x)$和$B(x)$可以通过$P,Q$推出来:
$$
A’[k] = \frac{P’[k] + Q’[k]}{2} \\
B’[k] = \frac{P’[k] - Q’[k]}{2i}
$$
哦,对了,当$k=0$时,由于不存在这一项,所以我们需要特判一下。并且由于我们的$P(x)$和$Q(x)$的实部和虚部都可以利用,所以我们对于七次$DFT$可以优化到$4$次$DFT$.

以下是共轭优化$FFT$的初号机:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <cmath>
#include <cstdio>
#include <iostream>

#define MAXN 423333
#define rr register
#define rep(i, a, b) for(i = a ; i <= b ; ++ i)

using namespace std ;
int f[MAXN], g[MAXN] ;
int R[MAXN], Len, x, Ans[MAXN] ;
const long double Pi = acos(-1) ;
struct node{
long double r, i ;
inline node (long double xx = 0, long double yy = 0){
r = xx, i = yy ;
}
inline node Conj() {
return node(r, -i);
}
}A[MAXN], B[MAXN], w[MAXN], t1[MAXN], t2[MAXN] ;
node H1[MAXN], H2[MAXN], H3[MAXN] ; int L1, L2, P, N, M ;
node operator + (node J, node Q){ return node(J.r + Q.r , J.i + Q.i); }
node operator - (node J, node Q){ return node(J.r - Q.r , J.i - Q.i); }
node operator * (node J, node Q){ return node(J.r * Q.r - J.i * Q.i , J.r * Q.i + J.i * Q.r) ; }
inline int qr(){
rr int k = 0 ; char c = getchar() ; while (!isdigit(c)) c = getchar() ;
while (isdigit(c)) k = (k << 1) + (k << 3) + c - 48, c = getchar() ; return k ;
}
void FFT(node *J, int flag){
rr int i, j, k, l ;
for(i = 0; i < N; i ++)
if(i < R[i]) swap(J[i], J[R[i]]) ;
for(j = 1; j < N; j <<= 1)
for(k = 0; k < N; k += (j << 1))
for(l = 0 ; l < j; ++ l){
node T = w[N / j * l] ; T.i *= flag ;
node Nx = J[k + l], Ny = T * J[k + j + l] ;
J[k + l] = Nx + Ny, J[k + j + l] = Nx - Ny ;
}
}
inline void Init(int L){
rr int i ; while (N <= L) N <<= 1, ++ Len ;
for (i = 0 ; i < N ; ++ i) A[i] = node(f[i] & 32767, f[i] >> 15) ;
for (i = 0 ; i < N ; ++ i) B[i] = node(g[i] & 32767, g[i] >> 15) ;
for (i = 0 ; i < N ; ++ i) w[i] = node(cos(Pi * i / N), sin(Pi * i / N)) ;
for (i = 0 ; i < N ; ++ i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (Len - 1)) ;
}
void MTT(){
node ia, ib, a1, a2, b1, b2 ;
rr int i, k, t = L1 + L2, q1, q2, q3 ; Init(t) ;
// for (i = 0 ; i <= N ; ++ i) cout << A[i].r << " "; cout << endl ;
// for (i = 0 ; i <= N ; ++ i) cout << B[i].r << " "; cout << endl ;
FFT(A, 1), FFT(B, 1) ;
for (i = 0 ; i < N ; ++ i){
k = (N - i) & (N - 1), ia = A[k].Conj(), ib = B[k].Conj() ;//(ia,ib) = Q
a1 = (ia + A[i]) * node(0.5, 0), a2 = (A[i] - ia) * node(0, -0.5) ;
b1 = (ib + B[i]) * node(0.5, 0), b2 = (B[i] - ib) * node(0, -0.5) ;
t1[i] = a1 * b1 + (a1 * b2 + b1 * a2) * node(0, 1), t2[i] = a2 * b2 ;
} FFT(t1, -1), FFT(t2, -1) ;
// for (i = 0 ; i <= N ; ++ i) cout << w[i].r << " "; cout << endl ;
// for (i = 0 ; i <= N ; ++ i) cout << A[i].r << " "; cout << endl ;
// for (i = 0 ; i <= N ; ++ i) cout << B[i].r << " "; cout << endl ;
for (i = 0 ; i < N ; ++ i){
q1 = (long long)(t1[i].r / N + 0.5) % P, q2 = (long long)(t1[i].i / N + 0.5) % P ;
q3 = (long long)(t2[i].r / N + 0.5) % P, Ans[i] = ((((long long)q3 << 30) % P + ((long long)q2 << 15) % P + q1) % P + P) % P ;
}
}
signed main(){
rr int i ;
cin >> L1 >> L2 >> P ; N = 1 ;
for (i = 0 ; i <= L1 ; ++ i) f[i] = qr() % P ;
for (i = 0 ; i <= L2 ; ++ i) g[i] = qr() % P ; MTT() ;
for (i = 0 ; i <= L1 + L2 ; ++ i) printf("%d ", Ans[i]) ; return 0 ;
}

$\rm{0x04\quad}$拼命卡常

好的,首先我们可以欣赏一下最初的版本(用小号交的拆系数$FFT$ + $O2$)

然后是大号的共轭优化$FFT$(不加$O2$,即上方代码):

特别的,以下是无共轭优化的拆系数$FFT$,不开$O2$:

好吧,他看起来没有快多少。毕竟都是同阶的复杂度,好像后者的常数更大那么一点……

于是考虑对共轭优化的进行大力卡常:

  • 多次使用的非全局变量。使用register修饰符。
  • 将$double$转换成为$long~double$ 。
  • 减少取模次数。
  • 从yjk那里偷来的$\rm{fread/fwrite}$

一番操作之后,我们成功地卡到了第五页上……

最后奉上最快的代码qwq:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
// luogu-judger-enable-o2
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")

#include <bits/stdc++.h>
#define MAXN 102333
#define MAXM 272333

#define rr register
#define ll long long
#define rep(i, a, b) for(i = a ; i <= b ; ++ i)

using namespace std ;

namespace IO{
const int ch_top=4e7+3;
char ch[ch_top],*now_r=ch-1,*now_w=ch-1;
inline int read(){
while(*++now_r<'0');
register int x=*now_r-'0';
while(*++now_r>='0')x=x*10+*now_r-'0';
return x;
}
inline void write(int x){
static char st[20];static int top;
while(st[++top]='0'+x%10,x/=10);
while(*++now_w=st[top],--top);
*++now_w=' ';
}
}

int R[MAXM], Len, x, Ans[MAXN << 1] ;
const double Pi = acos(-1) ;

struct node{
double r, i ;
inline node Conj() { return node(r, -i); }
inline node (double xx = 0, double yy = 0){
r = xx, i = yy ;
}//定义的Complex型是用来进行复数运算的
}A[MAXM], B[MAXM], w[MAXM], t1[MAXM], t2[MAXM] ; int L1, L2, P, N, M ;
inline node operator + (const node &J, const node &Q) { return node(J.r + Q.r , J.i + Q.i); }
inline node operator - (const node &J, const node &Q) { return node(J.r - Q.r , J.i - Q.i); }
inline node operator * (const node &J, const double &Q) { return (node) {J.r * Q, J.i * Q} ; }
inline node operator * (const node &J, const node &Q) { return node(J.r * Q.r - J.i * Q.i , J.r * Q.i + J.i * Q.r) ; }

inline void FFT(node *J){
rr node t ;
rr int i, j, k, l ;
for(i = 0; i < N; i ++)
if(i < R[i]) t = J[i], J[i] = J[R[i]], J[R[i]] = t ;
for(j = 1; j < N; j <<= 1)
for(k = 0; k < N; k += (j << 1))
for(l = 0 ; l < j; ++ l){
rr node T = w[N / j * l] ;
rr node Nx = J[k + l], Ny = T * J[k + j + l] ;
J[k + l] = Nx + Ny, J[k + j + l] = Nx - Ny ;
}
}
inline void IFFT(node *J){
reverse(J + 1, J + N), FFT(J) ;
rr int i ; rr double qwq = 1.0 / N ;
rep(i, 0, N - 1) J[i] = J[i] * qwq ;
}

using namespace IO ;
int main(){
fread(ch,1,ch_top,stdin);
rr int i, k, t ; rr node ia, ib, a1, a2, b1, b2 ;
t = ((L1 = read()) + (L2 = read())), P = read(), N = 1 ;
while (N <= t) N <<= 1, ++ Len ;
for (i = 0 ; i <= L1 ; ++ i) x = read(), A[i] = node(x & 32767, x >> 15) ;
for (i = 0 ; i <= L2 ; ++ i) x = read(), B[i] = node(x & 32767, x >> 15) ;
for (i = 0 ; i < N ; ++ i) w[i] = node(cos(Pi * i / N), sin(Pi * i / N)), R[i] = (R[i >> 1] >> 1) | ((i & 1) << (Len - 1)) ;
FFT(A), FFT(B) ;
for (i = 0 ; i < N ; ++ i){
k = (N - i) & (N - 1), ia = A[k].Conj(), ib = B[k].Conj() ;//(ia,ib) = Q
a1 = (ia + A[i]) * node(0.5, 0), a2 = (A[i] - ia) * node(0, -0.5) ;
b1 = (ib + B[i]) * node(0.5, 0), b2 = (B[i] - ib) * node(0, -0.5) ;
t1[i] = a1 * b1 + (a1 * b2 + b1 * a2) * node(0, 1), t2[i] = a2 * b2 ;
}
IFFT(t1), IFFT(t2) ;
for (i = 0 ; i <= t ; ++ i)
write(((ll)(t1[i].r + 0.5) + ((ll)(t1[i].i + 0.5) % P << 15) + ((ll)(t2[i].r + 0.5) % P << 30)) % P) ;
fwrite(ch,1,now_w-ch,stdout); return 0 ;
}

实践证明,以上代码不加任何优化(不开$\rm{O2/3/fast}$)甚至可以快$4ms$!

$\rm{0x00\quad Afterword}$

嗯,其实按道理来讲,不是特别难。但是这跟$HLPP$一样,都是打死都不会考的算法,所以学这些只是为了娱乐……听起来挺苍凉?

但似乎,从一开始就不应该把应付考点作为OI的初衷吧,虽然如果没有获利,没有多少人会去学,但是不沾染功利的OI,似乎可爱了那么几分呢……

并且在学的过程中顺便认识了一个巨佬CMXRYNP,嘿嘿,也算不亏啦。

本篇文章真实完稿时间是$2019/3/17$,因为太懒+太忙,鸽了一个月$\rm{qaq}$.

$\rm{Reference}$

-------------本文结束感谢您的阅读-------------