矩阵快速幂

计蒜客 Coin 附带矩阵加速模板

Posted on

大概是第5次写矩阵快速幂……
这道题在DP还是矩阵上都是属于简单范畴……,但没想到用DP就很烦……

感觉矩阵加速还是挺有必要学一下的……

顺带学了一下什么叫分数取模用逆元

题意:
让你抛 k 次 硬币,求上面朝上的次数为偶数的概率。
单次正面朝上的概率为 \( \dfrac {q} {p} \)

思路:
令 dp [ n ] [ 0 ] 表示抛了 n 次后,朝上次数为偶数的概率,dp[ n ] [ 1 ] 则表示奇数的概率。
得转移方程
$ dp[n][0] = \dfrac {p-q} {p} \times dp[n-1][0] + \dfrac {q} {p} \times dp[n-1][1] $
$ dp[n][1] = \dfrac {q} {p} \times dp[n-1][0] + \dfrac {p-q} {p} \times dp[n-1][1] $

然后矩阵加速一下就解决了……

#include <bits/stdc++.h>

#define each(i, n) for (int(i) = 0; (i) < (n); (i)++)
#define reach(i, n) for (int(i) = n - 1; (i) >= 0; (i)--)
#define range(i, st, en) for (int(i) = (st); (i) <= (en); (i)++)
#define rrange(i, st, en) for (int(i) = (en); (i) >= (st); (i)--)
#define fill(ary, num) memset((ary), (num), sizeof(ary))

using namespace std;

typedef long long ll;

const int inf = 0x3f3f3f3f;

const int max_n = 5;
const int mod = 1e9 + 7;

struct Matrix {
    int mat[max_n][max_n], n;
    Matrix(int _n = 1)
    {
        n = _n;
        memset(mat, 0, sizeof mat);
    }
    Matrix operator*(const Matrix& a) const
    {
        Matrix tmp(n);
        for (int i = 1; i < n; ++i)
            for (int j = 1; j < n; ++j)
                if (mat[i][j])
                    for (int k = 1; k < n; ++k)
                        tmp.mat[i][k] = (tmp.mat[i][k] + 1LL * mat[i][j] * a.mat[j][k] % mod) % mod;
        return tmp;
    }
    Matrix Pow(ll m)
    {
        Matrix ret(n), a(*this);
        for (int i = 1; i < n; i++)
            ret.mat[i][i] = 1;
        while (m) {
            if (m & 1)
                ret = ret * a;
            a = a * a;
            m >>= 1;
        }
        return ret;
    }
};

ll fastPow(ll n, ll m)
{
    ll ret = 1;
    while (m) {
        if (m & 1)
            ret = ret * n % mod;
        m >>= 1;
        n = n * n % mod;
    }
    return ret;
}

int main()
{
    int T;
    ll p, q, k;
    scanf("%d", &T);
    while (T--) {
        scanf("%lld %lld %lld", &p, &q, &k);
        ll up = q * fastPow(p, mod - 2) % mod;
        ll down = (p - q) * fastPow(p, mod - 2) % mod;
        if (k == 1) {
            printf("%lld\n", up);
            continue;
        }
        Matrix base(3), res(3);
        base.mat[1][1] = down, base.mat[1][2] = up;
        base.mat[2][1] = up, base.mat[2][2] = down;
        res = base.Pow(k);
        printf("%d\n", res.mat[1][1]);
    }
    return 0;
}
Havel

HDU 2454 Degree Sequence of Graph G

Posted on

Havel定理判定简单图

被上下界网络流烦傻了,暂时不想写代码太长的题。

虽然是一道水题,但也是建立在一定理论基础上的。

再科普一下简单图,没有重边没有自环的图就是简单图。区别有向图和无向图,也就是说,有向图的反向边不算重边。

给定一个非负整数序列{dn},若存在一个无向图使得图中各点的度与此序列一一对应,则称此序列可图化。进一步,若图为简单图,则称此序列可简单图化

可图化的判定:

\( \sum_{i=1}^{n}d_{i} mod 2 == 0 \) 。关于具体图的构造,我们可以简单地把奇数度的点配对,剩下的全部搞成自环。

可简单图化的判定(Havel定理):

