【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）
以此类推。
1 /* 3509 */ 2 #include <iostream> 3 #include <string> 4 #include <map> 5 #include <queue> 6 #include <set> 7 #include <stack> 8 #include <vector> 9 #include <deque> 10 #include <algorithm> 11 #include <cstdio> 12 #include <cmath> 13 #include <ctime> 14 #include <cstring> 15 #include <climits> 16 #include <cctype> 17 #include <cassert> 18 #include <functional> 19 #include <iterator> 20 #include <iomanip> 21 using namespace std; 22 //#pragma comment(linker,"/STACK:102400000,1024000") 23 24 // #define DEBUG 25 #define sti set<int> 26 #define stpii set<pair<int, int> > 27 #define mpii map<int,int> 28 #define vi vector<int> 29 #define pii pair<int,int> 30 #define vpii vector<pair<int,int> > 31 #define rep(i, a, n) for (int i=a;i<n;++i) 32 #define per(i, a, n) for (int i=n-1;i>=a;--i) 33 #define clr clear 34 #define pb push_back 35 #define mp make_pair 36 #define fir first 37 #define sec second 38 #define all(x) (x).begin(),(x).end() 39 #define SZ(x) ((int)(x).size()) 40 #define lson l, mid, rt<<1 41 #define rson mid+1, r, rt<<1|1 42 43 typedef struct mat_t { 44 __int64 m[55][55]; 45 46 mat_t() { 47 memset(m, 0, sizeof(m)); 48 } 49 } mat_t; 50 51 const int maxn = 55; 52 __int64 A[maxn], B[maxn], F1[maxn], F2[maxn]; 53 int mod, L; 54 55 mat_t matMult(mat_t a, mat_t b) { 56 mat_t c; 57 58 rep(k, 0, L) { 59 rep(i, 0, L) { 60 if (a.m[i][k]) { 61 rep(j, 0, L) { 62 if (b.m[k][j]) { 63 c.m[i][j] = (c.m[i][j] + a.m[i][k]*b.m[k][j]%mod)%mod; 64 } 65 } 66 } 67 } 68 } 69 70 return c; 71 } 72 73 mat_t matPow(mat_t a, int n) { 74 mat_t ret; 75 76 rep(i, 0, L) ret.m[i][i] = 1; 77 78 while (n) { 79 if (n & 1) 80 ret = matMult(ret, a); 81 a = matMult(a, a); 82 n >>= 1; 83 } 84 85 return ret; 86 } 87 88 int main() { 89 ios::sync_with_stdio(false); 90 #ifndef ONLINE_JUDGE 91 freopen("data.in", "r", stdin); 92 freopen("data.out", "w", stdout); 93 #endif 94 95 int t; 96 int f1, f2, a, b; 97 int n, k_; 98 __int64 ans; 99 mat_t e, tmp; 100 int i, j, k; 101 __int64 c; 102 int l, r, nc; 103 104 scanf("%d", &t); 105 while (t--) { 106 scanf("%d %d %d %d %d %d %d", &f1,&f2, &a,&b, &k_, &n, &mod); 107 L = k_ + 2; 108 A[0] = B[0] = F1[0] = F2[0] = 1; 109 for (i=1; i<L; ++i) { 110 A[i] = A[i-1] * a % mod; 111 B[i] = B[i-1] * b % mod; 112 F1[i] = F1[i-1] * f1 % mod; 113 F2[i] = F2[i-1] * f2 % mod; 114 } 115 116 memset(e.m, 0, sizeof(e.m)); 117 e.m[0][0] = e.m[1][0] = 1; 118 for (j=1,k=k_; j<L; ++j,--k) { 119 for (i=1,c=1,nc=k+1,r=k,l=1; nc; ++i,--nc,c=c*r/l,--r,++l) { 120 e.m[i][j] = (c % mod) * A[k+1-i] % mod * B[i-1] % mod; 121 } 122 } 123 #ifdef DEBUG 124 for (i=0; i<L; ++i) { 125 for (j=0; j<L; ++j) 126 printf("%I64d ", e.m[i][j]); 127 putchar('\n'); 128 } 129 #endif 130 131 tmp = matPow(e, n-1); 132 133 ans = F1[k_] * tmp.m[0][0] % mod; 134 for (j=1; j<L; ++j) { 135 ans = (ans + F2[k_+1-j]*F1[j-1]%mod*tmp.m[j][0]%mod)%mod; 136 } 137 138 printf("%I64d\n", ans); 139 } 140 141 #ifndef ONLINE_JUDGE 142 printf("time = %d.\n", (int)clock()); 143 #endif 144 145 return 0; 146 }
