HDU1588___矩阵快速幂and斐波拉契通项的矩阵算法

整个算法来自http://blog.sina.com.cn/s/blog_626631420100vsug.html

他本人写的非常的清楚，感谢。

题目中的公式：
f(0)=0 f(1)=1
f(n)=f(n-1)+f(n-2) (n>=2)
g(i)=k*i+b;
对于Fib序列:
如果用F表示上市中的矩阵就有 F(n+1) = AF(n) 是等比数列
g(i)=k*i+b 是等差数列
所以F(g(i)) = F(b) + F(b+k)+F(b+2k)+....+F(b+nk)
           = F(b) + (A^k)F(b) + (A^2k)F(b)+….+(A^nk)F(b)
提取公因式 F(b)
           = F(b) [ E +A^k + A^2k + ….+ A^nk]   (式中E表示的是单位矩阵)
令 K = A^k 后
                 E +A^k + A^2k + ….+ A^nk 变成 K^0+K^1+K^2+…+K^n
构造矩阵
              

我的代码如下：

#include<cstdio>#include<cstring>#include<iostream>using namespace std;#define LL long longstruct matrix{    LL data[2][2];}E,A,O;struct supMatrix{    matrix m[2][2];}SUPE,SUPA;int k,b,n,M;void init(){    E.data[0][0] = E.data[1][1] = 1;    E.data[0][1] = E.data[1][0] = 0;    A.data[0][0] = A.data[0][1] = A.data[1][0] = 1;    A.data[1][1] = 0;    O.data[0][0] = O.data[0][1] = O.data[1][0] = O.data[1][1] = 0;    SUPE.m[0][0] = SUPE.m[1][1] = E;    SUPE.m[0][1] = SUPE.m[1][0] = O;}matrix m_mul(matrix X,matrix Y){    matrix ans = O;    for(int i=0;i<2;i++)    {        for(int j=0;j<2;j++)        {            for(int k=0;k<2;k++)            {                ans.data[i][j]+=X.data[i][k]*Y.data[k][j];                ans.data[i][j]%=M;            }        }    }    return ans;}matrix m_power(int x,matrix st){    matrix X = st;    matrix Y = E;    if(x==0)        return Y;    while(x!=1)    {        if(x&1)        {            x--;            Y=m_mul(Y,X);        }        else        {            x=x>>1;            X=m_mul(X,X);        }    }    return m_mul(X,Y);}matrix m_add(matrix X,matrix Y){    matrix ret = O;    for(int i=0;i<2;i++)    {        for(int j=0;j<2;j++)        {            ret.data[i][j] = X.data[i][j]+Y.data[i][j];            ret.data[i][j] %= M;        }    }    return ret;}supMatrix sm_mul(supMatrix X,supMatrix Y){    supMatrix ret;    for(int i=0;i<2;i++)    {        for(int j=0;j<2;j++)        {            ret.m[i][j] = O;            for(int k=0;k<2;k++)            {                ret.m[i][j] = m_add(ret.m[i][j],m_mul(X.m[i][k],Y.m[k][j]));            }        }    }    return ret;}supMatrix sm_power(int x,supMatrix st){    supMatrix X = st;    supMatrix Y = SUPE;    if(x==0)        return Y;    while(x!=1)    {        if(x&1)        {            x--;            Y=sm_mul(Y,X);        }        else        {            x=x>>1;            X=sm_mul(X,X);        }    }    return sm_mul(X,Y);}int main(){    init();    while(~scanf("%d%d%d%d",&k,&b,&n,&M))    {        matrix Fb = m_power(b,A);        matrix K = m_power(k,A);        SUPA.m[0][0] = K;        SUPA.m[0][1] = SUPA.m[1][1] = E;        SUPA.m[1][0] = O;        supMatrix ans1 = sm_power(n,SUPA);        matrix KN = ans1.m[0][1];        matrix ans2 = m_mul(Fb,KN);        printf("%I64d\n",ans2.data[1][0]);    }    return 0;}