Lightning----HDU_4305----生成树的计数----矩阵树定理

题目地址：http://acm.hdu.edu.cn/showproblem.php?pid=4305

题目的大意是给你一些点的坐标，然后有一个距离限制R。如果两点之间的距离小于R且他们之间没有点与他们共线就可以连通。最后要你求连通图的个数。

这个题目让我又学到了一点，那就是用矩阵树定理来计算生成树的个数。

在这里我不就证明展开讨论，因为我证明不来，感兴趣的可以看看周冬《生成树的计数及其应用》。

我就直接说定理就好了。

Matrix-Tree定理是解决生成树计数问题最有力的武器之一。
它首先于1847年被Kirchhoff证明。在介绍定理之前，我们首先明确几个概念：
1、G的度数矩阵D[G]是一个n*n的矩阵，并且满足：当i≠j时,dij=0；当i=j时，dij等于vi的度数。
2、G的邻接矩阵A[G]也是一个n*n的矩阵， 并且满足：如果vi、vj之间有边直接相连，则aij=1，否则为0。
我们定义G的Kirchhoff矩阵(也称为拉普拉斯算子)C[G]为C[G]=D[G]-A[G]，
则Matrix-Tree定理可以描述为：G的所有不同的生成树的个数等于其Kirchhoff矩阵C[G]
任何一个n-1阶主子式的行列式的绝对值。所谓n-1阶主子式，就是对于r(1≤r≤n)，
将C[G]的第r行、第r列同时去掉后得到的新矩阵，用Cr[G]表示

 

那么我们要做的就是根据题目给出的数据来建立相应的矩阵，然后求出这个矩阵的n-1阶行列式的绝对值就OK了。

另外为了优化一下算法，我们可以先DFS，看这个图是否能够连通，如果不能连通，则字节输出-1，后面就没必要在运算了。

但是事实远比我们想象的要复杂。我个人认为这题目最麻烦的就是求行列式。

我采用的是通过高斯消元法将矩阵化成上三角矩阵，然后求对角线的乘积。

我们可以看到，题目中说数据很大，所以我们要进行取模运算。

而且在对矩阵的操作中，如果交换行，则行列式的值就变成了相反数

如果乘以了某个系数，则行列式的值也相应的乘以这个系数。

我的思路来自于cyberzhg的博客。就是用一个系数来记录乘以的系数，最后想办法除掉就行了。

另外在消元的时候我们不能按照以往我们在线性代数中的那样直接除，这样会产生浮点数。所以我们采用乘以逆元的做法来消元。

当然这种做法就不可避免的乘上系数，所以我们就要记录了。

在最后除以系数的时候也有个要注意的地方，在这里再次感谢cyberzhg。下面是他给我的解答。

 

    比如因为初等行变换，行列式的值被扩大了7倍，
    最终算出来行列式的值是1，这时的1明显是因为太大取余后得到的1，
    这时要除以7就会出现小数，但实际的值必然是一个整数。
    假设除以7后应该得到的取余后的数是x，那么x*7取余后一定等于1，
    这部做的就是枚举x，因为x也是取余后的值，不可能大于MOD。

 

所以在我的代码中才会在最后出现了那样的操作。至于换行的话，我们只要记住到底换了几次。如果是奇数次的话，就要取反，也就是用MOD-ans。

网上还有一些其它大神的代码，但是我觉得这个比较容易理解，适合小牛和水牛，比如说我。我就是水牛。下面上代码：

/*
Matrix-Tree定理是解决生成树计数问题最有力的武器之一。
它首先于1847年被Kirchhoff证明。在介绍定理之前，我们首先明确几个概念：
1、G的度数矩阵D[G]是一个n*n的矩阵，并且满足：当i≠j时,dij=0；当i=j时，dij等于vi的度数。
2、G的邻接矩阵A[G]也是一个n*n的矩阵， 并且满足：如果vi、vj之间有边直接相连，则aij=1，否则为0。
我们定义G的Kirchhoff矩阵(也称为拉普拉斯算子)C[G]为C[G]=D[G]-A[G]，
则Matrix-Tree定理可以描述为：G的所有不同的生成树的个数等于其Kirchhoff矩阵C[G]
任何一个n-1阶主子式的行列式的绝对值。所谓n-1阶主子式，就是对于r(1≤r≤n)，
将C[G]的第r行、第r列同时去掉后得到的新矩阵，用Cr[G]表示
*/

#include<iostream>
using namespace std;

#define MOD 10007
#define N 305

int x[N],y[N];
int map[N][N];
bool visit[N];      //用于DFS，如果不能够连通，则可以直接返回-1
int n,r;
bool flag;

int min(int a,int b)
{
    return a<b?a:b;
}

