51nod 1341(推公式、矩阵快速幂)

/*
    51nod 1341
    中文题
    设f(n)=a0*bn+a1*b(n-1)+……+an*b0
    f(n+1)=a0*b(n+1)+a1*bn+……+an*b1+a(n+1)*b0
    观察可得f(n+1)=q*f(n)+a(n+1)*b0
    然后构造矩阵
      {q 3 0}
    A={0 p r}
      {0 0 1}
      {f(n-1)} {f(n)  }
    A*{an    }={a(n+1)}
      {0     } {0     }
        {f(0)} {f(n)  }
    A^n*{a0  }={a(n+1)}
        {0   } {0     }
*/
#include <iostream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define mod 1000000007
#define MAX 4
#define ll long long
#define PI acos(-1)
using namespace std;


struct matrix
{
    ll m[MAX][MAX];
}mat;
matrix matrixmul(matrix a,matrix b)
{
    matrix c;
    for(int i=1; i<=3; i++)
        for(int j=1; j<=3; j++)
        {
            c.m[i][j]=0;
            for(int k=1; k<=3; k++)
                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;
            c.m[i][j]%=mod;
        }
    return c;
}
matrix quickpow(matrix m,int n)
{
    matrix b;
    memset(b.m,0,sizeof(b.m));
    for(int i=1;i<=3;i++)
        b.m[i][i]=1;
    while(n>=1)
    {
        if(n&1)
            b=matrixmul(b,m);
        n=n>>1;
        m=matrixmul(m,m);
    }
    return b;
}
int main()
{
    int p,q,r,n;
    cin>>p>>q>>r>>n;
    mat.m[1][1]=q;mat.m[1][2]=3;mat.m[1][3]=0;//构造A矩阵
    mat.m[2][1]=0;mat.m[2][2]=p;mat.m[2][3]=r;
    mat.m[3][1]=0;mat.m[3][2]=0;mat.m[3][3]=1;
    matrix c;
    c=quickpow(mat,n);
    int ans=((c.m[1][2]*r%mod+mod)%mod+c.m[1][3])%mod;
    cout<<ans<<endl;
    return 0;
}