HDU - 4686 Arc of Dream（5*5构造矩阵快速幂）

An Arc of Dream is a curve defined by following function:

where
a 0 = A0
a i = a i-1*AX+AY
b 0 = B0
b i = b i-1*BX+BY
What is the value of AoD(N) modulo 1,000,000,007?

Input

There are multiple test cases. Process to the End of File.
Each test case contains 7 nonnegative integers as follows:
N
A0 AX AY
B0 BX BY
N is no more than 10 18, and all the other integers are no more than 2×10 9.

Output

For each test case, output AoD(N) modulo 1,000,000,007.

Sample Input

1
1 2 3
4 5 6
2
1 2 3
4 5 6
3
1 2 3
4 5 6

Sample Output

4
134
1902

 构造5*5矩阵同时解决和和项的问题。

【   axbx    0       0        0     1

      axby    ax     0        0     0

     aybx    0       bx      0     0

       ayay    ay     by      1     0

       0          0        0        0     1   】

#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int SIZE=20;
const ll MOD=1000000007;
struct Matrix {
    int r, c;
    ll mat[SIZE][SIZE];
    Matrix() {}
    Matrix(int _r, int _c) {
        r = _r;
        c = _c;
        memset(mat, 0, sizeof(mat));
    }
};
Matrix operator + (Matrix a, Matrix b) {
    Matrix s(a.r, a.c);
    for(int i = 0; i < a.r; i++) {
        for(int j = 0; j < a.c; j++) {
            s.mat[i][j] = (a.mat[i][j] + b.mat[i][j]) % MOD;
        }
    }
    return s;
}
Matrix operator * (Matrix a, Matrix b) {
    Matrix s(a.r, b.c);
    for(int i = 0; i < a.r; i++) {
        for(int k = 0; k < a.c; k++) { //j, k交换顺序运行更快
            for(int j = 0; j < b.c; j++) {
                s.mat[i][j] = (s.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % MOD;
            }
        }
    }
    return s;
}
Matrix pow_mod(Matrix a, ll n) {
    Matrix ret(a.r, a.c);
    for(int i = 0; i < a.r; i++) {
        ret.mat[i][i]=1; //单位矩阵
    }
    Matrix tmp(a);
    while(n) {
        if(n&1) ret = ret * tmp;
        tmp = tmp * tmp;
        n >>= 1;
    }
    return ret;
}
int main()
{
	cin.tie(0);
	ios::sync_with_stdio(false);
	int n;
	while(cin>>n)
	{
		int A0,AX,AY,B0,BX,BY;
		cin>>A0>>AX>>AY>>B0>>BX>>BY;
		if(n == 0) 
		{
            cout << "0" << endl;
            continue;
        }
		Matrix A(5,5);
		A.mat[0][0] = A.mat[0][4] = AX * BX % MOD;
        A.mat[1][0] = A.mat[1][4] = AX * BY % MOD;
        A.mat[1][1] = AX;
        A.mat[2][0] = A.mat[2][4] = BX * AY % MOD;
        A.mat[2][2] = BX;
        A.mat[3][0] = A.mat[3][4] = AY * BY % MOD;
        A.mat[3][3] = 1;
        A.mat[3][1] = AY;
        A.mat[3][2] = BY;
        A.mat[4][4] = 1;
        Matrix ans=(pow_mod(A, n-1));
        Matrix F(1,5);
        F.mat[0][0] = F.mat[0][4] = A0 * B0 % MOD;
        F.mat[0][1] = A0 % MOD;
        F.mat[0][2] = B0 % MOD;
        F.mat[0][3] = 1;
        F = F * ans;	
		cout<<F.mat[0][4]<<endl;
	}
	return 0;
}