最新文章

[SDOI2015]序列统计 题解

[SDOI2015]序列统计 题解

题目地址:洛谷:【P3321】[SDOI2015]序列统计 – 洛谷、BZOJ 

【生物】高中生物选修三:胚胎工程

【生物】高中生物选修三:胚胎工程

[本文适用于2015、2016级人教版生物课本] [不建议读者在毫无基础的情况下阅读此文, 

[USACO10HOL]赶小猪Driving Out the Piggies 题解

[USACO10HOL]赶小猪Driving Out the Piggies 题解

题目地址:洛谷:【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;
}
[USACO09NOV]灯Lights 题解

[USACO09NOV]灯Lights 题解

题目地址:洛谷:【P2962】[USACO09NOV]灯Lights – 洛谷 

[POI2000]病毒 题解

[POI2000]病毒 题解

题目地址:洛谷:【P2444】[POI2000]病毒 – 洛谷、BZOJ:Pr 

[HNOI2008]GT考试 题解

[HNOI2008]GT考试 题解

题目地址:洛谷:【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;
}
[POI2006]OKR-Periods of Words 题解

[POI2006]OKR-Periods of Words 题解

题目地址:洛谷:【P3435】[POI2006]OKR-Periods of Words  

KS的OI模板库开始维护啦~

KS的OI模板库开始维护啦~

最近是省选复习阶段,就想边打板子边整理板子,于是就在GayHub上建了个仓库放模板,也分享 

Lucas定理及其扩展原理及应用

Lucas定理及其扩展原理及应用

Lucas定理(Lucas’s theorem)

内容

对于非负整数m和n和素数p,有
\mathrm{C}_n^m \equiv \prod_{i=0}^k \mathrm{C}_{n_i}^{m_i} \pmod{p}
其中,m = m_kp^k + m_{k-1}p^{k-1} + \cdots + m_1p + m_0n = n_kp^k + n_{k-1}p^{k-1} + \cdots + n_1p + n_0
也就是说,这个是m和n的p进制展开。
m < n时,\mathrm{C}_n^m = 0

求组合数的值

既然我们可以把m和n以p进制展开求值,那么我们可以设计递归来解决这个问题。
我们知道\mathrm{C}_n^0 = 1,那么以m=0为递归终点即可。代码如下。

inline LL lucas(LL n, LL m, LL p) {
    if(!m) return 1;
    return C(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

其中C就是直接求组合数值的函数,这个可以预处理杨辉三角或者预处理阶乘及阶乘逆元。

例题:【P3807】【模板】卢卡斯定理 – 洛谷

直接应用上述方法即可。代码如下。

// Code by KSkun, 2018/3
#include <cstdio>

#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 = 100005;

LL T, n, m, p, inv[MAXN], mul[MAXN];

inline LL C(LL n, LL m, LL p) {
    if(m > n) return 0;
    return mul[n] * inv[n - m] % p * inv[m] % p;
}

inline LL lucas(LL n, LL m, LL p) {
    if(!m) return 1;
    return C(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

inline void init() {
    mul[0] = inv[0] = mul[1] = inv[1] = 1;
    for(int i = 2; i <= std::min(n, p); i++) {
        inv[i] = -(p / i) * inv[p % i] % p;
    }
    for(int i = 2; i <= std::min(n, p); i++) {
        inv[i] = inv[i] * inv[i - 1] % p;
        mul[i] = mul[i - 1] * i % p;
    }
}

int main() {
    T = readint();
    while(T--) {
        n = readint(); m = readint(); p = readint(); n += m;
        init();
        printf("%lld\n", (lucas(n, m, p) + p) % p);
    }
    return 0;
}

扩展Lucas定理(exLucas)

内容

扩展Lucas定理与Lucas定理的区别是模数非质数,此时,我们考虑将模数质因数分解如下。
P = p_1^{k_1} \cdot p_2^{k_2} \cdots p_n^{k_n}
只要我们能求出组合数对于每个p_i^{k_i}的余数,我们就可以利用中国剩余定理(中国剩余定理及其扩展原理及应用 | KSkun’s Blog)解决这个问题。也就是说,我们将问题变为了这个:
\begin{cases} \mathrm{C}_n^m \equiv a_1 \pmod{p_1^{k_1}} \\ \mathrm{C}_n^m \equiv a_2 \pmod{p_2^{k_2}} \\ \vdots \\ \mathrm{C}_n^m \equiv a_n \pmod{p_n^{k_n}} \end{cases}
现在的问题是,上面每个式子的余数怎么求出。
如果我们使用组合数的定义,我们可以分别求出n! \bmod p_i^{k_i}, m! \bmod p_i^{k_i}, (n-m)! \bmod p_i^{k_i}
首先,我们知道阶乘中存在一些p_i因子,对于这些因子我们可以处理出来单独运算。对于提取一个因子后的阶乘,我们观察它们的规律。假如p=3,要求10!,提取后就变成了:
1 \times 2 \times 1 \times 4 \times 5 \times 2 \times 7 \times 8 \times 3 \times 10 \times 3^3
p_i的倍数组成了一个新的阶乘,这一部分我们递归下去处理。考虑剩下的部分,即1 \times 2 \times \cdots \times (p_i^{k_i}-1) \times (p_i^{k_i}+1) \times \cdots \times (2p_i^{k_i}-1) \times \cdots。我们发现\prod_{i=1, i \bmod p_i \neq 0}^{p_i^{k_i}-1} i \equiv \prod_{i=p_i^{k_i}+1, i \bmod p_i \neq 0}^{2p_i^{k_i}-1} i \pmod{p_i^{k_i}},也就是说,只需要求一个长为p_i^{k_i}的阶乘,并且计算有多少完整的一段,使用快速幂即可求出。至于多出去的那段,直接暴力乘进去即可。

例题:[国家集训队]礼物 题解 | KSkun’s Blog

代码参见该题题解文章。代码内有注释,可以对照上述内容理解。

参考资料

[SCOI2012]喵星球上的点名 题解

[SCOI2012]喵星球上的点名 题解

题目地址:洛谷:【P2336】[SCOI2012]喵星球上的点名 – 洛谷、B