HDU 1588Gauss Fibonacci（矩阵）

链接：http://acm.hdu.edu.cn/showproblem.php?pid=1588

题意很好理解，不多说。现在开始做矩阵专题，都是入门类型的，但是对于矩阵的构造还是很难想到，所以这个题目也想了两天。主要想说说我在做这题时构造矩阵的想法。

首先，斐波拉契数列是可以写成矩阵形式的，这应该都懂，所以无形之中给了矩阵构造的提示。下边就是我做完后整理的思路。原谅菜鸟的思路

代码：

#include <string.h>
#include <stdio.h>
//#define mod 2008
#define MAX_dimension 2 
typedef int MATRIX_TYPE;
typedef __int64 MAX_INT_TYPE; //矩阵乘法运算时和加法运算时，为避免溢出，在取模之前用这个类型暂存。
typedef MATRIX_TYPE Matrix[MAX_dimension][MAX_dimension]; //矩阵类型的定义
int ndim=MAX_dimension;

int mod;
Matrix A={ {1,1},
           {1,0}};
void m_zero(Matrix  x)
{
    memset(x, 0, sizeof(MATRIX_TYPE)*MAX_dimension*MAX_dimension);
}

void m_one(Matrix  x)
{
    int i;
    m_zero(x);
    for(i=0;i<ndim;++i)
        x[i][i]=1;
}

void m_copy(Matrix  src,Matrix  dest)
{
    memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_dimension*MAX_dimension);
}
//z=x+y;
void m_add(Matrix  x,Matrix  y,Matrix  z)
{
    int i,j;
    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
            if(mod<=1)z[i][j]=x[i][j]+y[i][j];
            else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;
}
//c=a*b
void m_multiple(Matrix  a,Matrix b,Matrix c)
{
    int i,j,k;
    MAX_INT_TYPE t;

    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
        {
            c[i][j]=0;
            if(mod<=1)
                for(k=0;k<ndim;k++) c[i][j]+=a[i][k]*b[k][j];
            else
                for(k=0;k<ndim;k++)
                {
                    c[i][j]+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod;
                    c[i][j]%=mod;
                }

        }
}
//y=x^n
void m_pow(Matrix x, unsigned int n, Matrix y)
{
    Matrix temp,temp_x;
    m_one(y);
    m_copy(x,temp_x);
    while(n)
    {
        if (n & 1)
        {
            m_multiple(temp_x,y,temp);
            m_copy(temp,y);
        }
        if ((n >>= 1))
        {
            m_multiple(temp_x,temp_x,temp);
            m_copy(temp,temp_x);
        }
    }
}

//根据已经算出来的x^1 x^2 x^4 ...去求矩阵的幂y=x^n
void m_pow_with_saved(Matrix  x_p[],unsigned int n, Matrix y)
{
    Matrix temp;
    m_one(y);
    for(int k=1;n;++k,n>>=1)
    {
        if (n & 1)
        {
            m_multiple(x_p[k],y,temp);
            m_copy(temp,y);
        }
    }
}

//计算根据已经算出来的x^1 x^2 x^4 ...去求矩阵的幂和y=sum(x^n)
void m_pow_sum1(Matrix  x_p[],unsigned int n, Matrix y)
{
    m_zero(y);
    if(n==0)return;
    if(n==1) m_copy(x_p[1],y);
    else
    {
        Matrix temp;
        
        m_pow_sum1(x_p,n>>1,temp);//计算 x^1+...+^(n/2)
        m_add(temp,y,y);
        
        //开始计算 x^(1+n/2)+...+x^n
        Matrix temp2;
        m_pow_with_saved(x_p,n>>1,temp2);
        Matrix temp3;
        m_multiple(temp,temp2,temp3);
        m_add(temp3,y,y);
        if(n&1)
        {
            m_multiple(temp2,temp2,temp3);//计算x^(n-1)
            
            m_multiple(temp3,x_p[1],temp2);//计算 x^n
            
            m_add(temp2,y,y);//add x^n
        }
    }

}
//计算x^1 x^2 x^4 ...，然后调用求幂和的地方
void m_pow_sum(Matrix x, unsigned int n, Matrix y)
{
    //calculate x^1 x^2 x^4 ... x^logn
    unsigned int i=0,logn,n2=n;
    for(logn=0,n2=n;n2;logn++,n2 >>= 1);
    Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)];
    m_one(pow_arry[0]);
    m_copy(x,pow_arry[1]);
    for(i=1;i<logn;i++)
    {
        m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]);
    }

    m_pow_sum1(pow_arry,n,y);
    delete []pow_arry;
}

int main()
{
    //freopen("in.txt","r",stdin);
    int k,b,n;
    Matrix ans,res,sum;

    while(scanf("%d%d%d%d",&k,&b,&n,&mod)!=EOF)
    {
        if(b)
         m_pow(A,b,ans);//计算A^b
        else
            m_one(ans);

        m_pow(A,k,res);//计算A^k

        m_zero(sum);
        m_pow_sum(res,n-1,sum);//计算A^(k*i)之和
        m_multiple(ans,sum,res);//计算A^b*(sum(A^(k*i))),1<=i<=n
        m_add(res,ans,sum);//加上i等于0的值
        printf("%d\n",sum[0][1]%mod);
    }
    return 0;
}