FFT与多项式、生成函数题目泛做

题目1 COGS 很强的乘法问题

高精度乘法用FFT加速

 1 #include <cstdlib>
 2 #include <iostream>
 3 #include <algorithm>
 4 #include <cstring>
 5 #include <cstdio>
 6 #include <cmath>
 7 #define Pi acos(-1)
 8 
 9 using namespace std;
10 const int N = 600000 + 5;
11 const int rad = 10;
12 
13 struct Complex {
14   double r, i;
15   Complex(double x = 0, double y = 0) { r = x; i = y; }
16   Complex operator + (Complex x) { return Complex(r + x.r, i + x.i); }
17   Complex operator - (Complex x) { return Complex(r - x.r, i - x.i); }
18   Complex operator * (Complex x) { return Complex(r * x. r - i * x.i, r * x.i + i * x.r); }
19 }a[N], b[N];;
20 
21 int n, m, l, r[N], c[N];
22 char ss[N];
23 
24 void dft(Complex *s, int f) {
25   for(int i = 0; i < n; ++ i)
26     if(i < r[i]) swap(s[i], s[r[i]]);
27   for(int i = 1; i < n; i <<= 1) {
28     Complex wn(cos(Pi / i), f * sin(Pi / i));
29     for(int j = 0; j < n; j += (i << 1)) {
30       Complex w(1, 0);
31       for(int k = 0; k < i; ++ k) {
32         Complex x = s[j + k], y = w * s[j + k + i];
33         s[j + k] = x + y;
34         s[j + k + i] = x - y;
35         w = w * wn;
36       }
37     }
38   }
39   if(f == -1)
40     for(int i = 0; i < n; ++ i)
41       s[i].r /= n;
42 }
43 
44 int main() {
45 #ifndef stone
46   freopen("bettermul.in", "r", stdin);
47   freopen("bettermul.out", "w", stdout);
48 #endif
49 
50   int len1, len2;
51   scanf("%s", ss); len1 = strlen(ss);
52   for(int i = 0; i < len1; ++ i) a[i] = ss[len1 - i - 1] - '0';
53   scanf("%s", ss); len2 = strlen(ss);
54   for(int i = 0; i < len2; ++ i) b[i] = ss[len2 - i - 1] - '0';
55   n = max(len1, len2); m = (n << 1) - 2;
56   for(n = 1; n <= m; n <<= 1) l ++;
57   for(int i = 0; i < n; ++ i)
58     r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
59   dft(a, 1); dft(b, 1);
60   for(int i = 0; i < n; ++ i) a[i] = a[i] * b[i];
61   dft(a, -1);
62   for(int i = 0; i <= m; ++ i) c[i] = (int) (a[i].r + 0.5);
63   for(int i = 0; i <= m; ++ i) {
64     if(c[i] >= rad) {
65       c[i + 1] += c[i] / rad;
66       c[i] %= rad;
67       if(i == m) m ++;
68     }
69   }
70   int pos = m;
71   while(c[pos] == 0) -- pos;
72   for(int i = pos; i >= 0; -- i) printf("%d", c[i]);
73 
74 #ifndef stone
75   fclose(stdin); fclose(stdout);
76 #endif
77   return 0;
78 }

COGS

题目2 BZOJ4259 残缺的字符串

算法讨论：（摘自Claris的博客，写得很清楚，我就是看这个才会做的）

