矩阵快速幂

矩阵乘法：

　　　a（n, m）, b(m, p) 为两个二维矩阵 相乘可得矩阵c(n, p)  a中每一行和b中每一列相对应数的乘积之和；

1、

　　hdoj 1575--Tr A   http://acm.hdu.edu.cn/showproblem.php?pid=1575　　

　　A为一个矩阵 求A^K中主对角线之和 ， 结果取模；

#define mod 9973
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
struct matrix
{
    int a[15][15];
}; int n;
matrix mul(matrix a, matrix b)
{
    matrix ret;
    memset(ret.a, 0, sizeof(ret.a));
    for(int i = 0; i < n; i++)
    {
        for(int k = 0; k < n; k++)
            if(a.a[i][k])
                for(int j = 0; j < n; j++)
                    if(b.a[k][j])
                        ret.a[i][j]=(a.a[i][k]*b.a[k][j]+ret.a[i][j])%mod;
    }
    return ret;
}
matrix mpower(matrix tem, int m)
{
    matrix I;
    for(int i = 0; i < 15; i++)
        for(int j = 0; j < 15; j++)
            I.a[i][j] = (i==j);
    while(m)
    {
        if(m&1)
            I = mul(I, tem);
            m >>= 1;
        tem = mul(tem, tem);
    }
    return I;
 }
