hdu 1588 Gauss Fibonacci 矩阵

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
LL maxn;
LL mod;
/*LL a[maxn][maxn]={
            {1,0,0,0},
            {1,3,1,0},
            {0,2,0,1},
            {0,7,0,0}
};*/
struct matrix{
    LL f[4][4];
};
matrix mul(matrix a,matrix b)
{
    matrix c;
    memset(c.f,0,sizeof(c.f));
    LL i,j,k;

    for(k=0;k<maxn;k++)
    {
        for(i=0;i<maxn;i++)
        {
            if(!a.f[i][k])continue;
            for(j=0;j<maxn;j++)
            {
                //cout<<mod<<endl;
                if(!b.f[k][j])continue;
                c.f[i][j]=(c.f[i][j]+a.f[i][k]*b.f[k][j])%mod;
            }
        }
    }
    return c;
}
matrix pow_mod(matrix a,LL b)
{
    matrix s;
    LL i,j;
    for(i=0;i<maxn;i++)
        for(j=0;j<maxn;j++)
        s.f[i][j]=(i==j?1:0);
    if(b<0)
    {
        s.f[0][0]=0;s.f[0][1]=1;
        s.f[1][0]=1;s.f[1][1]=-1;
        return s;
    }
    while(b)
    {

        if(b&1)
            s=mul(s,a);
        a=mul(a,a);
        b=b>>1;
    }
    return s;
}
int main()
{
    LL k,b,n;
    while(cin>>k>>b>>n>>mod)
    {
        if(k==0&&b==0)
        {
            cout<<0<<endl;
            continue;
        }
        matrix e,g,p,q,h;
        e.f[0][0]=e.f[0][1]=e.f[1][0]=1;
        e.f[1][1]=0;
        maxn=2;
        p=pow_mod(e,k);
        q=pow_mod(e,b-1);
        maxn=4;
        memset(g.f,0,sizeof(g.f));
        g.f[0][0]=g.f[1][1]=g.f[2][0]=g.f[3][1]=1;
        g.f[2][2]=p.f[0][0];g.f[2][3]=p.f[0][1];
        g.f[3][2]=p.f[1][0];g.f[3][3]=p.f[1][1];
        e=pow_mod(g,n);
        h.f[0][0]=e.f[2][0];h.f[0][1]=e.f[2][1];
        h.f[1][0]=e.f[3][0];h.f[1][1]=e.f[3][1];
        maxn=2;
        e=mul(q,h);
        LL ans=e.f[0][0];
        cout<<ans<<endl;
    }
    return 0;
}
/*
    |f[n] f[n-1]|=|f[1] f[0]|*|1 1|^n-1
                              |1 0|  
    令A=|1 1|,B=|1 1|^k;
        |1 0|   |1 0|
    sum=f[b]+f[k+b]+...+f[(n-1)*k+b]
    =|f[1] [0]|*(A^(b-1)+A^(k+b-1)+...+A^((n-1)*k)+b-1)
    =|1 0|*A^(b-1)*(B^0+B^1+...+B^(n-1));
    
    令b[n]=B^n,s[n]=b[0]+b[1]+...b[n-1];
    |s[n] b[n]|=|s[0] b[0]|*|I 0|^n(I 为单位矩阵)
                            |I B|
    求得s[n](这是个矩阵)
    sum=|1 0|*A^(b-1)*s[n];
    
    注意下A^(-1)是A^(-1)*A=I,A^(-1)=|0  1|
                                    |1 -1|
                                    
    这里汗一下，用|f[n] f[n-1]|用习惯了，忘了可以|f[n+1] f[n]|的，这样可以避免
A^(-1)的出现。
*/