中国（北方）大学生程序设计训练赛（第一周）(Problem E: Water Problem-矩阵快速幂)

已知f(1),f(2),n,
f(n+1)=f(n)+f(n-1)+sin(n*Pi/2),(n>=2)
求f(n)

矩阵快速幂，周期乘4个矩阵

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<functional>
#include<iostream>
#include<cmath>
#include<cctype>
#include<ctime>
using namespace std;
#define For(i,n) for(int i=1;i<=n;i++)
#define Fork(i,k,n) for(int i=k;i<=n;i++)
#define ForkD(i,k,n) for(int i=n;i>=k;i--)
#define Rep(i,n) for(int i=0;i<n;i++)
#define ForD(i,n) for(int i=n;i;i--)
#define RepD(i,n) for(int i=n;i>=0;i--)
#define Forp(x) for(int p=pre[x];p;p=next[p])
#define Forpiter(x) for(int &p=iter[x];p;p=next[p])  
#define Lson (x<<1)
#define Rson ((x<<1)+1)
#define MEM(a) memset(a,0,sizeof(a));
#define MEMI(a) memset(a,127,sizeof(a));
#define MEMi(a) memset(a,128,sizeof(a));
#define INF (2139062143)
#define F (1000000007)
#define eps (1e-3)
#define MAXN (5+10)
#define MAXM (5+10)
typedef long long ll;
ll mul(ll a,ll b){return (a*b)%F;}
ll add(ll a,ll b){return (a+b)%F;}
ll sub(ll a,ll b){return (a-b+llabs(a-b)/F*F+F)%F;}
void upd(ll &a,ll b){a=(a%F+b%F)%F;}

struct M  
{  
    int n,m;  
    ll a[MAXN][MAXN];  
    M(int _n=0){n=m=_n;MEM(a);} 
    M(int _n,int _m){n=_n,m=_m;MEM(a);}
    void mem (int _n=0){n=m=_n;MEM(a);}
    void mem (int _n,int _m){n=_n,m=_m;MEM(a);}
    void fibM(int t) {
        n=m=3; MEM(a)
        a[1][1]=a[1][2]=a[2][1]=a[3][3]=1;;
        a[1][3]=t;
    }
    friend M operator*(M a,M b)  
    {  
        M c(a.n,b.m);  
        For(k,a.m)
            For(i,a.n)  
                For(j,b.m)  
                    c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%F;  
        return c;     
    }  
    friend M operator+(M a,M b)  
    {  
        For(i,a.n)  
            For(j,a.m)  
                a.a[i][j]=(a.a[i][j]+b.a[i][j])%F;  
        return a;
    }  
    void make_I(int _n)  
    {  
        n=m=_n; MEM(a)
        For(i,n) a[i][i]=1;  
    }  
}A,B,C,D,I,ans,E;
M pow2(M a,ll b)  
{  
    M c;c.make_I(a.n);    
    static bool a2[1000000];    
    int n=0;while (b) a2[++n]=b&1,b>>=1;    
    For(i,n)    
    {    
        if (a2[i]) c=c*a;    
        a=a*a;    
    }    
    return c;    
}
int main()
{
    // freopen("E.in","r",stdin);

    A.fibM(1);
    B.fibM(0);
    C.fibM(-1);
    D=((A*B)*C)*B;
    ll f1,f2,n;
    while(scanf("%lld%lld%lld",&f1,&f2,&n)==3) {
        if (n==1) {
            printf("%lld\n",f1);   
        } else if (n==2) {
            printf("%lld\n",f2);   
        } else {
            ll c=n-2;
            ll p=c/4;
            I.make_I(3);
            if (p) I=pow2(D,p);
            c%=4;
            if (c) {
                I=B*I;c--;
            }if (c) {
                I=C*I;c--;
            }if (c) {
                I=B*I;c--;
            }
            ll ans=0;
           upd(ans,mul(I.a[1][1],f2));
           upd(ans,mul(I.a[1][2],f1));
           upd(ans,I.a[1][3]);
           ans=sub(ans,0);
           printf("%lld\n",ans);
        }

    }
    return 0;
}