矩阵快速幂

目的：快速计算矩阵A的b次幂，即计算A^b;

代码如下：

typedef struct
{
    int m[MAXN][MAXN];
}Matrix;
Matrix a,per;
int n,M;
void init()///给矩阵赋初值
{
    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        scanf("%d",&a.m[i][j]);
        a.m[i][j]%=M;
        per.m[i][j]=(i==j);///如果i==j,返回1,否则，返回0.
     }///per.m[][]为对角线元素为1，其余为0的矩阵
}
Matrix multi(Matrix a,Matrix b)
{
    Matrix c;
    int i,j,k;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        c.m[i][j]=0;
        for(k=0;k<n;k++)
        {
            c.m[i][j]+=a.m[i][k]*b.m[k][j];
        }
        c.m[i][j]%=M;
    }
    return c;
}
Matrix power(int k)///矩阵的快速幂
{
    Matrix c,p;
    Matrix ans=per;///p第一次与ans相乘后，得到的结果还是p本身不变
    p=a;
    while(k)///二分法
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        k>>=1;
        p=multi(p,p);
    }
    return ans;
}

例题一  Matrix  Power Series(pku 3233 )

给定一个n*n的矩阵和一个整数k，要求计算S=A+A^2+A^3.......+A^k;

输入：输入仅包含一组测试数据，第一行包含三个整数n(n<=30),k(k<=10^9),m(m<10^4);接下来的n行每行包含n个n个小于32 768的非负整数位矩阵n*n个元素

输出：S%m；

代码如下：

#include<stdio.h>///把G++改成C++就AC
#include<iostream>
#include<string.h>
#define MAXN 101
using namespace std;
typedef struct///自定义结构体
{
    int m[MAXN][MAXN];
}Matrix;
Matrix a,per;
int n,m;
void init()///初始化
{
     int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        scanf("%d",&a.m[i][j]);
         a.m[i][j]%=m;
        per.m[i][j]=(i==j);
    }
}
Matrix add(Matrix a,Matrix b)///矩阵相加
{
    Matrix c;
    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        c.m[i][j]=(a.m[i][j]+b.m[i][j])%m;
    }
    return c;
}
Matrix multi(Matrix a,Matrix b)///矩阵相乘
{
    Matrix c;
    int k,i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        c.m[i][j]=0;
        for(k=0;k<n;k++)
        {
            c.m[i][j]+=a.m[i][k]*b.m[k][j];
        }
        c.m[i][j]%=m;
    }
    return c;
}
Matrix power(int k)///矩阵的k次幂
{
    Matrix c,p,ans=per;
    p=a;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        k>>=1;
        p=multi(p,p);
    }
    return ans;
}
Matrix MatrixSum(int k)///矩阵之和，递归求解
{              ///求的是A+A^2+A^3...A^k
    if(k==1)
        return a;
    Matrix temp,b;
    temp=MatrixSum(k/2);
    if(k&1)///若k为奇数.
    {
        b=power(k/2+1);
        temp=add(temp,multi(temp,b));
        temp=add(temp,b);
    }
    else
    {
        b=power(k/2);
        temp=add(temp,multi(temp,b));
    }
    return temp;
}
int main()
{
    int i,j,k;
    while(scanf("%d%d%d",&n,&k,&m)!=EOF)
    {
        init();
        Matrix ans=MatrixSum(k);
        for(i=0;i<n;i++)
        {
            for(j=0;j<n-1;j++)
                printf("%d ",ans.m[i][j]);
            printf("%d\n",ans.m[i][j]);
        }
    }
    return 0;
}