HDU 5607 矩阵快速幂

HDU 5607
题目链接：
http://acm.hdu.edu.cn/showproblem.php?pid=5607
题意：
给一个n个点有向图(<=50)，有m条边(<=1000).
现在已知对于一个点，下次可以等可能的走到有出边的下一个点。
从一个点出发，求最后走到某个点的概率x/y.输出x * y^(1e9+5)
思路：
妥妥的用矩阵快速幂优化一个转移矩阵没什么好说的。
自己写的时候刚开始写了一发double爆了精度，然后分数类爆了时间。
最后看别人题解发现x * y^(1e9+5)是可以看成乘上y的逆元，然后这题就解决了。
源码：
写残版（但是AC了）

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <vector>
#include <utility>
using namespace std;
#define mod (1000000007)
#define LL long long
const int MAXN = 50 + 5;
LL gcd(LL a, LL b){return b == 0 ? a : gcd(b, a % b);}
struct Matrix
{
    LL a[MAXN][MAXN];
    int num;
    void init(int _num)
    {
        num = _num;
        for(int i = 1 ; i <= num ; i++){
            for(int j = 1 ; j <= num ; j++) a[i][j] = 0;
            a[i][i] = 1;
        }
    }
    void print()
    {
        for(int i = 1 ; i <= num ; i++){
            for(int j = 1 ; j <= num ; j++) printf("%10I64d ", a[i][j]);
            printf("\n");
        }
    }
    Matrix operator * ( Matrix &rbs)const{
        Matrix ans;
        ans.init(num);
        for(int i = 1 ; i <= num ; i++) ans.a[i][i] = 0;
//        printf("\nans\n");
//        ans.print();
//        printf("\nnow\n");
//        for(int i = 1 ; i <= num ; i++){
//            for(int j = 1 ; j <= num ; j++) printf("%10I64d ", a[i][j]);
//            printf("\n");
//        }
//        printf("\nrbs\n");
//        rbs.print();
        for(int i = 1 ; i <= num ; i++){
            for(int k = 1 ; k <= num ; k++){
                for(int j = 1 ; j <= num ; j++){
                    ans.a[i][j] = (ans.a[i][j] + a[i][k] * rbs.a[k][j] % mod) % mod;
//                    if(i == 1 && j == 1){
//                        printf("i = %d, j = %d, k = %d, a[i][k] = %I64d, a[j][k] = %I64d\n", i, j, k, a[i][k], rbs.a[k][j]);
//                    }
                }
            }
        }
//        printf("\nans\n");
//        ans.print();
//        printf("\n");
        return ans;
    }
}org, cur;
Matrix ppow(Matrix a, int v)
{
    Matrix ans;
    ans.init(a.num);
    while(v){
        if(v & 1)   ans = ans * a;
        a = a * a;
        v >>= 1;
    }
    return ans;
}
LL ppow(LL a, int v)
{
    LL ans = 1;
    while(v){
        if(v & 1)   ans = (ans * a) % mod;
        a = (a * a) % mod;
        v >>= 1;
    }
    return ans;
}
LL rev(LL a){return ppow(a, mod - 2);}
vector<int>lin[MAXN];
int n;
int outdegree[MAXN];
int main()
{
//    freopen("HDU 5607.in", "r", stdin);
    int m;
    while(scanf("%d%d", &n, &m) != EOF){
        for(int i = 1 ; i <= n ; i++)   lin[i].clear(), outdegree[i] = 0;
        int u, v;
        for(int i = 0 ; i < m ; i++){
            scanf("%d%d", &u, &v);
            outdegree[u]++;
            lin[v].push_back(u);
        }

        ///init_matrix
        org.num = n;
        for(int i = 1 ; i <= n ; i++)   for(int j = 1 ; j <= n ; j++) org.a[i][j] = 0;
        for(int i = 1 ; i <= n ; i++){
            for(int j = 0 ; j < (int)lin[i].size() ; j++){
                org.a[i][lin[i][j]] = rev(outdegree[lin[i][j]]);
            }
        }
        int q;
        scanf("%d", &q);
        while(q--){
            scanf("%d%d", &u, &v);
//            printf("\norg\n");
//            org.print();
            cur = ppow(org, v);
//            printf("\ncur\n");
//            cur.print();
            for(int i = 1 ; i <= n ; i++)   printf("%I64d ", cur.a[i][u]);
            printf("\n");
        }
    }
    return 0;
}

