HDU 4686 Arc of Dream（矩阵快速幂！）

首先，要知道对于任何一个递推式，我们都可以用矩阵法来优化，加快速度求出第n项或前n项和。

这道题开始想的费马小定理，但是发现矩阵不能用，orz。

后来瞄了下题解才发现就是推出矩阵递推式！！没错，就是这么简单暴力！

我们要求的是Sn,变量就有: Sn,Fn,Xn,Yn，由于后面还要加个常数，所以有个1，构建5*5的矩阵

先推出递推式关系：

Xn = (a1*Xn-1 + b1) //自己重新定义了题目的变量~

Yn = (a2*Yn-1 + b2)

Sn-1 = Sn-2 + Fn-1; Fn = Xn*Yn = (a1*Xn-1 + b1)(a2*Yn-1 + b2) = a1*a2*Yn-1*Xn-1 + a2b1Yn-1 + a1b2Xn-1 + b1b2;

之后就可以推出矩阵转移方程了。

Sn-1  1 1 0 0 0 Sn-2

Fn  0 a1*a2  a1b2   a2b1    b1b2 Fn-1

Xn         = 0 0 a1 0 b                                 * Xn-1

Yn 0 0 0 a2 b2 Yn-1

1 0 0 0 0 1 1

当n=0的时候需要特判为0；

#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;

#define LL long long
const LL mod = 1e9+7;

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

struct matrix {
    LL m[20][20];
    matrix () {
        memset(m,0,sizeof(m));
    }
};

matrix mul(matrix a,matrix b) {
    matrix tmp;
    for(int i = 1;i <= 5;i++)
        for(int k = 1;k <= 5;k++)
            if(a.m[i][k])
                for(int j = 1;j <= 5;j++)
                    tmp.m[i][j] = (tmp.m[i][j] + a.m[i][k] * b.m[k][j]) % mod;
    return tmp;
}

matrix powmul(matrix a,LL n) {
    matrix tmp;
    for(int i = 1;i <= 5;i++)
        tmp.m[i][i] = 1;
    while(n) {
        if(n & 1)
            tmp = mul(tmp,a);
        a = mul(a,a);
        n >>= 1;
    }
    return tmp;
}

int main() {
    while(~scanf("%lld%lld%lld%lld%lld%lld%lld",&n,&a0,&ax,&ay,&b0,&bx,&by)) {
        matrix base,ans;
        base.m[1][1] = base.m[1][2] = 1;
        base.m[2][2] = ax*bx % mod; base.m[2][3] = ax*by % mod; base.m[2][4] = bx*ay % mod; base.m[2][5] = ay*by % mod;
        base.m[3][3] = ax % mod; base.m[3][5] = ay % mod;
        base.m[4][4] = bx % mod; base.m[4][5] = by % mod;
        base.m[5][5] = 1;
        /*for(int i = 1;i <= 5;i++)
            for(int j = 1;j <= 5;j++)
                printf("%lld%c",base.m[i][j],j == 5 ? '\n' : ' ');*/
        if(n == 0) {
            printf("0\n");
            continue;
        }
        ans.m[1][1] = a0*b0 % mod;
        ans.m[2][1] = ((a0*ax + ay) % mod)*((b0*bx + by) % mod) % mod;
        ans.m[3][1] = (a0*ax + ay) % mod;
        ans.m[4][1] = (b0*bx + by) % mod;
        ans.m[5][1] = 1;
        /*for(int i = 1;i <= 5;i++)
            for(int j = 1;j <= 5;j++)
                printf("%lld%c",ans.m[i][j],j == 5 ? '\n' : ' ');*/
        base = powmul(base,n-1);
        ans = mul(base,ans);
        printf("%lld\n",ans.m[1][1]);
    }
}