HDU 4686 (推公式+矩陣快速冪)

題目連線：傳送門 

題意： 

給定ai,bi的推倒公式，求sigma(ai*bi) ( 0<= i < n)

分析：

_A0 = A0
A_i = A_i-1*Ax+Ay
B₀ = B0
B_i = B_i-1*Bx+By

設Fi表示Ai*Bi

Fi = Ai * Bi

===>    Fi = (Ai-1*Ax + Ay)*(Bi-1*Bx + By)

===>    Fi = Ax*Bx*(Ai-1*Bi-1) + Ax*By*Ai-1 + Bx*Ay*Bi-1

===>    Fi = Fi-1 + Ax*By*Ai-1 + Bx*Ay*Bi-1 + Ay*By

| F[i] |   | 1 , Ax*By , Bx*Ay , Ay*By |    | F[i-1] |

| Ai   | = | 0 ,   Ax  ,   0   ,   Ay  |  * |  Ai-1  |

| Bi   |   | 0 ,   0   ,   Bx  ,   By  |    |  Bi-1  |

| 1    |   | 0 ,   0   ,   0   ,    1  |    |    1   |

設sum[i] = sigma(ak*bk) ( 0<= k < i)

sum[i] = sum[i-1] + F[i]

sum[i] = sum[i-1] + Fi-1 + Ax*By*Ai-1 + Bx*Ay*Bi-1 + Ay*By

因此可以得到一個包含sum的狀態轉移矩陣

| sum[i] |   | 1 , 1 , Ax*By , Bx*Ay , Ay*By |    | sum[i-1] |

|  F[i]  |   | 0 , 1 , Ax*By , Bx*Ay , Ay*By |    |  F[i-1]  |

|    Ai        | =   | 0 ,     0 ,      Ax   ,         0   ,      Ay    |   *   |      Ai-1       |

|  Bi    |   | 0 , 0 ,   0   ,   Bx  ,   By  |    |   Bi-1   |

|  1     |   | 0 , 0 ,   0   ,   0   ,    1  |    |     1    |

這題需要注意的是n==0 輸出0，注意使用long long,防止中間運算溢位

具體程式碼如下：

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


const int N = 5;

typedef __int64 LL;

const LL mod = 1e9+7;

struct matrix{
    LL a[N][N];
    matrix(){
        memset(a,0,sizeof(a));
    }
}I;

void init(){
    for(int i=0;i<N;i++)
        I.a[i][i]=1;
}

matrix multi(matrix A, matrix B){
    matrix C ;
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            C.a[i][j]=0;
            for(int k=0;k<N;k++){
                C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j]%mod)%mod;
            }
        }
    }
    return C;
}

matrix quick(matrix A, LL n){
    matrix ans = I;
    while(n){
        if(n&1) ans = multi(A,ans);
        n>>=1;
        A=multi(A,A);
    }
    return ans;
}

int main()
{
    init();
    LL n,a0,b0,x1,y1,x2,y2;
    while(~scanf("%I64d",&n)){
        scanf("%I64d%I64d%I64d",&a0,&x1,&y1);
        scanf("%I64d%I64d%I64d",&b0,&x2,&y2);
        if(n==0){
            puts("0");
            continue;
        }
        x2=x2%mod;b0=b0%mod;y2=y2%mod;
        x1=x1%mod;a0=a0%mod;y1=y1%mod;
        matrix A=matrix();
        A.a[0][0]=1LL;
        A.a[0][1]=x1*x2%mod;
        A.a[0][2]=x1*y2%mod;
        A.a[0][3]=x2*y1%mod;
        A.a[0][4]=y1*y2%mod;
        A.a[1][1]=x1*x2%mod;
        A.a[1][2]=x1*y2%mod;
        A.a[1][3]=x2*y1%mod;
        A.a[1][4]=y1*y2%mod;
        A.a[2][2]=x1;
        A.a[2][4]=y1;
        A.a[3][3]=x2;
        A.a[3][4]=y2;
        A.a[4][4]=1LL;
        LL sum = 0;
        matrix ans = quick(A,n-1);
        sum =(sum + a0*b0%mod*ans.a[0][0]%mod)%mod;
        sum =(sum + a0*b0%mod*ans.a[0][1]%mod)%mod;
        sum =(sum + a0*ans.a[0][2]%mod)%mod;
        sum =(sum + b0*ans.a[0][3]%mod)%mod;
        sum =(sum + ans.a[0][4]%mod)%mod;
        printf("%I64d\n",sum);
    }
    return 0;
}