分数类版（但是TLE）
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <vector>
#include <utility>
using namespace std;
#define mod (1000000007)
#define LL long long
const int MAXN = 50 + 5;
LL gcd(LL a, LL b){return b == 0 ? a : gcd(b, a % b);}
struct Fraction
{
    LL a, b;
    Fraction(){}
    Fraction(LL _a, LL _b){a = _a, b = _b;}
    Fraction operator + (const Fraction &rbs)const{
        Fraction ans;
//        printf("%I64d %I64d, %I64d %I64d\n", a, b, rbs.a, rbs.b);
        ans.b = rbs.b * b / gcd(rbs.b, b);
        ans.a = a * (ans.b / b) + rbs.a * (ans.b / rbs.b);
        ans.a %= mod;
        ans.b %= mod;
        return ans;
    }
    Fraction operator * (const Fraction &rbs)const{
        Fraction ans;
        ans.a = a * rbs.a;
        ans.b = b * rbs.b;
        LL temp = gcd(ans.a, ans.b);
        if(temp){
            ans.a = ans.a / temp;
            ans.b = ans.b / temp;
            ans.a %= mod;
            ans.b %= mod;
        }
//        printf("in * %I64d %I64d, %I64d %I64d, %I64d %I64d\n", a, b, rbs.a, rbs.b, ans.a, ans.b);
        return ans;
    }
};
struct Matrix
{
    Fraction a[MAXN][MAXN];
    int num;
    void init(int _num)
    {
        num = _num;
        for(int i = 1 ; i <= num ; i++){
            for(int j = 1 ; j <= num ; j++) a[i][j] = Fraction(0, 1);
            a[i][i] = Fraction(1, 1);
        }
    }
    Matrix operator * (const Matrix &rbs)const{
        Matrix ans;
        ans.num = rbs.num;
        for(int i = 1 ; i <= num ; i++) for(int j = 1 ; j <= num ; j++) ans.a[i][j] = Fraction(0, 1);
        for(int i = 1 ; i <= num ; i++){
            for(int k = 1 ; k <= num ; k++){
                for(int j = 1 ; j <= num ; j++){
                    ans.a[i][j] = (ans.a[i][j] + a[i][k] * rbs.a[k][j]);
//                    printf("i = %d, j = %d, k = %d, a[i][k] = %I64d %I64d, rbs.a[k][j] = %I64d %I64d, ans.a[i][j] = %I64d %I64d\n", i, j, k, a[i][k].a, a[i][k].b, rbs.a[k][j].a, rbs.a[k][j].b, ans.a[i][j].a, ans.a[i][j].b);
                }
//                printf("ans.a[%d][%d] = %I64d %I64d\n", i, j, ans.a[i][j].a, ans.a[i][j].b);
            }
        }
        return ans;
    }
}org, cur;
Matrix ppow(Matrix a, int v)
{
//    printf("a.num = %d\n", a.num);
    Matrix ans;
    ans.init(a.num);
    while(v){
        if(v & 1)   ans = ans * a;
        a = a * a;
        v >>= 1;
    }
//    printf("a.num = %d, ans.num = %d\n", a.num, ans.num);
    return ans;
}
LL ppow(LL a, int v)
{
    LL ans = 1;
    while(v){
        if(v & 1)   ans = (ans * a) % mod;
        a = (a * a) % mod;
        v >>= 1;
    }
    return ans;
}
vector<int>lin[MAXN];
int n;
int outdegree[MAXN];
LL ansx[MAXN], ansy[MAXN];
void print(Matrix u)
{
    for(int i = 1 ; i <= u.num ; i++){
        for(int j = 1 ; j <= u.num ; j++)   printf("%I64d/%I64d ", u.a[i][j].a, u.a[i][j].b);
        printf("\n");
    }
}
int main()
{
//    freopen("HDU 5607.in", "r", stdin);
    int m;
    while(scanf("%d%d", &n, &m) != EOF){
        for(int i = 1 ; i <= n ; i++)   lin[i].clear(), outdegree[i] = 0;
        int u, v;
        for(int i = 0 ; i < m ; i++){
            scanf("%d%d", &u, &v);
            outdegree[u]++;
            lin[v].push_back(u);
        }

        ///init_matrix
        org.num = n;
        for(int i = 1 ; i <= n ; i++)   for(int j = 1 ; j <= n ; j++) org.a[i][j] = Fraction(0, 1);
        for(int i = 1 ; i <= n ; i++){
            for(int j = 0 ; j < (int)lin[i].size() ; j++){
                org.a[i][lin[i][j]] = Fraction(1, outdegree[lin[i][j]]);
            }
        }
//        printf("first\n");
        int q;
        scanf("%d", &q);
        while(q--){
            scanf("%d%d", &u, &v);
//            Matrix torg = org;
            cur = ppow(org, v);
//            printf("first\n");
//            print(cur);
//            printf("\n");
//            printf("second\n");
//            print(org);
//            org = torg;
            for(int i = 1 ; i <= n ; i++){
                printf("%I64d", cur.a[i][u].a * ppow(cur.a[i][u].b, 1000000005) % mod);
//                printf(" cur.a[i][u] = %I64d %I64d\n", cur.a[i][u].a, cur.a[i][u].b);
//                if(i == n)  printf("\n");
//                else    printf(" ");
                printf(" ");
            }
            printf("\n");
        }
    }
    return 0;
}