[SDOI2015]序列统计 题解
题目地址:洛谷:【P3321】[SDOI2015]序列统计 – 洛谷、BZOJ …
May all the beauty be blessed.
题目地址:洛谷:【P2973】[USACO10HOL]赶小猪Driving Out the Piggi… – 洛谷、BZOJ:Problem 1778. — [Usaco2010 Hol]Dotp 驱逐猪猡
一个无向图,节点1有一个炸弹,在每个单位时间内,有p/q的概率在这个节点炸掉,有1-p/q的概率随机选择一条出去的路到其他的节点上。问最终炸弹在每个节点上爆炸的概率。
输入格式:
* Line 1: Four space separated integers: N, M, P, and Q
* Lines 2..M+1: Line i+1 describes a road with two space separated integers: A_j and B_j
输出格式:
* Lines 1..N: On line i, print the probability that city i will be destroyed as a floating point number. An answer with an absolute error of at most 10^-6 will be accepted (note that you should output at least 6 decimal places for this to take effect).
输入样例#1:
2 1 1 2 1 2
输出样例#1:
0.666666667 0.333333333
我们考虑每次到达每个点爆炸是等可能的,因此在每个点爆炸的概率之比等于到达每个点的期望次数之比。而到达每个点的期望次数可以通过以下式子算出
\mathrm{E}(u) = \sum_{(u, v) \in E} \frac{\mathrm{E}(v) \times (1 - p / q)}{\mathrm{degree}[v]}
由于原图不一定是树,无法直接转移,我们考虑高斯消元的解法。对于除了1的节点,到达该点的期望次数等于右边这个求和,而对于节点1,由于一开始就在1位置,次数多1,方程等号右边就是1。
我们求出的是期望次数,和不为1,而概率和为1,需要加起来除一下才能得到答案。
// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>
#include <algorithm>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
const int MAXN = 305, EPS = 1e-10;
int n, m, u, v, p, q, gra[MAXN][MAXN], deg[MAXN];
double mat[MAXN][MAXN];
inline void gauss() {
for(int i = 1; i <= n; i++) {
int r = i;
for(int j = i + 1; j <= n; j++) {
if(std::fabs(mat[r][i]) < std::fabs(mat[j][i])) r = j;
}
if(r != i) {
for(int j = 1; j <= n + 1; j++) {
std::swap(mat[r][j], mat[i][j]);
}
}
for(int j = 1; j <= n; j++) {
if(j != i) {
double t = mat[j][i] / mat[i][i];
for(int k = i + 1; k <= n + 1; k++) {
mat[j][k] -= mat[i][k] * t;
}
}
}
}
for(int i = 1; i <= n; i++) {
mat[i][n + 1] /= mat[i][i];
}
}
int main() {
n = readint(); m = readint(); p = readint(); q = readint();
while(m--) {
u = readint(); v = readint();
gra[u][v] = gra[v][u] = 1;
deg[u]++; deg[v]++;
}
for(int i = 1; i <= n; i++) {
mat[i][i] = 1;
for(int j = 1; j <= n; j++) {
if(gra[i][j]) mat[j][i] = -(1 - double(p) / q) / deg[i];
}
}
mat[1][n + 1] = 1;
gauss();
double sum = 0;
for(int i = 1; i <= n; i++) {
sum += mat[i][n + 1];
}
for(int i = 1; i <= n; i++) {
printf("%.9lf\n", mat[i][n + 1] / sum);
}
return 0;
}
题目地址:洛谷:【P3193】[HNOI2008]GT考试 – 洛谷、BZOJ:Problem 1009. — [HNOI2008]GT考试
阿申准备报名参加GT考试,准考证号为N位数X1X2….Xn(0<=Xi<=9),他不希望准考证号上出现不吉利的数字。他的不吉利数学A1A2…Am(0<=Ai<=9)有M位,不出现是指X1X2…Xn中没有恰好一段等于A1A2…Am. A1和X1可以为0
输入格式:
第一行输入N,M,K.接下来一行输入M位的数。 N<=10^9,M<=20,K<=1000
输出格式:
阿申想知道不出现不吉利数字的号码有多少种,输出模K取余的结果.
输入样例#1:
4 3 100 111
输出样例#1:
81
首先计数问题要考虑动态规划,我们考虑枚举当前串后缀有几位是不吉利数字的前缀,也就是说dp[i][j]表示到i位,后缀有j位是不吉利数字前缀的方案数,注意这里的j就是指匹配上的最大长度,否则可能会计算重复。那么答案就是\sum_{i=0}^{m-1} dp[n][i]。
怎么转移呢?枚举后一位是什么数字然后看看它是不是不吉利数字接上去的那个?这种想法是有问题的,比如对于不吉利数字1312
,如果我们枚举到j=3,而后面接的是3,虽然不是不吉利数字第4位的2,却可以构成不吉利数字的前两位13
这个前缀。也就是说,枚举为3的时候实际上是向j=2转移的。我们令a[i][j]表示不吉利数字i位后缀转移到j位后缀的方案数,如果有了这个数组,我们的转移就可以一层一层来做了。
dp[i][j] = \sum_{k=0}^{m-1} dp[i-1][k] \cdot a[k][j]
初始值是dp[0][0] = 1。
接下来考虑这个a数组怎么求。我们发现其实它可以通过自我匹配来完成。我们考虑利用KMP算法,计算出fail数组后,枚举i位前缀后面接一个什么数,然后计算出此时的最长的与某一前缀相同的后缀长度j,此时i可以向j转移。
现在我们可以在O(m)时间内求a数组了,但是DP仍然是困难的,因为n的范围达到了惊人的1e9,这怎么做?
我们发现其实a数组是个矩阵,而且每次转移利用的a数组都相同,这是一个向量的线性变换呀!马上想到矩阵快速幂,计算转移矩阵a的n次幂,往初始向量上一乘,就得到了答案。这个转移是O(m^2 \log n)的。
// Code by KSkun, 2018/3
#include <cstdio>
#include <cstring>
#include <cctype>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF
: *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
inline char getdigit() {
char c;
while(!isdigit(c = fgc()));
return c;
}
const int MAXN = 30;
int n, m, p, fail[MAXN];
char str[MAXN];
struct Matrix {
LL a[MAXN][MAXN];
Matrix() {
memset(a, 0, sizeof(a));
}
inline Matrix operator*(const Matrix &rhs) const {
Matrix res;
for(int i = 0; i < m; i++) {
for(int j = 0; j < m; j++) {
for(int k = 0; k < m; k++) {
res.a[i][j] = (res.a[i][j] + a[i][k] * rhs.a[k][j]) % p;
}
}
}
return res;
}
inline Matrix& operator*=(const Matrix &x) {
return *this = *this * x;
}
};
Matrix t;
inline Matrix makeunit(int n) {
Matrix res;
for(int i = 0; i < n; i++) res.a[i][i] = 1;
return res;
}
inline Matrix fpow(Matrix n, int k) {
Matrix res = makeunit(m);
while(k) {
if(k & 1) res *= n;
n *= n;
k >>= 1;
}
return res;
}
inline void calfail() {
int i = 2, j = 0;
fail[1] = 0;
for(; i <= m; i++) {
while(j && str[j + 1] != str[i]) j = fail[j];
if(str[j + 1] == str[i]) j++;
fail[i] = j;
}
}
inline void kmp() {
for(int i = 0; i < m; i++) {
for(int j = 0; j <= 9; j++) {
int k = i;
while(k && str[k + 1] != j + '0') k = fail[k];
if(str[k + 1] == j + '0') k++;
if(k != m) t.a[i][k]++;
}
}
}
int main() {
n = readint(); m = readint(); p = readint();
for(int i = 1; i <= m; i++) {
str[i] = getdigit();
}
calfail();
kmp();
t = fpow(t, n);
Matrix res;
res.a[0][0] = 1;
res *= t;
LL ans = 0;
for(int i = 0; i < m; i++) {
ans = (ans + res.a[0][i]) % p;
}
printf("%lld", ans);
return 0;
}