飛飛兔家族

多項式

如何學習這些毒瘤？

背模板！......反正學了不久就忘光了，仍是背模板來的實在。

1.FFT

快速傅里葉變換。一切的基礎。

做用是在 點值式 與 係數式 之間進行轉換。

首先咱們要手寫複數。

而後背下來這個預處理的東西：

r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));

把n補全爲2的次冪，下界爲 n + m + 1，數組開4倍以上。

而後注意把n補全以後的各個地方n的寫法：

FFT中都是 < n，涉及要實際輸出的是 <= n 

最後別忘了 / n，輸出時 + 0.5，輸出 [0, n + m]

中間也要記清楚，Wn = (cos(pi / len), f * sin(pi / len));

手寫pi的精度高一些，3.1415926 5358979 32384626

FFT前記得置換。

上代碼：

 1 #include <cstdio>
 2 #include <algorithm>
 3 #include <cmath>
 4 const int N = 1000010;
 5 const double pi = 3.1415926535897932384626;
 6 
 7 struct cp {
 8     double x, y;
 9     cp(double tx = 0.0, double ty = 0.0) {
10         this->x = tx;
11         this->y = ty;
12     }
13     inline cp operator +(const cp &d) const {
14         return cp(x + d.x, y + d.y);
15     }
16     inline cp operator -(const cp &d) const {
17         return cp(x - d.x, y - d.y);
18     }
19     inline cp operator *(const cp &d) const { // 複數的乘法直接代數展開就能獲得公式 
20         return cp(x * d.x - y * d.y, x * d.y + y * d.x);
21     }
22 }A[N << 2], B[N << 2]; // 空間開單個多項式的4倍
23 
24 int r[N << 2], n, lm;
25 
26 inline void init() {
27     while((1 << lm) < n) {
28         lm++;
29     }
30     n = 1 << lm;
31     for(int i = 1; i < n; i++) { // 這裏從 0 或 1 開始,到 n 或 n - 1 都行 
32         r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1)); // 背誦! 
33     }
34     return;
35 }
36 
37 inline void FFT(int n, cp *a, int f) {
38     for(int i = 0; i < n; i++) {
39         if(i < r[i]) {
40             std::swap(a[i], a[r[i]]); // 先轉換 
41         }
42     }
43     for(int len = 1; len < n; len = len << 1) { // len 是倍增的 
44         cp Wn(cos(pi / len), f * sin(pi / len)); // 背誦!
45         for(int i = 0; i < n; i += (len << 1)) { // 加的是 (len << 1) !
46             cp w(1, 0);                          // 這裏三個 < 沒有 <= 
47             for(int j = 0; j < len; j++) {
48                 cp t = a[i + len + j] * w;
49                 a[i + len + j] = a[i + j] - t; /// error : t -> w
50                 a[i + j] = a[i + j] + t;
51                 w = w * Wn;
52             }
53         }
54     }
55     if(f == -1) {
56         for(int i = 0; i <= n; i++) {
57             a[i].x /= n; // 除n 
58         }
59     }
60     return;
61 }
62 
63 int main() {
64     int n1, n2;
65     scanf("%d%d", &n1, &n2);
66     for(int i = 0; i <= n1; i++) {
67         scanf("%lf", &A[i].x);
68     }
69     for(int i = 0; i <= n2; i++) {
70         scanf("%lf", &B[i].x);
71     }
72 
73     n = n1 + n2 + 1;
74     init();
75 
76     FFT(n, A, 1);
77     FFT(n, B, 1);
78     for(int i = 0; i <= n; i++) {
79         A[i] = A[i] * B[i];
80     }
81     FFT(n, A, -1);
82 
83     for(int i = 0; i <= n1 + n2; i++) {
84         printf("%d ", (int)(A[i].x + 0.5));
85     }
86 
87     return 0;
88 }

洛谷 P3803 AC代碼

FFT傳參1是把係數式轉爲點式，-1是逆運算。

用以解決多項式相乘，卷積等。

 2.NTT

快速數論變換，number theory transformation

大概就是把w改爲了g，在%%%意義下進行FFT。

與FFT的主要區別在於Wn那裏，以及各類取模逆元。

模板在下面一塊兒上。

3.多項式求逆元

反正就是一個極其毒瘤的東西...多項式不只有逆元還有ln，exp，開根...

如何求逆元?

首先在n = 1的時候逆元是常數項的逆元。

以後用倍增，列方程：

B(x) - B'(x) = 0 

平方同乘A(x)，解方程：

B(x) = 2 * B'(x) - A(x) * B'²(x)

轉換爲點值相乘：

B[i] = B'[i] * (2 - A[i] * B'[i]);

而後就能夠遞歸求解了。