注意：如果开第二个字符串的倍长的话，FFT会超时。还要注意统计答案的区间是什么。不要错了。

  1 #include <iostream>
  2 #include <algorithm>
  3 #include <cstdlib>
  4 #include <cstring>
  5 #include <cstdio>
  6 #include <cmath>
  7  
  8 using namespace std;
  9 typedef long long ll;
 10 const int N = 600000 + 5;
 11 #define Pi acos(-1)
 12  
 13 struct Cp {
 14   double r, i;
 15   Cp(double _r = 0, double _i = 0):
 16     r(_r), i(_i) {}
 17   Cp operator + (Cp x) {
 18     return Cp(r + x.r, i + x.i);
 19   }
 20   Cp operator - (Cp x) {
 21     return Cp(r - x.r, i - x.i);
 22   }
 23   Cp operator * (Cp x) {
 24     return Cp(r * x.r - i * x.i, r * x.i + i * x.r);
 25   }
 26 };
 27  
 28 char str1[N], str2[N];
 29 int len1, len2, n, m, l, cnt;
 30 int t[N], ans[N], r[N];
 31 Cp a[N], b[N], A[N], B[N], C[N], tp;
 32  
 33  
 34 void DFT(Cp *s, int f) {
 35   for(int i = 0; i < n; ++ i)
 36     if(i < r[i]) swap(s[i], s[r[i]]);
 37   for(int i = 1; i < n; i <<= 1) {
 38     Cp wn(cos(Pi / i), f * sin(Pi / i));
 39     for(int j = 0; j < n; j += (i << 1)) {
 40       Cp w(1, 0);
 41       for(int k = 0; k < i; ++ k) {
 42         Cp x = s[j + k], y = w * s[j + k + i];
 43         s[j + k] = x + y;
 44         s[j + k + i] = x - y;
 45         w = w * wn;
 46       }
 47     }
 48   }
 49   if(f == -1)
 50     for(int i = 0; i < n; ++ i)
 51       s[i].r /= n;
 52 }
 53  
 54 #define stone
 55  
 56 int main() {
 57 #ifndef stone
 58    
 59   freopen("bitch.in", "r", stdin);
 60   freopen("bitch.out", "w", stdout);
 61  
 62 #endif
 63    
 64   scanf("%d%d", &len1, &len2);
 65   scanf("%s%s", str1, str2);
 66   reverse(str1, str1 + len1);
 67   for(int i = 0; i < len1; ++ i) {
 68     if(str1[i] == '*') a[i].r = 0;
 69     else a[i].r = (str1[i] - 'a' + 1);
 70   }
 71   for(int i = 0; i < len2; ++ i) {
 72     if(str2[i] == '*') b[i].r = 0;
 73     else b[i].r = (str2[i] - 'a' + 1);
 74   }
 75   m = len1 + len2;
 76   for(n = 1; n <= m; n <<= 1) ++ l;
 77   for(int i = 0; i < n; ++ i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
 78  
 79   for(int i = 0; i < n; ++ i) {
 80     A[i].r = a[i].r * a[i].r * a[i].r;
 81     B[i].r = b[i].r;
 82   }
 83   DFT(A, 1); DFT(B, 1);
 84   for(int i = 0; i < n; ++ i) { C[i] = C[i] + A[i] * B[i]; }
 85   memset(A, 0, sizeof A); memset(B, 0, sizeof B);
 86   for(int i = 0; i < n; ++ i) {
 87     A[i].r = a[i].r * a[i].r;
 88     B[i].r = b[i].r * b[i].r;
 89   }
 90   DFT(A, 1); DFT(B, 1); tp = (Cp) {2, 0};
 91   for(int i = 0; i < n; ++ i) { C[i] = C[i] - A[i] * B[i] * tp; }
 92   memset(A, 0, sizeof A); memset(B, 0, sizeof B);
 93   for(int i = 0; i < n; ++ i) {
 94     A[i].r = a[i].r;
 95     B[i].r = b[i].r * b[i].r * b[i].r;
 96   }
 97   DFT(A, 1); DFT(B, 1);
 98   for(int i = 0; i < n; ++ i) {  C[i] = C[i] + A[i] * B[i]; }
 99   DFT(C, -1);
100   for(int i = len1 - 1; i < len2; ++ i)
101     if(C[i].r < 0.5) {
102       ans[++ cnt] = i - len1 + 2;
103     }
104   printf("%d
", cnt);
105   for(int i = 1; i <= cnt; ++ i) {
106     printf("%d ", ans[i]);
107   }
108  
109 #ifndef stone
110  
111   fclose(stdin); fclose(stdout);
112  
113 #endif
114    
115   return 0;
116 }

4259

题目3 BZOJ 万径人踪灭

算法讨论：

用FFT算出所有回文串（相邻的和不相邻的），然后用Manacher计算出所有相邻的回文串，然后做差统计答案即可。

用FFT算出所有回文串的时候，是有len + len - 1个位置要计算的，分别是在字符的位置上和在两个字符之间的位置。

  1 #include <cstdlib>
  2 #include <algorithm>
  3 #include <cstdio>
  4 #include <cstring>
  5 #include <iostream>
  6 #include <cmath>
  7  
  8 using namespace std;
  9 const int N = 400000 + 5;
 10 typedef long long ll;
 11 const int mod = 1e9 + 7;
 12 #define Pi acos(-1)
 13  
 14 struct Cp {
 15   double r, i;
 16   Cp(double _r = 0, double _i = 0) :
 17     r(_r), i(_i){}
 18   Cp operator + (Cp x) {
 19     return Cp(r + x.r, i + x.i);
 20   }
 21   Cp operator - (Cp x) {
 22     return Cp(r - x.r, i - x.i);
 23   }
 24   Cp operator * (Cp x) {
 25     return Cp(r * x.r - i * x.i, r * x.i + i * x.r);
 26   }
 27 };
 28  
 29 char ss[N], Ma[N << 1];
 30 int n, m, l, r[N], lens;
 31 Cp a[N], b[N];
 32 ll Mp[N << 1], t1[N], t2[N], p[N];
 33  
 34 void DFT(Cp *s, int f) {
 35   for(int i = 0; i < n; ++ i)
 36     if(i < r[i]) swap(s[i], s[r[i]]);
 37   for(int i = 1; i < n; i <<= 1) {
 38     Cp wn(cos(Pi / i), f * sin(Pi / i));
 39     for(int j = 0; j < n; j += (i << 1)) {
 40       Cp w(1, 0);
 41       for(int k = 0; k < i; ++ k) {
 42         Cp x = s[j + k], y = w * s[j + k + i];
 43         s[j + k] = x + y;
 44         s[j + k + i] = x - y;
 45         w = w * wn;
 46       }
 47     }
 48   }
 49   if(f == -1)
 50     for(int i = 0; i < n; ++ i)
 51       s[i].r /= n;
 52 }
 53  
 54 ll Manacher(char *s, int len) {
 55   int l = 0; ll res = 0;
 56   Ma[l ++] = '$';
 57   Ma[l ++] = '#';
 58   for(int i = 0; i < len; ++ i) {
 59     Ma[l ++] = s[i];
 60     Ma[l ++] = '#';
 61   }
 62   Ma[l] = '!';
 63   int Mx = 0, id = 0;
 64   for(int i = 0; i < l; ++ i) {
 65     if(Mx > i) {
 66       Mp[i] = min(Mp[2 * id - i], (ll)Mx - i);
 67     }
 68     else {
 69       Mp[i] = 1;
 70     }
 71     while(Ma[i + Mp[i]] == Ma[i - Mp[i]]) Mp[i] ++;
 72     if(i + Mp[i] > Mx) {
 73       Mx = i + Mp[i];
 74       id = i;
 75     }
 76     res += (Mp[i] >> 1); res %= mod;
 77   }
 78   return res;
 79 }
 80  
 81 int main() {
 82   ll ans = 0, tmp;
 83   scanf("%s", ss);
 84   n = strlen(ss); lens = strlen(ss);
 85   for(int i = 0; i < n; ++ i) {
 86     if(ss[i] == 'a')
 87       a[i].r = 1, b[i].r = 0;
 88     else a[i].r = 0, b[i].r = 1;
 89   }
 90   m = (n << 1) - 2;
 91   for(n = 1; n <= m; n <<= 1) l ++;
 92   for(int i = 0; i < n; ++ i) r[i] = (r[i >> 1] >> 1) |((i & 1) << (l - 1));
 93   DFT(a, 1); DFT(b, 1);
 94   for(int i = 0; i < n; ++ i) {
 95     a[i] = a[i] * a[i];
 96     b[i] = b[i] * b[i];
 97   }
 98   DFT(a, -1); DFT(b, -1);
 99   for(int i = 0; i <= n; ++ i) {
100     t1[i] = (ll) (a[i].r + 0.5);
101     t2[i] = (ll) (b[i].r + 0.5);
102   }
103   for(int i = 0; i < n; ++ i) {
104     t1[i] += t2[i];
105   }
106   tmp = Manacher(ss, strlen(ss));
107   p[0] = 1;
108   for(int i = 1; i <= lens; ++ i)
109     p[i] = p[i - 1] * 2 % mod;
110   lens = lens + lens - 1;
111   for(int i = 0; i < lens; ++ i) {
112     if(i & 1) {
113       ans = (ans + p[t1[i] >> 1] - 1) % mod;
114     }
115     else {
116       ans = (ans + (p[(t1[i] - 1) >> 1] << 1) - 1) % mod;
117     }
118   }
119   ans = (ans - tmp) % mod;
120   while(ans < 0) {
121     ans += mod;
122     ans %= mod;
123   }
124   printf("%lld
", ans);
125   return 0;
126 }

万径人踪灭

题目4 BZOJ3028

首先很高兴的表示这里的汉堡是承德的。然后表示如果会生成函数的话这就是个输入输出题。

系数是种类，指数是数量性质。

如下图：

然后就没有然后了。知道答案就考虑输入和输出了。表示0ms过没有压力。这题如果你要是写个搜索。。。。。

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <cstdlib>
 5 #include <algorithm>
 6 #include <cctype>
 7  
 8 using namespace std;
 9 const int N = 500 + 5;
10 const int mod = 10007;
11 typedef long long ll;
12  
13 ll n;
14  
15 inline ll read() {
16   ll x = 0;
17   char c = getchar();
18   while(!isdigit(c)) c = getchar();
19   while(isdigit(c)) {
20     x = x * 10 + c - '0';
21     x %= mod;
22     c = getchar();
23   }
24   return x;
25 }
26  
27 int main() {
28   n = read();
29   n = n * (n + 1) * (n + 2);
30   n /= 6;
31   n %= mod;
32   printf("%lld
", n);
33   return 0;
34 }

3028