HDU 4686 矩阵快速幂

题解

1. 拿到手, 想到是递推转矩阵快速幂， 但是我想的复杂了
   刚开始想的是用a[n] * b[n]的递推构造矩阵+ 分块 + 快速幂， 不过这时候开的矩阵是10 * 10一组最大多O（size^3） * log n
   size为矩阵阶数
2. 最早思路：
   $a_ib_i=a_{i-1}b_{i-1}AxBx+a_{i-1}AxBy+b_{i-1}AyBx+AyBy$
   所以:
   $$\left[ \begin{matrix} AxBx & AxBy & AyBx& By & 0 \\ 0 & Ax & 0 & 1 & 0\\ 0 & 0 & Bx & 0 &1\\ 0 &0&0&1&0\\ 0&0&0&0&1\\ \end{matrix} \right] * \left[ \begin{matrix} a_{i-1}b_{i-1} \\a_{i-1} \\ b_{i-1} \\ Ay \\ By\\ \end{matrix} \right] = \left[ \begin{matrix} a_ib_i \\ a_i\\ b_i\\ Ay\\ By \\ \end{matrix} \right]$$
   $记B=\left[ \begin{matrix} AxBx & AxBy & AyBx& By & 0 \\ 0 & Ax & 0 & 1 & 0\\ 0 & 0 & Bx & 0 &1\\ 0 &0&0&1&0\\ 0&0&0&0&1\\ \end{matrix} \right]\ S_i = \left[ \begin{matrix} a_ib_i \\ a_i\\ b_i\\ Ay\\ By \\ \end{matrix} \right]$
   $B*S_{i-1}=S_i，那么Sum（n） = \sum_{i=0}^{n-1}S_i, AOD(n)是Sum(n)的首项$
3. $Sum(n) = \sum_{i=0}^{n-1}S_i=\sum_{i=0}^{n-1}B^i*S_0=(\sum_{i=0}^{n-1}B^i) * S_0$
   令$T_n = \sum_{i=0}^{n-1}B^i, T_1 = E$
   那么$T_n = E + B*T_{n-1}$
   构造分块矩阵递推：
   $$\left[ \begin{matrix} T_n & E\\ \end{matrix} \right]=\left[ \begin{matrix} T_{n-1} & E\\ \end{matrix} \right]*\left[ \begin{matrix} B & 0\\ E & E \\ \end{matrix} \right]=\left[ \begin{matrix} T_1 & E\\ \end{matrix} \right]*\left[ \begin{matrix} B & 0\\ E & E \\ \end{matrix} \right]^{n-1}$$
   所以$Sum(n) = \left[ \begin{matrix} E & E\\ \end{matrix} \right]*\left[ \begin{matrix} B & 0\\ E & E \\ \end{matrix} \right]^{n-1}*S_0$
   矩阵10*10, T了。。。
4. 然后直接在递推时候考虑 AOD(n)
   即: AOD(n) = AOD(n-1) + a_ib_i,再次递推简化5 * 5, 细节处理A, 还是太弱了啊， 矩阵练了这么多。。

code

#include <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
const ll mod = 1000000007;

ll n, a0, ax, ay, b0, bx, by;

struct  Matrix{
    ll mat[5][5];
    void clr(){
        memset(mat, 0, sizeof mat);
    }
    void E(){
        clr();
        for(int i = 0; i < 5; ++i) mat[i][i] = 1;
    }
    Matrix operator * (const Matrix &b){
        Matrix res;
        res.clr();

        for(int i = 0; i < 5; ++i)
            for(int j = 0; j < 5; ++j)
                for(int k = 0; k < 5; ++k)
                    res.mat[i][j] = (res.mat[i][j] + mat[i][k] * b.mat[k][j]) % mod;
        return res;
    }

    Matrix operator ^ (ll b){
        Matrix res, tmp;
        res.E();
        memcpy(tmp.mat, mat, sizeof mat);

        while(b){
            if(b & 1) res = res * tmp;
            tmp = tmp * tmp;
            b >>= 1;
        }
        return res;
    }
};
int main(){

    while(cin >> n){
        cin >> a0 >> ax >> ay >> b0 >> bx >> by;
        if(n == 0) {
            cout << "0" << endl;
            continue;
        }
        a0 %= mod;
        ax %= mod;
        ay %= mod;
        b0 %= mod;
        bx %= mod;
        by %= mod;
        Matrix init, tmp;
        init.clr();
        tmp.clr();
        init.mat[0][0] = a0 * b0 % mod;
        init.mat[0][1] = a0 * b0 % mod;
        init.mat[0][2] = a0;
        init.mat[0][3] = b0;
        init.mat[0][4] = 1;
        /*tmp底数**/
        tmp.mat[0][0] = 1;
        tmp.mat[1][0] = ax * bx % mod;
        tmp.mat[2][0] = ax * by % mod;
        tmp.mat[3][0] = ay * bx % mod;
        tmp.mat[4][0] = ay * by % mod;
        tmp.mat[1][1] = ax * bx % mod;
        tmp.mat[2][1] = ax * by % mod;
        tmp.mat[3][1] = ay * bx % mod;
        tmp.mat[4][1] = ay * by % mod;
        tmp.mat[2][2] = ax;
        tmp.mat[4][2] = ay;
        tmp.mat[3][3] = bx;
        tmp.mat[4][3] = by;
        tmp.mat[4][4] = 1;

        init = init * (tmp ^ (n - 1));
        cout << init.mat[0][0] << endl;

    }

    return 0;
}

— 待修改