下面這個是 NTT + 求逆 的模板。

  1 #include <cstdio>
  2 #include <algorithm>
  3 #include <cstring>
  4 typedef long long LL;
  5 const int N = 1000010;
  6 const LL MO = 998244353, g = 3;
  7 
  8 int r[N << 2];
  9 LL A[N << 2], B[N << 2], C[N << 2];
 10 
 11 inline LL qpow(LL a, LL b) {
 12     a %= MO;
 13     LL ans = 1;
 14     while(b) {
 15         if(b & 1) {
 16             ans = ans * a % MO;
 17         }
 18         a = a * a % MO;
 19         b = b >> 1;
 20     }
 21     return ans;
 22 }
 23 
 24 inline void NTT(int np, LL *a, int f) {
 25     int n = 1;
 26     while(n < np) {
 27         n = n << 1;
 28     }
 29     for(int i = 0; i < n; i++) {
 30         if(i < r[i]) {
 31             std::swap(a[i], a[r[i]]);
 32         }
 33     }
 34 
 35     for(int len = 1; len < n; len = len << 1) {
 36         LL Wn = qpow(g, (MO - 1) / (len << 1)); // 注意這兩行! 
 37         if(f == -1) {
 38             Wn = qpow(Wn, MO - 2); // 注意這兩行! 
 39         }
 40         for(int i = 0; i < n; i += (len << 1)) {
 41             LL w = 1;
 42             for(int j = 0; j < len; j++) {
 43                 LL t = a[i + len + j] * w % MO;
 44                 a[i + len + j] = (a[i + j] - t + MO) % MO;
 45                 a[i + j] = (a[i + j] + t) % MO;
 46                 w = w * Wn % MO;
 47             }
 48         }
 49     }
 50 
 51     if(f == -1) {
 52         LL inv = qpow(n, MO - 2);
 53         for(int i = 0; i <= np; i++) {
 54             a[i] = a[i] * inv % MO;
 55         }
 56     }
 57     return;
 58 }
 59 
 60 void getinv(LL *A, LL *B, int len) {
 61     if(len == 1) {
 62         B[0] = qpow(A[0], MO - 2);
 63         return;
 64     }
 65     getinv(A, B, (len + 1) >> 1); // len 要加 1 
 66     int n = 1, lm = 0;
 67     while(n < (len << 1)) {
 68         n = n << 1;
 69         lm++;
 70     }
 71     for(int i = 1; i < n; i++) {
 72         r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));
 73     }
 74 
 75     memcpy(C, A, len * sizeof(LL));
 76     memset(C + len, 0, (n - len) * sizeof(LL)); // 兩種寫法 
 77     /*for(int i = 0; i < len; i++) {
 78         C[i] = A[i];
 79     }
 80     for(int i = len; i < n; i++) {
 81         C[i] = 0;
 82     }*/
 83 
 84     NTT(n, C, 1);
 85     NTT(n, B, 1); // 這裏對B的兩次NTT最好不要去掉，目測會WA 
 86     for(int i = 0; i < n; i++) {
 87         B[i] = (2 - C[i] * B[i] % MO) * B[i] % MO;
 88         if(B[i] < 0) {
 89             B[i] += MO;
 90         }
 91     }
 92     NTT(n, B, -1);
 93     for(int i = len; i < n; i++) {
 94         B[i] = 0; // 別忘了賦0 
 95     }
 96     return;
 97 }
 98 
 99 int main() {
100     int n;
101     scanf("%d", &n);
102     for(int i = 0; i < n; i++) {
103         scanf("%lld", &A[i]);
104     }
105     getinv(A, B, n);
106     for(int i = 0; i < n; i++) {
107         printf("%lld ", B[i]);
108     }
109 
110     return 0;
111 }

洛谷 P4238 AC代碼

4.多項式除法 & 取模

利用多項式求逆實現。

