【HDOJ】3509 Buge's Fibonacci Number Problem

快速矩阵幂，系数矩阵由多个二项分布组成。
第1列是（0，（a+b）^k）
第2列是（0，（a+b）^(k-1),0）
第3列是（0，（a+b）^(k-2),0，0）
以此类推。

/* 3509 */
#include <iostream>
#include <string>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <vector>
#include <deque>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cstring>
#include <climits>
#include <cctype>
#include <cassert>
#include <functional>
#include <iterator>
#include <iomanip>
using namespace std;
//#pragma comment(linker,"/STACK:102400000,1024000")

// #define DEBUG
#define sti                set<int>
#define stpii            set<pair<int, int> >
#define mpii            map<int,int>
#define vi                vector<int>
#define pii                pair<int,int>
#define vpii            vector<pair<int,int> >
#define rep(i, a, n)     for (int i=a;i<n;++i)
#define per(i, a, n)     for (int i=n-1;i>=a;--i)
#define clr                clear
#define pb                 push_back
#define mp                 make_pair
#define fir                first
#define sec                second
#define all(x)             (x).begin(),(x).end()
#define SZ(x)             ((int)(x).size())
#define lson            l, mid, rt<<1
#define rson            mid+1, r, rt<<1|1

typedef struct mat_t {
    __int64 m[55][55];
    
    mat_t() {
        memset(m, 0, sizeof(m));
    }
} mat_t;

const int maxn = 55;
__int64 A[maxn], B[maxn], F1[maxn], F2[maxn];
int mod, L;

mat_t matMult(mat_t a, mat_t b) {
    mat_t c;
    
    rep(k, 0, L) {
        rep(i, 0, L) {
            if (a.m[i][k]) {
                rep(j, 0, L) {
                    if (b.m[k][j]) {
                        c.m[i][j] = (c.m[i][j] + a.m[i][k]*b.m[k][j]%mod)%mod;
                    }
                }
            }
        }
    }
    
    return c;
}

mat_t matPow(mat_t a, int n) {
    mat_t ret;
    
    rep(i, 0, L)    ret.m[i][i] = 1;
    
    while (n) {
        if (n & 1)
            ret = matMult(ret, a);
        a = matMult(a, a);
        n >>= 1;
    }
    
    return ret;
}

int main() {
    ios::sync_with_stdio(false);
    #ifndef ONLINE_JUDGE
        freopen("data.in", "r", stdin);
        freopen("data.out", "w", stdout);
    #endif
    
    int t;
    int f1, f2, a, b;
    int n, k_;
    __int64 ans;
    mat_t e, tmp;
    int i, j, k;
    __int64 c;
    int l, r, nc;
    
    scanf("%d", &t);
    while (t--) {
        scanf("%d %d %d %d %d %d %d", &f1,&f2, &a,&b, &k_, &n, &mod);
        L = k_ + 2;
        A[0] = B[0] = F1[0] = F2[0] = 1;
        for (i=1; i<L; ++i) {
            A[i] = A[i-1] * a % mod;
            B[i] = B[i-1] * b % mod;
            F1[i] = F1[i-1] * f1 % mod;
            F2[i] = F2[i-1] * f2 % mod;
        }
        
        memset(e.m, 0, sizeof(e.m));
        e.m[0][0] = e.m[1][0] = 1;
        for (j=1,k=k_; j<L; ++j,--k) {
            for (i=1,c=1,nc=k+1,r=k,l=1; nc; ++i,--nc,c=c*r/l,--r,++l) {
                e.m[i][j] = (c % mod) * A[k+1-i] % mod * B[i-1] % mod;
            }
        }
        #ifdef DEBUG
        for (i=0; i<L; ++i) {
            for (j=0; j<L; ++j)
                printf("%I64d ", e.m[i][j]);
            putchar('
');
        }
        #endif
        
        tmp = matPow(e, n-1);
        
        ans = F1[k_] * tmp.m[0][0] % mod;
        for (j=1; j<L; ++j) {
            ans = (ans + F2[k_+1-j]*F1[j-1]%mod*tmp.m[j][0]%mod)%mod;
        }
        
        printf("%I64d
", ans);
    }
    
    #ifndef ONLINE_JUDGE
        printf("time = %d.
", (int)clock());
    #endif
    
    return 0;
}