hdu-2157-How many ways?? （矩阵）（dp优化加速）

题目链接：

题意：HDU2157How many ways??
有一个有向图，有N个点，M条边，点的编号从0~N-1，然后有T组询问，每组询问是(A,B,k)，问你从A点出发到B点恰好走k步的方案数。结果对1000取余。
(N<=20,M<=100,A,B<N,k<20,T<=100)

思路：

•有时候dp的递推式子也可以表示成常系数线性递推数列的形式，这时候就可以用矩阵快速幂进行优化了。

一开始用dp的思路来思考，设dp[s][t][k]表示从s点到t点共k步需要的方案数

              s(k)                                               =       sk-1                             *                  ans(行i能否到达列j，能到达就是1）

（dp[1][1][k]  dp[2][1][k] ......dp[n][1][k]）   （dp[1][1][k]  dp[2][1][k] ......dp[n][1][k]）    (0    0  .........1   )

（dp[1][2][k]  dp[2][2][k] .......dp[n][2][k]）   （dp[1][2][k]  dp[2][2][k] .......dp[n][2][k]）   (1    0  .........1   )

（dp[1][3][k] dp[2][3][k] .......dp[n][3][k]）= （dp[1][3][k] dp[2][3][k] .......dp[n][3][k]）* (1    0  .........1   )

   ....

（dp[1][n][k] dp[2][n][k] .......dp[n][n][k]）    （dp[1][n][k] dp[2][n][k] .......dp[n][n][k]）  (1    0  .........0  )

最后总结为

（dp[1][1][k]  dp[2][1][k] ......dp[n][1][k]）  

（dp[1][2][k] dp[2][2][k] .......dp[n][2][k]）  

（dp[1][3][k] dp[2][3][k] .......dp[n][3][k]）=  ans^k;

   ....

（dp[1][n][k] dp[2][n][k] .......dp[n][n][k]）   

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
using namespace std;
#define ll long long
#define N 50
ll mod=1000;
ll n,m;
struct Matrix
{
    ll r,c;
    ll m[N][N];
    Matrix(){}
    Matrix(ll r,ll c):r(r),c(c){}
    Matrix operator *(const Matrix& B)//乘法
    {
        Matrix T(r,B.c);
        for(int i=1;i<=T.r;i++)
        {
            for(int j=1;j<=T.c;j++)
            {
                ll tt = 0;
                for(int k=1;k<=c;k++)
                    tt +=( m[i][k]*B.m[k][j]) % mod;
                T.m[i][j] = tt % mod;
            }
        }
        return T;
    }
	Matrix operator =(const Matrix& B)//复制
    {
        for(int i=1;i<=r;i++)
        {
            for(int j=1;j<=c;j++)
            {
                m[i][j] = B.m[i][j];
            }
        }
    }
    Matrix operator +(const Matrix& B)//加法
    {
        for(int i=1;i<=r;i++)
        {
            for(int j=1;j<=c;j++)
            {
                m[i][j]+= B.m[i][j];
            }
        }
    }
    Matrix Unit(ll h) // 对角线矩阵
    {
        Matrix T(h,h);
        memset(T.m,0,sizeof(T.m));
        for(int i=1;i<=h;i++)
            T.m[i][i] = 1;
        return T;
    }
    Matrix Pow(ll n)  //矩阵幂
    {
        Matrix P = *this,Res = Unit(r);
        while(n!=0)
        {
            if(n&1)
                Res =Res*P;
            P = P*P;
            n >>= 1;
        }
        return Res;
    }
    void Print()//输出
    {
        for(int i=1;i<=r;i++)
        {
            for(int j=1;j<=c;j++)
                printf("%d ",m[i][j]);
            printf("\n");
        }
    }
}Single;
int main(){
	int u,v,w;
	int t;
	while(~scanf("%lld%lld",&n,&m))
	{
		
		if(n==0&&m==0)break;
		Matrix a(n,n),b(n,n);
		for(int i=1;i<=n;i++)
		{
			for(int j=1;j<=n;j++)
			{
				a.m[i][j]=0;
			}
		}
		for(int i=0;i<m;i++)
		{
			scanf("%d%d",&u,&v);
			u++,v++;
			a.m[u][v]=1;
		}

		scanf("%d",&t);
		for(int i=0;i<t;i++)
		{
			scanf("%d%d%d",&u,&v,&w);
			u++,v++;
			b=a;
			b=a.Pow(w);
			printf("%lld\n",b.m[u][v]);
		}
	}
	return 0;
}