int max(int a,int b)
{
    return a>b?a:b;
}

bool dfs(int x)
{
    if(!visit[x])
    {
        visit[x] = true;
        for(int i=0;i<n;i++)
        {
            if(map[x][i]) //如果x可以到i
                dfs(i);
        }
    }
    for(int i=0;i<n;i++)
        if(!visit[i])
            return false;
    return true;
}

void init()
{
    int i,j,k;
    scanf("%d%d",&n,&r);
    for(i=0;i<n;i++)
    {
        scanf("%d%d",&x[i],&y[i]);
    }
    memset(map,0,sizeof(map));
    memset(visit,false,sizeof(visit));
    r = r*r;

    //先构建这个图的邻接矩阵
    for(i=0;i<n;i++)
    {
        for(j=i+1;j<n;j++)
        {
            if(((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j])) <= r)
            {
                map[i][j] = map[j][i] = 1;
                for(k=0;k<n;k++)
                {
                    //如果存在有点在它们连线的中间，这不能连
                    if(k!=i && k!=j && x[k] > min(x[i],x[j]) && x[k] < max(x[i],x[j])
                       && (x[k]-x[i])*(y[j]-y[i]) == (x[j]-x[i])*(y[k]-y[i]))
                    {
                        map[i][j] = map[j][i] = 0;
                        break;
                    }
                }
            }
        }
    }

    //在建立Kirchhoff矩阵之前，我们得先判断这个是否连通，只有连通了建立这个矩阵才有意义

    flag = dfs(0);
    if(!flag)
        return ;

    //下面由邻接矩阵来构造Kirchhoff矩阵
    for(i=0;i<n;i++)
    {
        for(j=i+1;j<n;j++)
        {
            if(map[i][j]) //如果i,j连通，则度数加一，且其余的为-1,但是我们知道这个最终是要模的明所以可以在-1上加MOD
            {
                ++map[i][i];
                ++map[j][j];
                map[i][j] = map[j][i] = MOD-1;
            }
        }
    }
}

//当把Kirchhoff矩阵建立好之后，我们就要来求n-1阶的行列式了
//思路是利用高斯消元法将矩阵化成上三角
//然后求对角线的积
//但是在化简矩阵的时候，我们会交换以及乘以某些系数
//由于我们求的是绝对值，如果是因为交换行得到了相反数的绝对值，就可以取反加MOD
//我们只需要求乘的系数的乘积
//因为当你乘以X，行列式的值也要相应的乘以X

int Gauss_Det(int row,int clo)
{
    bool flag2 = false;
    int xishu = 1;
    int i=0;
    int j,k;
    int index;
    int result;
    for(j=0;j<clo;j++) //从第0列开始消元
    {
        index = i;
        for(k=i;k<row;k++)  //在第j列找到一个不是0的，然后和i行进行交换
        {
            if(map[k][j] >0 )
            {
                index = k;
                break;
            }
        }
        if(map[index][j])
        {
            if(index != i)
            {
                for(k=j;k<clo;k++)
                {
                    int t = map[i][k];
                    map[i][k] = map[index][k];
                    map[index][k] = t;
                }
                flag2 = !flag2;
            }
            //消元
            for(k=i+1;k<row;k++)
            {
                if(map[k][j]>0)
                {
                    xishu = (xishu*map[i][j])%MOD;
                    for(int l = clo-1;l>=j;l--)  //因为我们还要用到map[k][j]，所以我们在最后更新
                    {
                        map[k][l] = (map[k][l]*map[i][j] - map[i][l]*map[k][j])%MOD; //这里不用除法是为了避免误差
                        if(map[k][l]<0)
                            map[k][l]+=MOD;
                    }
                }
            }
            i++;
        }

    }



    result = 1;
    for(i=0;i<row;i++)
    {
        result = (result*map[i][i])%MOD;
    }
    /*
    比如因为初等行变换，行列式的值被扩大了7倍，
    最终算出来行列式的值是1，这时的1明显是因为太大取余后得到的1，
    这时要除以7就会出现小数，但实际的值必然是一个整数。
    假设除以7后应该得到的取余后的数是x，那么x*7取余后一定等于1，
    这部做的就是枚举x，因为x也是取余后的值，不可能大于MOD。
    */

    for(i=1;i<MOD;++i)
    {
        if((xishu * i) % MOD == result)
        {
            break;
        }
    }
    result = i;

    if(flag2)       //这是在取反后的操作
        result = MOD - result;
    return result;
}

int main()
{
    int t;
    cin>>t;
    while(t--)
    {
        init();
        if(!flag)
            cout<<-1<<endl;
        else
            cout<<Gauss_Det(n-1,n-1)<<endl;
    }
    return 0;
}