矩阵高速幂模板篇

转载请注明出处：http://blog.csdn.net/u012860063

计算f[n] = f[n-1] + 2 * f[n-2] + c；

输入：n和mod和c。

代码例如以下：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
struct Matrix
{
    LL m[4][4];
} I,A,B,T;
LL mod, n, c;
int ssize = 3;

Matrix Mul(Matrix a,Matrix b)
{
    int i, j, k;
    Matrix c;
    for(i = 1; i <= ssize; i++)
    {
        for(j = 1; j <= ssize; j++)
        {
            c.m[i][j]=0;
            for(k = 1; k <= ssize; k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j]);
                c.m[i][j]%=mod;
            }
        }
    }
    return c;
}

Matrix quickpagow(LL n)
{
    Matrix m = A, b = I;
    while(n > 0)
    {
        if(n & 1)
            b = Mul(b,m);
        n = n >> 1;
        m = Mul(m,m);
    }
    return b;
}
int main()
{
    while(~scanf("%lld%lld%lld",&n,&mod,&c))
    {
        memset(I.m,0,sizeof(I.m));
        memset(A.m,0,sizeof(A.m));
        memset(B.m,0,sizeof(B.m));
        for(int i = 1; i <= ssize; i++)
        {
            //单位矩阵
            I.m[i][i]=1;
        }
        A.m[1][1] = 1;//初始化等比矩阵
        A.m[1][2] = 2;
        A.m[1][3] = A.m[2][1] = A.m[3][3] = 1;

        B.m[1][1] = 2;
        B.m[2][1] = 1;
        B.m[3][1] = c;
        if(n==1)
        {
            printf("%I64d\n",1%mod);
            continue;
        }
        if(n==2)
        {
            printf("%I64d\n",2%mod);
            continue;
        }
        T = quickpagow(n-2);//注意n-1为负的情况
        T = Mul(T,B);
        printf("%lld\n",T.m[1][1]%mod);
    }
    return 0;
}
/*
1 99 1
2 99 1
3 99 1
4 99 1
5 99 1
6 99 1
*/