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 long

struct 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;
}