计蒜客 嘟嘟数列——构造矩阵

传送门

嘟嘟最近发现了一个数列，将其命名为嘟嘟数列，其中第 n 项的值 dudu ( n ) = $a_0$ × $b_0$ + $a_1$ × $b_1$ + … + $a_{n-1}$ × $b_{n-1}$ .

已知 $a_0$ 和 $b_0$ ，$a_i$ = $a_{i-1}$ × AX + AY ， $b_i$ = $b_{i-1}$ × BX + BY 。

求 dudu ( n ) 模 1000000007 的值。

输入格式

输入包含多组测试数据，对于每组测试数据：

第一行包含包含一个整数 n ( n ≤ $10^{18}$ ) 。

第二行包含三个整数 $a_0$ AX AY ($a_0$, AX, AY ≤ 2 × $10^9$ )。

第三行包含三个整数 $b_0$ BX BY ($b_0$, BX, BY ≤ 2 × $10^9$ )。

输出格式

对于每组测试数据，输出一个整数。

样例输入

2 

1 2 3 

1 2 3 

8 

1 2 3 

5 6 7 

样例输出

26 

994233556 

解题思路：
看到这个题目之后，$n$ 的数据范围太大，暴力肯定超时，所以应该想一 个 $log(n)$ 的复杂度，而且这还是一个函数，所以自然就想到了构造矩阵，然后通过矩阵快速幂来得到答案，最关键的还是如何构造这个矩阵。首先我们令

$$f(n)=a_i*b_i=(a_{i-1}*Ax+Ay)*(b_{i-1}*Bx+By)$$
我们将这个式子展开后得到： $$f(n)=AxBx*f(n-1)+AxBy*a_{i-1}+AyBx*b_{i-1}+AyBy$$
我们令 $f_a(i)=a_{i}, f_b(i)=b_{i}$, 所以得到： $f(n)=AxBx*f(n-1)+AxBy*f_a(i)+AyBx*f_b(i)+AyBy$
由此我们大体上可以想出矩阵会有四项
因为 $f_a(n) = Ax*f_a(n-1)+Ay$ ， $f_b(n) = Bx*f_b(n-1)+By$，所以就会增加两项， $Ay,\ \ By$
又因为我们是要求的 $ans(n)=\sum_{i=0}^{n-1}f(i)=ans(n-1)+f(n-1)$，
所以有 $\{ans(n),f(n-1),f_a(n-1),f_b(n-1),AyBy,Ay,By\}*A=\{ans(n+1),f(n),f_a(n),f_b(n),AyBy,Ay,By\}$
可以构造出矩阵 $A$ 如下：
$$\begin{matrix} 1 & 0 & 0 &0 & 0 &0 &0\\ AxBx & AxBx & 0 &0 & 0 &0 &0\\ AxBy & AxBy & Ax &0 & 0 &0 &0\\ AyBx & AyBx & 0 &Bx & 0 &0 &0\\ 1 & 1 & 0 &0 & 1 &0 &0\\ 0 & 0 & 1 &0 & 0 &1 &0\\ 0 & 0 & 0 &1 & 0 &0 &1\\ \end{matrix}$$

构造出矩阵来之后剩下的就简单了，就是矩阵快速幂了。
代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
using namespace std;
typedef long long LL;
const LL MOD = 1000000007;
const int N = 7;
typedef struct{
    LL mat[N][N];
}Matrix;

Matrix Multi(Matrix a, Matrix b){
    Matrix c;
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            c.mat[i][j] = 0;
            for(int k=0; k<N; k++){
                c.mat[i][j] = (c.mat[i][j] + (a.mat[i][k] * b.mat[k][j]))%MOD;
                ///c.mat[i][j] %= MOD;
            }
        }
    }
    return c;
}

Matrix Pow(Matrix a, LL b){
    Matrix p;
    memset(p.mat, 0, sizeof(p));
    for(int i=0; i<N; i++) p.mat[i][i] = 1;
    while(b){
        if(b & 1)
            p = Multi(p, a);
        b>>=1;
        a = Multi(a, a);
    }
    return p;
}
int main(){
    LL n;
    while(~scanf("%lld",&n)){
        LL a0, ax, ay, b0, bx, by;
        scanf("%lld%lld%lld%lld%lld%lld",&a0, &ax, &ay, &b0, &bx, &by);
        if(n == 0){
            puts("0");
            continue;
        }
        Matrix p;
        memset(p.mat, 0, sizeof(p.mat));
        p.mat[0][0] = 1;
        p.mat[1][0] = ax*bx%MOD, p.mat[1][1] = ax*bx%MOD;
        p.mat[2][0] = ax*by%MOD, p.mat[2][1] = ax*by%MOD, p.mat[2][2] = ax;
        p.mat[3][0] = ay*bx%MOD, p.mat[3][1] = ay*bx%MOD, p.mat[3][3] = bx;
        p.mat[4][0] = 1, p.mat[4][1] = 1, p.mat[4][4] = 1;
        p.mat[5][2] = 1, p.mat[5][5] = 1;
        p.mat[6][3] = 1, p.mat[6][6] = 1;
        p = Pow(p, n-1);
        LL ans = (a0*b0) % MOD;
        ans = (ans + (a0*b0%MOD)*p.mat[1][0] % MOD) %MOD;
        ans = (ans + a0*p.mat[2][0]%MOD) % MOD;
        ans = (ans + b0*p.mat[3][0]%MOD) % MOD;
        ans = (ans + (ay*by%MOD)*p.mat[4][0]%MOD) % MOD;
        ans = (ans + ay*p.mat[5][0]%MOD) % MOD;
        ans = (ans + by*p.mat[6][0]%MOD) % MOD;
        printf("%lld\n",ans);
    }
    return 0;
}