把序列排成不增序,即 \( d_1 \geq d_2 \geq \ldots \geq d_n \),则d可简单图化当且仅当
\( d' = \left{ d_2-1 , d_3-1 \ldots d_{d_1+1}-1 , d_{d_1+2} , d_{d_1+3} , \ldots d_n \right} \) 可简单图化。简单的说,把d排序后,找出度最大的点( 设度为\( d_1\) ),把它与度次大的d1个点之间连边,然后这个点就可以不管了,一直继续这个过程,直到建出完整的图,或出现负度等明显不合理的情况。

题意:
题目很长……没怎么看,就是给你几个点,再给你每个点的度数,问你能否构成简单图。

思路:
裸题,按照定理模拟一下即可。

#include <bits/stdc++.h>

using namespace std;
const int inf = 0x3f3f3f3f;
const int maxn = 1005;

int deg[maxn];

int main()
{
    int T, n;
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        bool flag = true;
        int sum = 0;
        for (int i = 0; i < n; i++) {
            scanf("%d", deg + i);
            sum += deg[i];
            if (deg[i] >= n)
                flag = false;
        }
        if (sum & 1)
            flag = false;
        if (!flag) {
            puts("no");
            continue;
        }
        for (int i = 0; i < n; i++) {
            sort(deg + i, deg + n, greater<int>());
            int d = deg[i];
            if (d < 0 || i + d >= n) {
                flag = false;
                break;
            }
            for (int j = 1; j <= d; j++)
                deg[i + j]--;
        }
        puts(flag ? "yes" : "no");
    }
    return 0;
}
AC自动机

POJ 2778 DNA Sequence

Posted on

做的第一道AC自动机非模板题,感觉真的是神了。
一些AC自动机的性质自己完全不知道……

废话不多说,直接上题解。

题意:
给你\( n \)种病毒,要你求出长度为 \( m \) 的不含这 \( n \) 中病毒的数量。

思路:
这道题可以转化成图论解决,当然只能是思路而已
对于AC自动机上的每一个结点,他都可以看成是一个字符串的一部分前缀,也可以认为是一个字符串的一部分后缀。
作为一部分后缀,如果我们在一直不含病毒串的同时,长度达到了要求,我们就可以将其认为是一种可行方案。

解题关键 : 离散数学上的矩阵阶乘定理

构建一个矩阵,矩阵值 \( Mat\left[ i , j \right] \) 代表的是 \( i \) 到 \( j \) 的边数。对其求 k 次幂,所得新矩阵为 \( Mat\left[ i , j \right] \) 代表的是 \( i \) 到 \( j \) 长度为 \( k \) 的路径数量。

应用到这题来,其实我们对于AC自动机上的每一个结点,我们只要求出它的邻接矩阵,再通过矩阵快速幂,累加和,就是我们要的答案。

下面是求邻接矩阵的方法与说明

作为一部分前缀,如果这个结点它选择走向了另一个结点从而形成了我们的病毒串,我们认为这是不可行的。
这是很显然的,比如说有病毒串 ATC ,而我 当前结点及其它的所有前结点构成的字符串为 *****AT,如果这时候我们走向了 T 的孩子 C结点 ,这时候就构成了病毒串,我这里暂时称它为病毒结点。

另一个不允许我们跑动的状态是跑 fail 指针的时候。当我们已知我们的fail指针是指向的病毒结点的时候,我们自然不能通过这个fail指针跑向病毒结点。并且我们的当前结点也成为了病毒结点。因为他们的后缀相同。
举个例子就比如 R -> ATC ,T指向了 R -> T ,其中 T 为病毒结点,那么 R -> AT 中的 T 也成为了病毒结点。

初学者 当然包括我啦 有个很自然的困惑就是这样的情况 R -> ATCT , R -> T 病毒结点的时候,R -> ATCT 的第一个 T 能不能有效的被 “感染” 为病毒结点。
答案当然是显然的,因为我们的遍历过程是宽度优先。如果是深度优先,恐怕不行。

AC Code

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>

using namespace std;

typedef long long ll;
const int max_n = 120;
const int max_c = 4;
const int mod = 1e5;
char buf[max_n];

int getChar(char& ch)
{
    if (ch == 'A')
        return 0;
    else if (ch == 'C')
        return 1;
    else if (ch == 'G')
        return 2;
    return 3;
}

struct Matrix {
    int mat[max_n][max_n], n;
    Matrix(int _n = 1)
    {
        n = _n;
        memset(mat, 0, sizeof mat);
    }
    Matrix operator*(const Matrix& a) const
    {
        Matrix tmp(n);
        for (int i = 1; i < n; ++i)
            for (int j = 1; j < n; ++j)
                if (mat[i][j])
                    for (int k = 1; k < n; ++k)
                        tmp.mat[i][k] = (tmp.mat[i][k] + 1LL * mat[i][j] * a.mat[j][k] % mod) % mod;
        return tmp;
    }
};

struct Aho {
    int next[max_n][max_c], fail[max_n], end[max_n];
    int root, size;
    queue<int> que;

    int newNode()
    {
        for (int i = 0; i < max_c; i++)
            next[size][i] = 0;
        end[size++] = 0;
        return size - 1;
    }

    inline void init()
    {
        size = 1;
        root = newNode();
    }