說多了都是淚啊...先貼個初級模板上來。

  1 #include <cstdio>
  2 #include <algorithm>
  3 #include <cstring>
  4 typedef long long LL;
  5 const int N = 100010;
  6 const LL MO = 998244353, g = 3;
  7 
  8 int r[N << 2];
  9 LL a[N << 2], b[N << 2], C[N << 2], invB[N << 2], f[N << 2], q[N << 2];
 10 
 11 inline LL qpow(LL a, LL b) {
 12     LL ans = 1;
 13     a %= MO;
 14     while(b) {
 15         if(b & 1) {
 16             ans = ans * a % MO;
 17         }
 18         a = a * a % MO;
 19         b = b >> 1;
 20     }
 21     return ans;
 22 }
 23 
 24 inline void NTT(int np, LL *a, int f) {
 25     int n = 1;
 26     while(n < np) {
 27         n = n << 1;
 28     }
 29     for(int i = 1; i < np; i++) {
 30         if(i < r[i]) {
 31             std::swap(a[i], a[r[i]]);
 32         }
 33     }
 34 
 35     for(int len = 1; len < n; len = len << 1) {
 36         LL Wn = qpow(g, (MO - 1) / (len << 1));
 37         if(f == -1) {
 38             Wn = qpow(Wn, MO - 2);
 39         }
 40         for(int i = 0; i < n; i += (len << 1)) {
 41             LL w = 1;
 42             for(int j = 0; j < len; j++) {
 43                 LL t = w * a[i + len + j] % MO;
 44                 a[i + len + j] = (a[i + j] - t + MO) % MO;
 45                 a[i + j] = (a[i + j] + t) % MO;
 46                 w = w * Wn % MO;
 47             }
 48         }
 49     }
 50 
 51     if(f == -1) {
 52         LL inv = qpow(n, MO - 2);
 53         for(int i = 0; i <= n; i++) {
 54             a[i] = a[i] * inv % MO;
 55         }
 56     }
 57     return;
 58 }
 59 
 60 void getinv(LL *A, LL *B, int len) {
 61     if(len == 1) {
 62         B[0] = qpow(A[0], MO - 2);
 63         return;
 64     }
 65     getinv(A, B, (len + 1) >> 1);
 66 
 67     int n = 1, lm = 0;
 68     while(n < (len << 1)) {
 69         n = n << 1;
 70         lm++;
 71     }
 72     for(int i = 1; i < n; i++) {
 73         r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));
 74     }
 75 
 76     memcpy(C, A, len * sizeof(LL));
 77     memset(C + len, 0, (n - len) * sizeof(LL));
 78     NTT(n, C, 1);
 79     NTT(n, B, 1);
 80     for(int i = 0; i < n; i++) {
 81         B[i] = (2 - C[i] * B[i] % MO) * B[i] % MO;
 82         if(B[i] < 0) {
 83             B[i] += MO;
 84         }
 85     }
 86     NTT(n, B, -1);
 87     for(int i = len; i < n; i++) {
 88         B[i] = 0;
 89     }
 90 
 91     return;
 92 }
 93 
 94 inline void div(LL *A, LL *B, LL *F, LL *R, int np, int m) {
 95     std::reverse(A, A + np + 1);
 96     std::reverse(B, B + m + 1);
 97     getinv(B, invB, np - m + 1);
 98 
 99     int len = np - m;
100     int n = 1, lm = 0;
101     while(n < (len << 1)) {
102         n = n << 1;
103         lm++;
104     }
105     for(int i = 1; i < n; i++) {
106         r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));
107     }
108 
109     for(int i = 0; i <= len; i++) {
110         C[i] = A[i];
111     }
112     for(int i = len + 1; i < n; i++) {
113         C[i] = 0;
114         invB[i] = 0;
115     }
116     NTT(n, C, 1);
117     NTT(n, invB, 1);
118     for(int i = 0; i < n; i++) {
119         F[i] = C[i] * invB[i] % MO;
120     }
121     NTT(n, F, -1);
122     std::reverse(F, F + len + 1);
123     for(int i = len + 1; i < n; i++) {
124         F[i] = 0;
125     }
126 
127     // get F
128 
129     std::reverse(A, A + np + 1);
130     std::reverse(B, B + m + 1);
131     n = 1, lm = 0;
132     while(n <= np) {
133         n = n << 1;
134         lm++;
135     }
136     for(int i = 1; i < n; i++) {
137         r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lm - 1));
138     }
139     for(int i = 0; i <= np; i++) {
140         C[i] = F[i];
141     }
142     for(int i = np + 1; i < n; i++) {
143         C[i] = 0;
144     }
145     for(int i = m + 1; i < n; i++) {
146         B[i] = 0;
147     }
148     NTT(n, B, 1);
149     NTT(n, C, 1);
150     for(int i = 0; i < n; i++) {
151         C[i] = C[i] * B[i] % MO;
152     }
153     NTT(n, C, -1);
154     for(int i = 0; i < m; i++) {
155         R[i] = (A[i] - C[i] + MO) % MO;
156     }
157 
158     return;
159 }
160 
161 int main() {
162     int n, m;
163     scanf("%d%d", &n, &m);
164     for(int i = 0; i <= n; i++) {
165         scanf("%lld", &a[i]);
166     }
167     for(int i = 0; i <= m; i++) {
168         scanf("%lld", &b[i]);
169     }
170     div(a, b, f, q, n, m);
171     for(int i = 0; i <= n - m; i++) {
172         printf("%lld ", f[i]);
173     }
174     puts("");
175     for(int i = 0; i < m; i++) {
176         printf("%lld ", q[i]);
177     }
178 
179 
180     return 0;
181 }

洛谷 P4512 AC代碼