int main()
{
    int T;
    scanf("%d", &T);
    while(T--)
    {
        int m;
        matrix num;
        scanf("%d%d", &n, &m);
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                scanf("%d", &num.a[i][j]);
        num = mpower(num, m);
        int sum = 0;
        for(int i = 0; i < n; i++)
        {
            sum = (sum+num.a[i][i])%mod;
        }
        printf("%d
", sum);
    }
    return 0;
}

2、

　　Poj 3070--Fibonacci      http://poj.org/problem?id=3070

 　  奇妙, 二维矩阵a{0, 1, 1, 1}^n = F(n) && (n != 0); Fib数列矩阵求法 ； {a,  b} * {0, 1， 1， 1}  = {b, a+b};  

#include <cmath>
#include <cstdio>
#include <cstring>
const int MOD = 1e4;
struct matrix
{
    int a[2][2];
};
matrix Matrixmul(matrix a, matrix b)
{
   // int cnt = 0;
    matrix ret;
    memset(ret.a, 0, sizeof(ret.a));
    for(int i = 0; i < 2; i++)
        for(int k = 0; k < 2; k++)
        {
                if(a.a[i][k])
                for(int j = 0; j < 2; j++)
                {
                    if(b.a[k][j])
                        ret.a[i][j]=(a.a[i][k]*b.a[k][j]+ret.a[i][j])%MOD;
                }
             /*   cnt++;
                for(int L = 0; L < 2; L++){
                    for(int M = 0; M < 2; M++)
                        printf("%d ", ret.a[L][M]);
                    printf("
");
                }

                printf("%d


", cnt);*/
        }
        return ret;
}
matrix Matrixpow(matrix a, int n)
{
    matrix I;
    for(int i = 0; i < 2; i++)
        for(int j = 0; j < 2; j++)
            I.a[i][j] = (i==j);
    while(n)
    {
        if(n&1)
            I = Matrixmul(I, a);
        n >>= 1;
        a = Matrixmul(a, a);
    }
    return I;
}
int main()
{
    int n;
    while(scanf("%d", &n)!= EOF, n != -1)
    {
        if(n == 0)
            printf("0
");
        else
        {
            matrix num, tem;
            num.a[0][0] = 0;
            num.a[0][1] = 1;
            num.a[1][0] = 1;
            num.a[1][1] = 1;
            tem = Matrixpow(num, n);
            printf("%d
", tem.a[0][1]);
        }
    }
    return 0;
}

3、 

　　Poj3233--Matrix Power Series (递推+矩阵) re++!!

　　S = A + A² + A³ + … + A^{k  ；}A为矩阵；

　　例子：

　　　　（1）k = 6 有： S(6) = (1 + A^3) * (A + A^2 + A^3) = (1 + A^3) * S(3)。
　　　　（2）k = 7 有： S(7) = A + (A + A^4) * (A + A^2 + A^3) = A + (A + A^4) * S(3)。

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
int n, m, k; int MOD;
struct mat
{
    int  a[41][41];
};
mat d;
struct
{
    mat expo(mat a, int k)
    {
        /*if(k == 1)
            return a;*/
        mat e;
        memset(e.a, 0, sizeof(e.a));
        for(int i = 0; i < n; i++)
            e.a[i][i] = 1;
        if(k == 0)
            return e;
        while(k)
        {
            if(k&1)
                e=mul(e, a);
            k>>=1;
            a=mul(a, a);
        }
        return e;
    }
    mat mul(mat a, mat b)
    {
        mat ret;
        memset(ret.a, 0, sizeof(ret.a));
        for(int i = 0; i < n; i++)
            for(int k = 0; k < n; k++)
                if(a.a[i][k])
                    for(int j = 0; j < n; j++)
                        if(b.a[k][j])
                        {
                            ret.a[i][j]=a.a[i][k]*b.a[k][j]+ret.a[i][j];
                            if(ret.a[i][j]>=MOD)
                                ret.a[i][j]%=MOD;
                        }
        return ret;
    }
}MmM;//Matrix_muti_Modern;
void print(mat a)
{
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
            printf(j==n-1?"%d
":"%d ", a.a[i][j]);
    }
}
mat add(mat a, mat b)
{
    mat t;
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
        {
            t.a[i][j]=a.a[i][j]+b.a[i][j];
            if(t.a[i][j] >= m)
                t.a[i][j]%=MOD;
        }
    return t;
}
//（1）k = 6 有： S(6) = (1 + A^3) * (A + A^2 + A^3) = (1 + A^3) * S(3)。
//（2）k = 7 有： S(7) = A + (A + A^4) * (A + A^2 + A^3) = A + (A + A^4) * S(3)。
mat sum(int k)
{
    if(k == 1)
        return d;
    if(k&1)
        return add(sum(k-1), MmM.expo(d, k));
    else
    {
        mat s=sum(k>>1);
        return add(s, MmM.mul(s, MmM.expo(d, k>>1))); 
        //S(n>>1)+S(n>>1)*A^(n>>1)=S(n>>1)(1+A^(n>>1)&&S(n>>1)=(A1+...+A(n>>1)); 
    }
}
int main()
{
    mat tem;
    while(scanf("%d%d%d", &n, &k, &m) != EOF)
    {
        MOD = m;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
            {
                scanf("%d", &d.a[i][j]);
                if(d.a[i][j] >= MOD)
                    d.a[i][j]%=MOD;
            }
        tem = sum(k);
        print(tem);
    }
    return 0;
}

4、

　　Poj3735--Training little cats

　　n, m, k; Ai(i =1,2, 3)代表操作； n只猫， k条操作， 重复m次；  主要是构造矩阵；

　　A1 a (a食物数+1)；  A2 a, b(a食物给b); A3 a(吃完所有)；

　　构造矩阵Mat(0--n, 0--n);  A1时Mat[0][a]+1; A2时交换a, b列; A3时a列置0；

　　题解：http://www.cppblog.com/y346491470/articles/157284.html　

#include <cstdio>
#include <cstring>
#include <iostream>
#define LL long long
using namespace std;
struct mat
{
    LL a[101][101];
};
int n, m, k;
mat matinit(mat a)
{
    memset(a.a, 0, sizeof(a.a));
    for(int i = 0; i <= n; i++)
        a.a[i][i] = 1;
    return a;
}
mat matmul(mat a, mat b)
{
    mat res;
    memset(res.a, 0, sizeof(res.a));
    for(int i = 0; i <= n; i++)
        for(int k = 0; k <= n; k++)
            if(a.a[i][k])
                for(int j = 0; j <= n; j++)
                    if(b.a[k][j])
                        res.a[i][j] += a.a[i][k]*b.a[k][j];
    return res;
}
mat matpow(mat a, int k)
{
    mat I;
    memset(I.a, 0, sizeof(I.a));
    for(int i = 0; i <= n; i++)
        I.a[i][i] = 1;
   /* if(k == 0)
        return I; */
    while(k)
    {
        if(k&1)
            I=matmul(I, a);
        k>>=1;
        a=matmul(a, a);
    }
    return I;
}
int main()
{
    while(scanf("%d%d%d", &n, &m, &k) != EOF ,n ,m, k)
    {
        mat ret, num;
        char str[2];  int a, b;
        num = matinit(num);
        /*for(int i = 0; i <= n; i++)
            for(int j = 0; j <= n; j++)
                printf("%d
", num.a[i][j]);*/
        for(int i = 0; i < k; i++)
        {
            cin >> str;
            if(str[0] == 'g')
            {
                scanf("%d", &a);
                num.a[0][a] +=1;
            }
            if(str[0] == 'e')
            {
                scanf("%d", &a);
                for(int i = 0; i <= n; i++)
                    num.a[i][a] = 0;
            }
            if(str[0] == 's')
            {
                scanf("%d%d", &a, &b);
                for(int i = 0; i <= n; i++)
                {
                    LL temp = num.a[i][a];
                    num.a[i][a] = num.a[i][b];
                    num.a[i][b] = temp;
                }
            }
        }
        ret = matpow(num, m);
        for(int i = 1; i <= n; i++)
            printf(i==n?"%lld
":"%lld ", ret.a[0][i]);
    }
    return 0;
}

5、

　　Poj 3150 有点神，