    inline void insert(char str[])
    {
        int len = strlen(str), now = root;
        for (int i = 0; i < len; i++) {
            int c = getChar(str[i]);
            if (!next[now][c])
                next[now][c] = newNode();
            now = next[now][c];
        }
        end[now] = 1;
    }

    void build()
    {
        fail[root] = root;
        for (int i = 0; i < max_c; i++)
            if (!next[root][i])
                next[root][i] = root;
            else {
                fail[next[root][i]] = root;
                que.push(next[root][i]);
            }
        while (!que.empty()) {
            int now = que.front();
            que.pop();
            if (end[fail[now]])
                end[now] = 1;
            for (int i = 0; i < max_c; i++)
                if (!next[now][i])
                    next[now][i] = next[fail[now]][i];
                else {
                    fail[next[now][i]] = next[fail[now]][i];
                    que.push(next[now][i]);
                }
        }
    }

    Matrix getMatrix()
    {
        Matrix a(size);
        for (int i = 1; i < size; i++)
            for (int j = 0; j < 4; j++)
                if (!end[next[i][j]])
                    a.mat[i][next[i][j]]++;
        return a;
    }
} aho;

Matrix Pow(Matrix a, ll m)
{
    Matrix ans(a.n);
    for (int i = 1; i < a.n; i++)
        ans.mat[i][i] = 1;
    while (m) {
        if (m & 1)
            ans = ans * a;
        a = a * a;
        m >>= 1;
    }
    return ans;
}

int main()
{
    int n;
    ll m;
    while (scanf("%d %lld", &n, &m) != EOF) {
        aho.init();
        for (int i = 0; i < n; i++) {
            scanf("%s", buf);
            aho.insert(buf);
        }
        aho.build();
        Matrix a = aho.getMatrix();
        a = Pow(a, m);
        int res = 0;
        for (int i = 1; i < a.n; i++)
            res = (res + a.mat[1][i]) % mod;
        printf("%d\n", res);
    }
    return 0;
}
Cayley

BZOJ 4766 文艺计算姬

Posted on

完全二分图的生成树计数问题。

这里是作为找规律的题目出现的……
对于这张图的一边 n 个节点,另一边 m 个节点。 其结果为

$$n^{m-1} \times m^{n-1}$$

因为数据非常大,为 1e18 ,long long 也会爆炸
用到了快速乘,快速幂

题意:
计算完全二分图的生成树数量

思路:
公式已给。证明无。
……

AC Code

#include <cstdio>

using namespace std;
typedef unsigned long long ull;

ull fastMul(ull a, ull b, ull mod)
{
    ull ret = 0;
    for (; a; a >>= 1, b = (b + b) % mod)
        if (a & 1)
            ret = (ret + b) % mod;
    return ret;
}

ull fastPow(ull x, ull a, ull mod)
{
    ull ret = 1;
    for (; a; a >>= 1, x = fastMul(x, x, mod))
        if (a & 1)
            ret = fastMul(ret, x, mod);
    return ret;
}
int main()
{
    ull n, m, p;
    scanf("%lld %lld %lld", &n, &m, &p);
    printf("%lld\n", fastMul(fastPow(n, m - 1, p), fastPow(m, n - 1, p), p));
    return 0;
}
Cayley

BZOJ 1430 小猴打架

Posted on

在搜最小生成树计数的题目的时候不小心发现了一个有趣的东西

51nod 1601

完全图的最小生成树计数问题……题解说需要01字典树,不是很会,先放一放
然后就找到了这个资料 Prüfer编码与Cayley公式

简单说完全图的生成树数量为 \( n^{n-2} \)
这个运用矩阵树也可以计算的

题意:
n只猴子打架,两只猴子打架后就会成为好朋友,就不会打架,朋友关系可传递,问打架的序列方案数。

思路:
n个点的生成树数量为 \( n^{n-2} \) ,顺序为 (n-1)!

AC Code

#include <cstdio>
#define each(i, n) for (int(i) = 0; (i) < (n); (i)++)
#define reach(i, n) for (int(i) = n - 1; (i) >= 0; (i)--)
#define range(i, st, en) for (int(i) = (st); (i) <= (en); (i)++)
#define rrange(i, st, en) for (int(i) = (en); (i) >= (st); (i)--)
const int mod = 9999991;
int n;
long long ans = 1;
int main()
{
    scanf("%d", &n);
    range(i, 1, n - 2) ans = ans * n % mod;
    range(i, 1, n - 1) ans = ans * i % mod;
    printf("%lld", ans);
    return 0;
}

Matrix-Tree

BZOJ 4031 小Z的房间

Posted on

矩阵树定理第二题,用于有取模处理的运算。
中间一直超时,还一位bzoj会卡宏定义……真是石乐志

