HLG 1905 f(N) 矩阵快速幂

链接：http://acm.hrbust.edu.cn/index.php?m=ProblemSet&a=showProblem&problem_id=1915

Description：

定义函数f(N) = ∑ai*bi (i从0到N-1)。

并且

a0 = A0

ai = a(i-1)*AX+AY

b0 = B0

bi = b(i-1)*BX+BY

现在求f(N)对1000000007求余的值。

Input:

本题有多组测试数据，对于每组测试数据包含7个正整数，格式如下

N

A0 AX AY

B0 BX BY

N的取值范围不超过10^18，并且其他的数都不超过10^9。

Output:

对于每组测试数据，输出f(N)对1000000007求余的值

代码如下：

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#define MOD 1000000007
#define INF 0x7f
#define RST(N)memset(N, 0, sizeof(N))

typedef unsigned long long ULL;
typedef long long LL;

LL xx[5][5];

void mul1(LL a[5][5], LL b[5][5], LL c[5][5])
{
    for(LL i=0; i<5; i++) {
        for(LL j=0; j<5; j++) {
            c[i][j] = 0;
            for(LL x=0;x<5;x++) {
                c[i][j] += ((a[i][x]%MOD)*(b[x][j]%MOD))%MOD;
                c[i][j] %= MOD;
            }
        }
    }
}

void mul2(LL a[1][5], LL b[5][5], LL c[1][5])
{
    for(LL i=0; i<1; i++) {
        for(LL j=0; j<5; j++) {
            c[i][j] = 0;
            for(LL x=0; x<5; x++) {
                c[i][j] += ((a[i][x]%MOD)*(b[x][j]%MOD))%MOD;
                c[i][j] %= MOD;
            }
        }
    }
}

void Mul(LL x, LL a[5][5])
{
    LL b[5][5];
    if(x == 1) {
        for(LL i=0; i<5; i++) {
            for(LL j=0; j<5; j++) {
                a[i][j] = xx[i][j];
            }
        }
        return ;
    }
    if(x%2 == 0) {
        Mul(x/2, b);
        mul1(b, b, a);
    }else if(x%2 == 1) {
        Mul(x-1, b);
        mul1(b, xx, a);
    }
}

int main()
{
    LL n, Ax, Ay, Bx, By, A0, B0;
    LL start[1][5], mid[5][5], end[1][5];

    while(~scanf("%lld", &n))  {
        RST(xx);
        scanf("%lld %lld %lld", &A0, &Ax, &Ay);
        scanf("%lld %lld %lld", &B0, &Bx, &By);

        A0%=MOD, Ax%=MOD, Ay%=MOD, B0%=MOD, Bx%=MOD, By%=MOD;
        if(n == 1) printf("%lld\n", (A0*B0)%MOD);
        else{
            long long int a1=(A0*Ax%MOD+Ay)%MOD;
            long long int b1=(B0*Bx%MOD+By)%MOD;
            start[0][0] = (A0*B0)%MOD, start[0][1] = (a1*b1)%MOD;
            start[0][2] = a1, start[0][3] = b1, start[0][4] = 1;
            xx[0][0] = xx[1][0] = xx[4][4] = 1;
            xx[1][1] = (Ax*Bx)%MOD, xx[2][1] = (Ax*By)%MOD;
            xx[3][1] = (Ay*Bx)%MOD, xx[4][1] = (Ay*By)%MOD;
            xx[2][2] = Ax, xx[4][2] = Ay;
            xx[3][3] = Bx, xx[4][3] = By;
            Mul(n-1, mid), mul2(start, mid, end);
            printf("%lld\n", end[0][0]%MOD);
        }
    }
    return 0;
}