您的位置：首页 > 其它

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

2013-04-15 15:26 357 查看
整个算法来自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;
}
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 点击这里给我发消息
标签： 
相关文章推荐
新的分享
章节导航