题意:
给你一个n×m矩阵,矩阵中的' . '表示成一个房间,' * ‘ 表示柱子,在相邻的房间中可以连边,要求每个房间都相连,求方案数。

思路:
算是比较裸的矩阵树定理应用了,对每一个点,标号,再对相邻都是点的格子建边。

AC Code

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <queue>

using namespace std;

#define each(i, n) for (int(i) = 0; (i) < (n); (i)++)
#define reach(i, n) for (int(i) = n - 1; (i) >= 0; (i)--)
#define range(i, st, en) for (int(i) = (st); (i) <= (en); (i)++)
#define rrange(i, st, en) for (int(i) = (en); (i) >= (st); (i)--)

typedef long long ll;
const int maxn = 105;
const int mod = 1e9;

ll c[maxn][maxn];
char mat[maxn][maxn];
int G[maxn][maxn];
int dx[] = { -1, 0, 1, 0 };
int dy[] = { 0, -1, 0, 1 };

ll getDet(ll a[][maxn], int n)
{
    range(i, 1, n) range(j, 1, n) a[i][j] = (a[i][j] + mod) % mod;
    ll ret = 1;
    range(i, 1, n - 1)
    {
        range(j, i + 1, n - 1) while (a[j][i])
        {
            ll t = a[i][i] / a[j][i];
            range(k, i, n - 1) a[i][k] = (a[i][k] - a[j][k] * t % mod + mod) % mod;
            range(k, i, n - 1) swap(a[i][k], a[j][k]);
            ret = -ret;
        }
        if (a[i][i] == 0)
            return 0;
        ret = ret * a[i][i] % mod;
    }
    return (ret + mod) % mod;
}

int main()
{
    int n, m, num = 0;
    scanf("%d %d", &n, &m);
    each(i, n) scanf("%s", mat[i]);
    each(i, n) each(j, m) if (mat[i][j] == '.') G[i][j] = ++num;
    each(i, n) each(j, m) if (G[i][j])
    {
        int u = G[i][j];
        each(k, 4)
        {
            int x = i + dx[k], y = j + dy[k];
            if (x < 0 || x >= n || y < 0 || y >= m || mat[x][y] != '.')
                continue;
            int v = G[x][y];
            c[u][u]++, c[u][v]--;
        }
    }
    printf("%lld\n", getDet(c, num));
    return 0;
}
Maths

SPOJ 104 Highways

Posted on

早上在学最小生成树计数。
最小生成树算法基于 Matrix-Tree 定理,Matrix-Tree定理是被广泛应用与求生成树计数中。
其具体流程为构造基尔霍夫矩阵,再对其求 n-1 阶行列式即可。

  • 基尔霍夫矩阵 为 其度数矩阵D[G] - 邻接矩阵A[G]
  • 度数矩阵D[G] 满足,当 i≠ j 时,dij=0;当i=j时,dij等于vi的度数。
  • 邻接矩阵A[G] 满足,如果vi、vj之间有边直接相连,则aij=1,否则为0。

不得不说SPOJ的数据之强,一个小问题找了我大半天才找着,SPOJ倒是我想做验证模板的不二之选。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <queue>
#define ll long long

using namespace std;

const int inf = 0x3f3f3f3f;
const int maxn = 105;

ll c[maxn][maxn];
ll degree[maxn];

inline bool scan_d(int& num)
{
    char in;
    in = getchar();
    if (in == EOF)
        return false;
    while (in < '0' || in > '9')
        in = getchar();
    num = in - '0';
    while (in = getchar(), in >= '0' && in <= '9')
        num *= 10, num += in - '0';
    return true;
}

ll getDet(ll a[][maxn], int n)
{
    ll ret = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = i + 1; j <= n; j++)
            while (a[j][i]) {
                ll t = a[i][i] / a[j][i];
                for (int k = i; k <= n; k++)
                    a[i][k] = (a[i][k] - a[j][k] * t);
                for (int k = i; k <= n; k++)
                    swap(a[i][k], a[j][k]);
                ret = -ret;
            }
        if (a[i][i] == 0)
            return 0;
        ret = ret * a[i][i];
    }
    return ret < 0 ? -ret : ret;
}

int main()
{
    int T, u, v, n, m;
    scanf("%d", &T);
    while (T--) {
        scanf("%d %d", &n, &m);
        memset(c, 0, sizeof c);
        memset(degree, 0, sizeof degree);
        while (m--) {
            scan_d(u);
            scan_d(v);
            c[u][v] = c[v][u] = -1;
            degree[u]++, degree[v]++;
        }
        for (int i = 1; i <= n; i++)
            c[i][i] = degree[i]; //因为不存在自环呀
        printf("%lld\n", getDet(c, n - 1));
    }
    return 0;
}