[BZOJ 3640]JC的小苹果
逆矩阵文明,,,
很显然我们可以定义一个状态\(f[i][j]\)表示当前血量为\(i\)走到\(j\)的概率,然后肥肠爆歉这个东西没法DP(可能会有的点的伤害为0,这样可以凿出来环)。考虑高消,这个东西有个好处是不可能从血量低的到血量高的状态,所以可以从大到小枚举血量,这样各层事独立的,复杂度比直接高消降低了很少。可惜复杂度为\(O(sn^3)\)(设血量为\(s\)),会T掉。
考虑转移的过程,转移时等价于解这样一个方程:
\[
\begin{aligned}
a_{11}x_{1} + a_{12}x_{2} + \ldots + a_{1n}x_{n} &= c_1\\
a_{21}x_{1} + a_{22}x_{2} + \ldots + a_{2n}x_{n} &= c_2\\
&\ldots\\
a_{n1}x_{1} + a_{n2}x_{2} + \ldots + a_{nn}x_{n} &= c_n
\end{aligned}
\]
其中的未知数\(x\)事我们要求的东西,\(c\)表示从高血量状态转移过来的概率(这个可以视作常数)。根据友矩阵那一套理论,这一系列方程等价于以下等式:
\[
\begin{bmatrix}
a_{11} & a_{12} & \ldots & a_{1n}\\
a_{21} & a_{22} & \ldots & a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n1} & a_{n2} & \ldots & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_{1}\\ x_{2} \\ \vdots \\ x_{n}
\end{bmatrix}
=
\begin{bmatrix}
c_1\\ c_2\\ \vdots\\ c_n
\end{bmatrix}
\]
不妨将之简记为\(AB = C\),其中我们只有\(B\)不知道,要求的也是\(B\)。我们除一下就可以得到\(B = A^{-1}C\),然后我们还发现每一层的\(A\)都是一样的,所以每一层的\(A^{-1}\)也都是一样的,预处理即可。这样转移部分的复杂度就变成了\(O(sn^2)\)(矩阵乘法在这里事方阵乘列向量)。
至于逆矩阵的求法,我们知道对矩阵做初等变化也就等价于乘上另一个矩阵。因此,我们将一个矩阵\(A\)用类似于高消的手段消为单位阵\(I\),所做的初等变换也就等价于乘上\(A^{-1}\)。我们对一个单位阵\(I\)作用上一样的操作,也就等于给这个单位阵乘上了\(A^{-1}\),这样我们就得到了\(A^{-1}\)。
代码:
/************************************************************** Problem: 3640 User: danihao123 Language: C++ Result: Accepted Time:8708 ms Memory:13416 kb ****************************************************************/ #include <cstdio> #include <cstring> #include <cstdlib> #include <cassert> #include <cmath> #include <algorithm> #include <functional> #include <utility> #include <vector> #include <queue> #include <set> #include <map> int n, m, hp; const int maxn = 151, maxm = 5005; const int maxh = 10005; typedef double R; typedef R Mat[maxn][maxn]; Mat D, F; void calc_inv() { for(int i = 1; i <= n; i ++) { F[i][i] = 1; } for(int i = 1; i <= n; i ++) { int r = i; for(int j = i + 1; j <= n; j ++) { if(fabs(D[j][i]) > fabs(D[r][i])) { r = j; } } // assert(r > -1); if(r != i) { for(int j = 1; j <= n; j ++) { std::swap(D[r][j], D[i][j]); std::swap(F[r][j], F[i][j]); } } R bs = D[i][i]; for(int j = 1; j <= n; j ++) { D[i][j] /= bs; F[i][j] /= bs; } for(int k = 1; k <= n; k ++) { if(k != i) { R rate = D[k][i]; for(int j = 1; j <= n; j ++) { D[k][j] -= rate * D[i][j]; F[k][j] -= rate * F[i][j]; } } } } } void matrix_mul(const Mat &A, const Mat &B, int a, int b, int c, Mat &res) { static Mat C; for(int i = 1; i <= a; i ++) { for(int j = 1; j <= c; j ++) { C[i][j] = 0; } } for(int i = 1; i <= a; i ++) { for(int j = 1; j <= c; j ++) { for(int k = 1; k <= b; k ++) { C[i][j] += A[i][k] * B[k][j]; } } } for(int i = 1; i <= a; i ++) { for(int j = 1; j <= c; j ++) { res[i][j] = C[i][j]; } } } int first[maxn], deg[maxn]; int next[maxm << 1], to[maxm << 1]; int gcnt = 0; void add_edge(int u, int v) { gcnt ++; next[gcnt] = first[u]; first[u] = gcnt; to[gcnt] = v; } void ins_edge(int u, int v) { deg[u] ++; add_edge(u, v); if(u != v) { deg[v] ++; add_edge(v, u); } } int atk[maxn]; R f[maxh][maxn]; R solve() { for(int i = 1; i <= n; i ++) { D[i][i] = 1.00; } for(int i = 1; i < n; i ++) { for(int j = first[i]; j; j = next[j]) { int v = to[j]; if(!atk[v]) { D[v][i] -= 1.00 / (R(deg[i])); } } } calc_inv(); R ans = 0; static Mat C; f[hp][1] = 1.00; for(int i = hp; i >= 1; i --) { for(int j = 1; j <= n; j ++) { C[j][1] = f[i][j]; } matrix_mul(F, C, n, n, 1, C); for(int j = 1; j < n; j ++) { for(int e = first[j]; e; e = next[e]) { int v = to[e]; if(atk[v] > 0 && i - atk[v] > 0) { f[i - atk[v]][v] += C[j][1] / (R(deg[j])); } } } ans += C[n][1]; } return ans; } int main() { scanf("%d%d%d", &n, &m, &hp); for(int i = 1; i <= n; i ++) { scanf("%d", &atk[i]); } for(int i = 1; i <= m; i ++) { int u, v; scanf("%d%d", &u, &v); ins_edge(u, v); } printf("%.8lf\n", solve()); return 0; }
[BZOJ 3328]PYXFIB
又学了个新套路呢qwq
题面要求的其实是:
\[\sum_{i = 0}^n [k|i]\binom{n}{i} F_i\]
不考虑那个\([k|i]\),式子是非常像一个二项式展开的……
考虑构造矩阵的二项式展开,可以发现(其中\(A\)是斐波那契数列的二项式展开):
\[(I + A)^n = \sum_{i = 0}^n \binom{n}{i} A^i\]
然后\(k=1\)的情况我们就会做辣!接下来的问题是,如何取出一个生成函数中次数为\(k\)倍数的项求和?
然后我们发现单位根有个非常优良的性质:
\[\frac{1}{k} \sum_{i=1}^{k} \xi_k^{ni} = [k|n]\]
这里的\(n\)是一个常数,但是其实可以对应到原式里的某一项的次数。至于证明……满足\(k|n\)的情况左边显然是一堆1取平均值,其他情况可以通过等比数列求和证明左边为0。
于是可以构造多项式\((I+xA)^n\),把所有\(k\)次单位根做为\(x\)(求这个的话,可以利用原根)带入,对结果取平均值即可。
代码:
/************************************************************** Problem: 3328 User: danihao123 Language: C++ Result: Accepted Time:9080 ms Memory:832 kb ****************************************************************/ #include <cstdio> #include <cstring> #include <cstdlib> #include <algorithm> #include <utility> #include <cmath> #include <vector> typedef long long ll; typedef ll Mat[2][2]; ll p; int k; ll pow_mod(ll a, ll b) { ll ans = 1LL, res = a % p; while(b) { if(1LL & b) ans = (ans * res) % p; res = (res * res) % p; b /= 2LL; } return ans; } ll inv(ll x) { return pow_mod(x, p - 2LL); } void factor(int x, std::vector<int> &V) { int bd = sqrt(x + 0.5); for(int i = 2; i <= bd; i ++) { if(x % i == 0) { V.push_back(i); while(x % i == 0) x /= i; } } if(x > 1) V.push_back(x); } int get_phi() { if(p == 2LL) return 1LL; std::vector<int> V; factor(p - 1LL, V); for(int i = 2; i <= (p - 1); i ++) { bool flag = true; for(int j = 0; j < V.size(); j ++) { int up = (p - 1) / V[j]; if(pow_mod(i, up) == 1LL) { flag = false; break; } } #ifdef LOCAL if(flag) printf("xi : %d\n", i); #endif if(flag) return i; } } void mat_mul(const Mat &A, const Mat &B, int n, int m, int t, Mat &C) { Mat D; memset(D, 0, sizeof(D)); for(int i = 0; i < n; i ++) { for(int j = 0; j < t; j ++) { for(int k = 0; k < m; k ++) { D[i][j] = (D[i][j] + (A[i][k] * B[k][j]) % p) % p; } } } memcpy(C, D, sizeof(D)); } void mat_pow(const Mat &A, int n, ll b, Mat &ret) { Mat C, res; memset(C, 0, sizeof(C)); for(int i = 0; i < n; i ++) C[i][i] = 1LL; memset(res, 0, sizeof(res)); memcpy(res, A, sizeof(res)); while(b) { if(1LL & b) mat_mul(C, res, n, n, n, C); mat_mul(res, res, n, n, n, res); b /= 2LL; } memcpy(ret, C, sizeof(C)); } int main() { int T; scanf("%d", &T); while(T --) { ll n; scanf("%lld%d%lld", &n, &k, &p); ll xi = pow_mod(get_phi(), (p - 1) / k); ll w = 1LL; ll ans = 0; for(int i = 0; i < k; i ++) { Mat X; memset(X, 0, sizeof(X)); X[0][0] = X[1][0] = X[0][1] = w; X[0][0] = (X[0][0] + 1LL) % p; X[1][1] = (X[1][1] + 1LL) % p; mat_pow(X, 2, n, X); #ifdef LOCAL puts("X : "); for(int i = 0; i < 2; i ++) { for(int j = 0; j < 2; j ++) { printf("%lld ", X[i][j]); } puts(""); } #endif ans = (ans + X[0][0]) % p; w = (w * xi) % p; } ans = (ans * inv(k)) % p; printf("%lld\n", ans); } return 0; }
[BZOJ 1009][HNOI2008]GT考试
好劲啊这题……
首先先想DP的做法,我们用[tex]p[i][j][/tex]表示若字符串已匹配前缀[tex][1,i][/tex],有多少种方案使得追加上一个数后匹配[tex][1,j][/tex],这个用KMP可以很方便的求出(事实上暴力也可以)。
然后再来看DP的做法,转移方程大致是这样的:
[tex]d[i][j]=\Sigma_{k=0}^{m-1}(d[i-1][k]\times p[k][j])[/tex]
但是n非常大,这么做注定要TLE。
不过看这个转移方程,是不是和矩阵乘法的定义很像?是的。
我们可以用一个[tex]1\times m[/tex]的矩阵[tex]D[/tex]来表示某一个阶段的DP值,我们可以把[tex]p[/tex]组织成一个[tex]m\times m[/tex]矩阵[tex]P[/tex]。然后我们可以发现:
[tex]D_i\times P=D_{i+1}[/tex]
由于矩阵乘法满足结合律,所以:
[tex]D_n=D_0\times P^n[/tex]
很明显可以快速幂搞出来。
代码:
/************************************************************** Problem: 1009 User: danihao123 Language: C++ Result: Accepted Time:80 ms Memory:820 kb ****************************************************************/ #include <cstdio> #include <cstring> #include <algorithm> using namespace std; typedef unsigned long long ull; typedef ull Matrix[22][22]; ull MOD; #define REP(i,n) for(i=0;i<(n);i++) #define CL_ARR(x,v) memset(x,v,sizeof(x)) #define CP_ARR(from,to) memcpy(to,from,sizeof(from)) inline void MatrixMul(Matrix A,Matrix B,int m,int n,int p,Matrix& res){ register int i,j,k; Matrix C; CL_ARR(C,0); REP(i,m){ REP(j,p){ REP(k,n){ C[i][j]=(C[i][j]+(A[i][k]*B[k][j])%MOD)%MOD; } } } CP_ARR(C,res); } void MatrixPow(Matrix A,int m,ull p,Matrix& res){ Matrix a,temp; register int i; CL_ARR(a,0); memcpy(a,A,sizeof(a)); CL_ARR(temp,0); REP(i,m){ temp[i][i]=1; } while(p){ if(1&p){ MatrixMul(temp,a,m,m,m,temp); } p>>=1; MatrixMul(a,a,m,m,m,a); } CP_ARR(temp,res); } char buf[25]; int f[25]; int main(){ register int i,u,j; register ull ret=0; ull n; int m; Matrix ans,P; scanf("%llu%d%llu",&n,&m,&MOD); scanf("%s",buf+1); CL_ARR(ans,0); ans[0][0]=1; CL_ARR(P,0); P[0][0]=9;P[0][1]=1; // f[0]=f[1]=0; for(i=1;i<m;i++){ u=f[i]; while(u && buf[u+1]!=buf[i+1]){ u=f[u]; } if(buf[u+1]==buf[i+1]){ f[i+1]=u+1; }else{ f[i+1]=0; } for(j=0;j<10;j++){ u=i; while(u && buf[u+1]!=j+'0'){ u=f[u]; } if(buf[u+1]==j+'0'){ P[i][u+1]=(P[i][u+1]+1)%MOD; }else{ P[i][0]=(P[i][0]+1)%MOD; } } } MatrixPow(P,m,n,P); MatrixMul(ans,P,1,m,m,ans); for(i=0;i<m;i++){ ret=(ret+ans[0][i])%MOD; } printf("%llu\n",ret); return 0; }