【基于连通性的状态压缩DP】【NOI2007】生成树计数

最新推荐文章于 2020-11-07 18:13:09 发布

Whjpji

最新推荐文章于 2020-11-07 18:13:09 发布

阅读量2.2k

点赞数

分类专栏： OI 文章标签： matrix permutation date 存储优化

本文链接：https://blog.csdn.net/whjpji/article/details/7617796

版权

OI 专栏收录该内容

177 篇文章 0 订阅

订阅专栏

问题描述
最近,小栋在无向连通图的生成树个数计算方面有了惊人的进展,他发现:
·n 个结点的环的生成树个数为 n。
·n 个结点的完全图的生成树个数为 n^(n-2)。
这两个发现让小栋欣喜若狂,由此更加坚定了他继续计算生成树个数的想法,他要计算出各种各样图的生成树数目。
一天,小栋和同学聚会,大家围坐在一张大圆桌周围。小栋看了看,马上想到了生成树问题。如果把每个同学看成一个结点,邻座(结点间距离为 1)的同学间连一条边,就变成了一个环。可是,小栋对环的计数已经十分娴熟且不再感兴趣。于是,小栋又把图变了一下:不仅把邻座的同学之间连一条边,还把相隔一个座位(结点间距离为 2)的同学之间也连一条边,将结点间有边直接相连的这两种情况统称为有边相连,如图 1 所示。

小栋以前没有计算过这类图的生成树个数,但是,他想起了老师讲过的计算任意图的生成树个数的一种通用方法:构造一个 n×n 的矩阵 A={aij}

其中 di 表示结点 i 的度数。
与图 1 相应的 A 矩阵如下所示。为了计算图 1 所对应的生成数的个数,只要去掉矩阵 A 的最后一行和最后一列,得到一个(n-1)×(n-1)的矩阵 B,计算出矩阵B 的行列式的值便可得到图 1 的生成树的个数。

所以生成树的个数为 B = 3528 。小栋发现利用通用方法,因计算过于复杂而很难算出来,而且用其他方法也难以找到更简便的公式进行计算。于是,他将图做了简化,从一个地方将圆桌断开,这样所有的同学形成了一条链,连接距离为 1 和距离为 2 的点。例如八个点的情形如下:

这样生成树的总数就减少了很多。小栋不停的思考,一直到聚会结束,终于找到了一种快捷的方法计算出这个图的生成树个数。可是,如果把距离为 3 的点也连起来,小栋就不知道如何快捷计算了。现在,请你帮助小栋计算这类图的生成树的数目。
输入文件
输入文件中包含两个整数 k, n,由一个空格分隔。k 表示要将所有距离不超过 k(含 k)的结点连接起来,n 表示有 n 个结点。
输出文件
输出文件输出一个整数,表示生成树的个数。由于答案可能比较大,所以你只要输出答案除 65521 的余数即可。
输入样例
35
输出样例
75
样例说明
样例对应的图如下:

朴素模拟代码：

