快速矩阵幂,系数矩阵由多个二项分布组成。
第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[][]; mat_t() {
memset(m, , sizeof(m));
}
} mat_t; const int maxn = ;
__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, , L) {
rep(i, , L) {
if (a.m[i][k]) {
rep(j, , 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, , L) ret.m[i][i] = ; while (n) {
if (n & )
ret = matMult(ret, a);
a = matMult(a, a);
n >>= ;
} 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_ + ;
A[] = B[] = F1[] = F2[] = ;
for (i=; i<L; ++i) {
A[i] = A[i-] * a % mod;
B[i] = B[i-] * b % mod;
F1[i] = F1[i-] * f1 % mod;
F2[i] = F2[i-] * f2 % mod;
} memset(e.m, , sizeof(e.m));
e.m[][] = e.m[][] = ;
for (j=,k=k_; j<L; ++j,--k) {
for (i=,c=,nc=k+,r=k,l=; nc; ++i,--nc,c=c*r/l,--r,++l) {
e.m[i][j] = (c % mod) * A[k+-i] % mod * B[i-] % mod;
}
}
#ifdef DEBUG
for (i=; i<L; ++i) {
for (j=; j<L; ++j)
printf("%I64d ", e.m[i][j]);
putchar('\n');
}
#endif tmp = matPow(e, n-); ans = F1[k_] * tmp.m[][] % mod;
for (j=; j<L; ++j) {
ans = (ans + F2[k_+-j]*F1[j-]%mod*tmp.m[j][]%mod)%mod;
} printf("%I64d\n", ans);
} #ifndef ONLINE_JUDGE
printf("time = %d.\n", (int)clock());
#endif return ;
}
05-08 14:49