HDU 5607 graph（概率DP+矩阵快速幂）

Description
在一个N个点(标号1~n),M条边的有向图上,一开始我在点u,每一步我会在当前点的出边中等概率的选一条走过去,求走了恰好K步后走到每个点的概率.
Input
第一行两个正整数N,M表示点数和边数.
接下来M行,每行两个正整数X,Y表示一条X向Y的一条有向边(保证没有重边和自环).
接下来一个正整数Q,表示询问个数.
接下来Q行,每行两个正整数u,K表示开始的点和步数.
N≤50,M≤1000,Q≤20,u≤n,K≤10^9
每个点保证至少有一个出边
Output
Q行,每行N个数字,用空格隔开,第i个数字表示从u开始走K步到i的概率
考虑到输出的答案可能会有精度问题,经过一定的分析后可以发现答案一定可以被表示成X/Y的形式,你只需输出X×Y^(10^9+5) mod (10^9+7)的值即可.在每行后面多输出一个空格,否则可能会使你PE.
Sample Input
3 2
1 2
1 3
1
1 1
Sample Output
0 500000004 500000004
Solution
用dp[i][j][step]表示从i开始走step步到达j的概率，那么轻易得到以下转移方程
dp[i][j][step]=sum{dp[j][k][step-1]/cnt[i],i->j有边}，其中cnt[i]为i的出度
由于步数是10^9级别，故要用矩阵快速幂加速，将dp[i][j]看作一个矩阵，那么每次转移就是将这个矩阵左乘一个矩阵A，其中A[i][j]=1/cnt[i] (当i到j有边时) or 0(当i到j无边时)，所以对于每次查询u,k,只需要输出A^k中第u行即可
Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
#define maxn 55
#define mod 1000000007ll
struct Mat
{
    ll mat[maxn][maxn];//矩阵 
    int row,col;//矩阵行列数 
};
Mat mod_mul(Mat a,Mat b,ll p)//矩阵乘法 
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i=0;i<ans.row;i++)      
        for(int k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(int j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat Mod_pow(Mat a,ll k,ll p)//矩阵快速幂 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(int i=0;i<a.row;i++)
        for(int j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
ll mod_pow(ll a,ll b,ll p)
{
    a%=p;
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
int n,m,q,u,cnt[maxn];
ll k;
Mat A;
int main()
{
    while(~scanf("%d%d",&n,&m))
    {
        memset(cnt,0,sizeof(cnt));
        A.row=A.col=n;
        while(m--)
        {
            int a,b;
            scanf("%d%d",&a,&b);
            a--,b--;
            cnt[a]++;
            A.mat[a][b]=1;
        }
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                A.mat[i][j]=A.mat[i][j]*mod_pow(cnt[i],mod-2,mod)%mod;
        scanf("%d",&q);
        while(q--)
        {
            scanf("%d%I64d",&u,&k);
            Mat B=Mod_pow(A,k,mod);
            for(int i=0;i<n;i++)
                printf("%I64d ",B.mat[u-1][i]);
            printf("\n");
        }
    }
    return 0;
}