/***************************\
 * @prob: NOI2007 count    *
 * @auth: Wang Junji       *
 * @stat: TLE: 40          *
 * @date: May. 29th, 2012  *
 * @memo: 模拟              *
\***************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <string>
#include <cstring>

const int maxN = 1010, MOD = 65521;
typedef long long int64;
long long a[maxN][maxN], ans;
int ord[maxN], n, K;

int main()
{
    freopen("count.in", "r", stdin);
    freopen("count.out", "w", stdout);
    scanf("%d%d", &K, &n);
    for (int i = 1; i < n; ++i)
    for (int j = 1; j < n; ++j)
    if (i == j) a[i][j] = std::min(n, i + K) - std::max(1, i - K);
    else a[i][j] = (abs(i - j) <= K) ? -1 : 0;
    for (int i = 1; i < n; ++i) ord[i] = i;
    ans = 1;
    for (int i = 1; i < n; ++i) (ans *= a[i][i]) %= MOD;
    for (;;)
    {
        std::next_permutation(ord + 1, ord + n);
        if (ord[1] > K + 1) break;
        bool br = 1, cn = 0;
        for (int i = 1; i < n; ++i)
        if (ord[i] > i + K || ord[i] < i - K)
            {cn = 1; break;}
        if (cn) continue;
        for (int i = 1; i < n; ++i)
            if (ord[i] - i) {br = 0; break;}
        if (br) break;
        long long tmp = 1;
        for (int i = 1; i < n; ++i)
            if (!((tmp *= a[i][ord[i]]) %= MOD)) break;
        do
        {
            if (!tmp) break;
            for (int i = 1; i < n; ++i)
            for (int j = i + 1; j < n; ++j)
                if (ord[i] > ord[j]) tmp = -tmp;
        } while (0);
        (ans += tmp) %= MOD;
    }
    printf("%d\n", (int)((ans + MOD) % MOD));
    return 0;
}

但此题正解并不是用行列式来进行计算，而是基于连通性的状态压缩DP。
用f[i][j]表示前i个结点中，倒数K个结点的连通性为j的方案总数（其中用最小表示法来描述倒数K个结点的连通性）。
首先用一个递归来找出所有的状态。（如K = 3时有111，112，121，122，123五种状态；K = 4时有15种状态；K = 5时有52种状态。）
程序中为了方便，一律使用八进制倒序存储每一个状态。

我们在每一种状态中，要维护出该状态对应的实际情况有多少种，可以用题目中给出的行列式算出。（比如状态111对应3种不同的生成树

；状态1111对应16种不同的生成树

；状态11111对应125种不同的生成树

。）

那么在递推的时候，要保证以下两点：
1) 由当前状态推出的下一状态不能有环；
2) 由当前状态推出的下一状态必须与当前状态所有的点都连通。
而我们可以用并查集来判断并求出所有满足条件的转移数，存储在g数组中, 方便转移。

于是转移的方程为：

，
最后只需要输出f[n][0]（即所有点都连通，由于状态111或1111或11111总是第一个生成出来，序号为1）。
代码：

/******************************\
 * @prob: NOI2007 count       *
 * @auth: Wang Junji          *
 * @stat: RE: 80              *
 * @date: May. 30th, 2012     *
 * @memo: 基于连通性的状压DP     *
\******************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>

typedef long long int64;
const int maxN = 10010, maxS = 60, MOD = 65521;
int64 f[maxN][maxS], n;
int status[maxS], ord[060000], g[maxS][maxS], F[10], K, cnt_s;

inline bool check(int ths)
{
    int tmp = 2;
    for (int i = 3; i < K * 3; i += 3)
    {
        for (int j = 1; j < ((ths >> i) & 7); ++j)
            if (!(tmp & (1 << j))) return 0;
        tmp |= 1 << ((ths >> i) & 7);
    } //用最小表示法的定义来判断状态ths是否合法。
    return 1;
}

void Find(int x, int ths)
{
    if (x >= K)
    {
        if (check(ths)) ord[ths] = cnt_s, status[cnt_s++] = ths;
        return;
    }
    for (int i = 1; i < K + 1; ++i) Find(x + 1, ths | (i << (3 * x)));
    return;
} //递归生成每一种状态。

static int Find(int x) {return F[x] == x ? x : (F[x] = Find(F[x]));}

int main()
{
    freopen("count.in", "r", stdin);
    freopen("count.out", "w", stdout);
    scanf("%d%lld", &K, &n);
    Find(1, 1);
    for (int i = 0; i < cnt_s; ++i)
    {
        f[K][i] = 1;
        static int cnt[10];
        memset(cnt, 0, sizeof cnt);
        for (int j = 0; j < K; ++j) ++cnt[(status[i] >> (j * 3)) & 7];
        for (int j = 1; j < K + 1; ++j)
        {
            if (cnt[j] == 3) f[K][i] = 3;
            if (cnt[j] == 4) f[K][i] = 16;
            if (cnt[j] == 5) f[K][i] = 125;
        }
        int ths = status[i];
        for (int s = 0; s < (1 << K); ++s)
        //01枚举下一个点与当前状态中所有点的连通关系。
        {
            for (int j = 0; j < K + 1; ++j) F[j] = j; //并查集数组初始化。
            for (int j = 0; j < K; ++j)
            for (int k = j + 1; k < K; ++k)
            if (((ths >> (j * 3)) & 7) == ((ths >> (k * 3)) & 7))
                F[Find(j)] = Find(k);
            //将当前状态中的连通关系维护出来。
            bool flag = 1;
            for (int j = 0; j < K; ++j) if (s & (1 << j))
            {
                if (Find(K) - Find(j)) F[F[j]] = F[K];
                else {flag = 0; break;}
            }
            //将下一个点与当前状态中相连的点的连通关系维护出来并排除环。
            if (!flag) continue; flag = 0;
            for (int j = 1; j < K + 1; ++j)
                if (Find(0) == Find(j)) {flag = 1; break;}
            //保证当前状态的第一个点与后面的至少一个点相连，
            //否则整个图将会从此处断开。
            if (!flag) continue;
            int use = 0, nxt = 0;
            for (int j = 0; j < K; ++j)
            if (!(nxt & (7 << (j * 3))))
            {
                nxt |= ++use << (j * 3);
                for (int k = j + 1; k < K; ++k)
                if (Find(j + 1) == Find(k + 1))
                    nxt |= use << (k * 3);
            } //根据最小表示法的规则生成下一个状态。
            ++g[i][ord[nxt]];
            //由于多种不同的连通性可能对应同一种状态，
            //所以要把所有的情况都算进来。
        }
    }
    //找出每个状态能够推出那些状态，并求出转移数
    //（即从i～K+i-1的连通性推出i+1~K+i）
    for (int i = K + 1; i < n + 1; ++i)
    for (int j = 0; j < cnt_s; ++j)
    for (int k = 0; k < cnt_s; ++k)
        if (g[j][k]) (f[i][k] += g[j][k] * f[i - 1][j]) %= MOD;
    printf("%d\n", (int)f[n][0]);
    return 0;
}

直接递推的效率太低！
当我们发现所有转移都是相同的的时候，就应该想到——矩阵！矩阵乘法可以帮助我们快速转移！

就把之前用来转移的g数组转置（交换行列坐标）过后作为基底，将基底自乘n - K次之后在乘以在K处的所有状态即得到了最终的所有状态。
代码：

/********************************************\
 * @prob: NOI2007 count                     *
 * @auth: Wang Junji                        *
 * @stat: Accepted.                         *
 * @date: May. 30th, 2012                   *
 * @memo: 基于连通性的状压DP、矩阵乘法优化       *
\********************************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>

typedef long long int64;
const int maxN = 10010, maxS = 60, MOD = 65521;
int64 vect[maxS], res, n;
int status[maxS], ord[060000], F[10], K, cnt_s;

struct Matrix
{
    int64 a[maxS][maxS];
    Matrix() {memset(a, 0, sizeof a);}
    Matrix &operator=(const Matrix &b)
    {memcpy(a, b.a, sizeof b.a); return *this;}
    Matrix &operator*=(const Matrix &b)
    {
        Matrix ans;
        for (int i = 0; i < cnt_s; ++i)
        for (int j = 0; j < cnt_s; ++j)
        for (int k = 0; k < cnt_s; ++k)
            (ans.a[i][j] += a[i][k] * b.a[k][j]) %= MOD;
        memcpy(a, ans.a, sizeof ans.a);
        return *this;
    }
} g;

inline bool check(int ths)
{
    int tmp = 2;
    for (int i = 3; i < K * 3; i += 3)
    {
        for (int j = 1; j < ((ths >> i) & 7); ++j)
            if (!(tmp & (1 << j))) return 0;
        tmp |= 1 << ((ths >> i) & 7);
    }
    return 1;
}

void Find(int x, int ths)
{
    if (x >= K)
    {
        if (check(ths)) ord[ths] = cnt_s, status[cnt_s++] = ths;
        return;
    }
    for (int i = 1; i < K + 1; ++i) Find(x + 1, ths | (i << (3 * x)));
    return;
}

static int Find(int x) {return F[x] == x ? x : (F[x] = Find(F[x]));}

inline void pow(Matrix &a, int64 n)
{
    Matrix ans;
    for (int i = 0; i < cnt_s; ++i) ans.a[i][i] = 1;
    for (; n; n >>= 1, a *= a) if (n & 1) ans *= a;
    a = ans; return;
}

int main()
{
    freopen("count.in", "r", stdin);
    freopen("count.out", "w", stdout);
    scanf("%d%lld", &K, &n);
    Find(1, 1);
    for (int i = 0; i < cnt_s; ++i)
    {
        vect[i] = 1; static int cnt[10];
        memset(cnt, 0, sizeof cnt);
        for (int j = 0; j < K; ++j) ++cnt[(status[i] >> (j * 3)) & 7];
        for (int j = 1; j < K + 1; ++j)
        {
            if (cnt[j] == 3) vect[i] = 3;
            if (cnt[j] == 4) vect[i] = 16;
            if (cnt[j] == 5) vect[i] = 125;
        }
        int ths = status[i];
        for (int s = 0; s < (1 << K); ++s)
        {
            for (int j = 0; j < K + 1; ++j) F[j] = j;
            for (int j = 0; j < K; ++j)
            for (int k = j + 1; k < K; ++k)
            if (((ths >> (j * 3)) & 7) == ((ths >> (k * 3)) & 7))
                F[Find(j)] = Find(k);
            bool flag = 1;
            for (int j = 0; j < K; ++j) if (s & (1 << j))
            {
                if (Find(K) - Find(j)) F[F[j]] = F[K];
                else {flag = 0; break;}
            }
            if (!flag) continue; flag = 0;
            for (int j = 1; j < K + 1; ++j)
                if (Find(0) == Find(j)) {flag = 1; break;}
            if (!flag) continue;
            int use = 0, nxt = 0;
            for (int j = 0; j < K; ++j) if (!(nxt & (7 << (j * 3))))
            {
                nxt |= ++use << (j * 3);
                for (int k = j + 1; k < K; ++k)
                if (Find(j + 1) == Find(k + 1))
                    nxt |= use << (k * 3);
            }
            ++g.a[ord[nxt]][i];
        }
    }
    pow(g, n - K);
    for (int i = 0; i < cnt_s; ++i)
        (res += g.a[0][i] * vect[i]) %= MOD;
    printf("%d\n", (int)res);
    return 0;
}