标签：ng2 ems sizeof clear 还需要 names bsp class 还需
给定一个由小写字母构成的字符串$s$，要求添加$n(n\le 10^9)$个小写字母，求构成回文串的数目。
$n$辣么大，显然要矩阵快速幂嘛。
考虑从两端开始构造以s ss为子串的回文串，该回文串长度为$N=n+s$，每次添加相同的字符，则需要$(N+1)/2$次，则用dp来计算并使用矩阵乘法来优化转移会得到一个$O(|s|^6\log N)$的算法，显然是不可接受的。
考虑这个dp做法，设$f[i][j][k]$表示从两端添加了$k$次字符，原来的$s$的子串$s_{ij}$任然需要被构造出来时的方案数，我们发现会存在两种转移，一种是添加了一个字符，$s_{ij}$转移到$s_{i+1,j-1}$，这是因为$s_i=s_j$，添加一个字符可以消去$s$中两个字符，记为n25转移；另一种是$s_i\not = s_j$只消去一个字符，记为n25转移。同时要注意到，每种转移都有到自身的转移（即在回文串两端添加字符后$s_{ij}$不变，n25转移有25种添加方式转移到自身，n24转移有24种添加方式转移到自身）
我们发现其实dp可以拆成两部分来计算，一部分是计算存在$k$个n24转移的转移路径方案数，这样一来，n25的转移数为$(|s|-k+1)/2$，那么还需要添加$(N+1)/2 - k - (|s|-k+1)/2$个到自身的转移，这么一来，n24转移和n25转移的顺序就没有意义了，不妨把n24转移放到前面，n25转移放到后面，然后建图，从start到end的转移方案数就是这个图的临接矩阵的$(N+1)/2$次方。这么一来，复杂度是$O(|s|^4 \log N)$，依旧不可接受。但注意到这$k$种建图是可以合并成一个有$O(|s|)$个起点，$O(|s|)$个终点的图的，所以只需要在这一个图上做矩阵快速幂就好，复杂度为$O(|s|^3\log N)$，问题解决。还有一个优化是，该矩阵是上三角矩阵，可以利用这点来加速矩阵乘法。以上讨论基于$N$是偶数，当$N$是奇数时，注意到只有$s_{ii}$出发的n25转移是不合法的，减去即可。
1 #include <bits/stdc++.h> 2 using namespace std; 3 namespace my_header { 4 #define pb push_back 5 #define mp make_pair 6 #define pir pair<int, int> 7 #define vec vector<int> 8 #define pc putchar 9 #define clr(t) memset(t, 0, sizeof t) 10 #define pse(t, v) memset(t, v, sizeof t) 11 #define bl puts("") 12 #define wn(x) wr(x), bl 13 #define ws(x) wr(x), pc(‘ ‘) 14 const int INF = 0x3f3f3f3f; 15 typedef long long LL; 16 typedef double DB; 17 inline char gchar() { 18 char ret = getchar(); 19 for(; (ret == ‘\n‘ || ret == ‘\r‘ || ret == ‘ ‘) && ret != EOF; ret = getchar()); 20 return ret; } 21 template<class T> inline void fr(T &ret, char c = ‘ ‘, int flg = 1) { 22 for(c = getchar(); (c < ‘0‘ || ‘9‘ < c) && c != ‘-‘; c = getchar()); 23 if (c == ‘-‘) { flg = -1; c = getchar(); } 24 for(ret = 0; ‘0‘ <= c && c <= ‘9‘; c = getchar()) 25 ret = ret * 10 + c - ‘0‘; 26 ret = ret * flg; } 27 inline int fr() { int t; fr(t); return t; } 28 template<class T> inline void fr(T&a, T&b) { fr(a), fr(b); } 29 template<class T> inline void fr(T&a, T&b, T&c) { fr(a), fr(b), fr(c); } 30 template<class T> inline char wr(T a, int b = 10, bool p = 1) { 31 return a < 0 ? pc(‘-‘), wr(-a, b, 0) : (a == 0 ? (p ? pc(‘0‘) : p) : 32 (wr(a/b, b, 0), pc(‘0‘ + a % b))); 33 } 34 template<class T> inline void wt(T a) { wn(a); } 35 template<class T> inline void wt(T a, T b) { ws(a), wn(b); } 36 template<class T> inline void wt(T a, T b, T c) { ws(a), ws(b), wn(c); } 37 template<class T> inline void wt(T a, T b, T c, T d) { ws(a), ws(b), ws(c), wn(d); } 38 template<class T> inline T gcd(T a, T b) { 39 return b == 0 ? a : gcd(b, a % b); } 40 41 }; 42 using namespace my_header; 43 const int maxS = 200 + 20; 44 const int maxM = maxS * 2; 45 const int mod = 10007; 46 char str[maxS]; 47 int n, N, s, ans, f[maxS][maxS][maxS]; 48 inline void add(int &a, int b) { 49 (a += b) %= mod; 50 } 51 inline void sub(int &a, int b) { 52 (a -= b) %= mod; 53 if (a < 0) a += mod; 54 } 55 int dp(int l, int r, int k) { 56 int &ret = f[l][r][k]; 57 if (ret != -1) 58 return ret; 59 ret = 0; 60 if (l > r) 61 return 0; 62 if (l == r) 63 return ret = k == 0; 64 if (str[l] == str[r]) { 65 if (l == r - 1) 66 return ret = k == 0; 67 add(ret, dp(l + 1, r - 1, k)); 68 } else { 69 add(ret, dp(l + 1, r, k - 1)); 70 add(ret, dp(l, r - 1, k - 1)); 71 } 72 return ret; 73 } 74 struct Matrix { 75 int d[maxM][maxM]; 76 int n; 77 int* operator [] (int a) { 78 return d[a]; 79 } 80 void clear() { 81 memset(d, 0, sizeof d); 82 } 83 void init() { 84 memset(d, 0, sizeof d); 85 for (int i = 1; i <= n; ++i) 86 d[i][i] = 1; 87 } 88 }; 89 Matrix operator * (Matrix a, Matrix b) { 90 Matrix c; 91 c.n = a.n; 92 c.clear(); 93 for (int i = 1; i <= a.n; ++i) 94 for (int j = i; j <= a.n; ++j) 95 for (int k = i; k <= j; ++k) 96 add(c[i][j], a[i][k] * b[k][j]); 97 return c; 98 } 99 template<class T> inline T fpm(T b, int i) { 100 T r; 101 r.n = b.n; 102 r.init(); 103 for (; i; i >>= 1, b = b * b) 104 if (i & 1) 105 r = r * b; 106 return r; 107 } 108 int main() { 109 #ifdef lol 110 freopen("E.in", "r", stdin); 111 freopen("E.out", "w", stdout); 112 #endif 113 cin >> (str + 1) >> n; 114 s = strlen(str + 1); 115 memset(f, -1, sizeof f); 116 for (int i = 0; i <= s; ++i) 117 f[1][s][i] = dp(1, s, i); 118 N = (s + n + 1) >> 1; 119 int n25 = (s + 1) / 2; 120 int num = s + n25; 121 Matrix x, res; 122 x.clear(); 123 x.n = num + n25; 124 x[1][1] = 24; 125 for (int i = 2; i <= s; ++i) { 126 x[i][i] = 24; 127 x[i - 1][i] = 1; 128 } 129 130 x[s][s + 1] = 1; 131 x[s + 1][s + 1] = 25; 132 x[s + 1][num + 1] = 1; 133 x[num + 1][num + 1] = 26; 134 for (int i = s + 2; i <= s + n25; ++i) { 135 x[i][i] = 25; 136 x[i - 1][i] = 1; 137 x[i][num + i - s] = 1; 138 x[num + i - s][num + i - s] = 26; 139 } 140 res = fpm(x, N); 141 for (int k = 0; k <= s - 1; ++k) { 142 int a = s - k + 1; 143 int b = s + n25 + (s - k + 1) / 2; 144 add(ans, dp(1, s, k) * res[a][b]); 145 } 146 if ((s + n) % 2 == 1) { 147 res = fpm(x, N - 1); 148 for (int k = 0; k <= s - 1; ++k) 149 if ((s - k) % 2 == 0) { 150 int a = s - k + 1; 151 int b = s + (s - k + 1) / 2; 152 sub(ans, dp(1, s, k) * res[a][b]); 153 } 154 } 155 cout << ans << endl; 156 return 0; 157 }
标签：ng2 ems sizeof clear 还需要 names bsp class 还需
原文地址：http://www.cnblogs.com/ichn/p/6387779.html
