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

问题描述
最近,小栋在无向连通图的生成树个数计算方面有了惊人的进展,他发现:
·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;}