矩阵快速幂

1.引入

//整数快速幂

为了引入矩阵的快速幂，我先看了求整数的幂，即x^n

例如求x^8  就是x*x*x*x*x*x*x*x

正常的运算方式是，x的值一个个往上乘上去，乘法运算运行7次

（x  x)(x x) (x x)(x x)

也可以采用这种运算方式，先进行乘法得到x^2再对x^2进行三次乘法。这种运算要明显比第一种情况要快

所以对于整数快速幂，也是结合了这种思想

（x^m)*(x^n)=x^(m+n)

  x^19=(x^16)(x^2)(x^1)

//整数快速幂
int QuickPow(int x,int N)
{
    int res = x;
    int ans = 1;
    while(N)
    {
        if(N&1)
        {
            ans = ans * res;
        }
        res = res*res;
        N = N>>1;
    }
    return ans;
}

我们用 x^19来走这个函数

举例：x^19   

19的二进制为：1 0 0 1 1

初始：

ans = 1 ; res = x ;

则10011最后1位是1，为奇数

ans = res *ans =x ;

res =res * res = x ^ 2 ;

然后右移一位，1001

最后一位为1，为奇数

ans = res*ans = x*(x^2) = x^3；

res = res*res = x^2*x^2 = x^4 ；

然后右移一位， 100

为偶数

res = res*res = x^4*x^4 = x^8；

然后右移一位，10

为偶数

res=res*res=x^8*x^8=x^16;

然后右移一位，1

为奇数

ans = ans*res = (x^3)*（x^16） = x^19；

res = res*res = x^32

res=x^m，m始终是与二进制位置上的权值是相互对应的。

当二进制为0时候，我们只让res*res使得幂指数乘2，对应下一个二进制权值；

当二进制为1时候，我们让ans=ans*res,乘上了x幂次。

2.矩阵快速幂

看完了整数快速幂之后，再来学习一下，矩阵快速幂，

设有一个n*n的方阵A。

给出一个数M，即计算矩阵A的M次幂 A^M

struct Matrix  //结构体，矩阵类型
{
  int m[maxn][maxn];
}ans,res;

//计算矩阵乘法的函数，参数是矩阵a，矩阵b和一个n，代表这两个矩阵是几阶方阵
Maxtrix Mul(Maxtrix A,Maxtrix B,int n)
{
  Matrix tmp;//定义一个临时的矩阵，存放A*B的结果
  for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
    tmp.m[i][j]=0;
  for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
     for(int k=1;k<=n;k++)
     tmp.m[i][j]+=A.m[i][k]*B.m[k][j];
  return tmp;
}

//快速幂算法，求矩阵res的N次幂
void QuickPower(int N,int n)
{
  //对于整数求幂来说，ans初始化为1，对于矩阵乘法来说，ans应初始化为单位阵
  //对于单位矩阵来说，任何矩阵A*E=A
  for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
    {
     if(i==j) ans.m[i][j]=1;
      else ans.m[i][j]=0;
    }
  while(N)
  {
    if(N&1)
      ans=Mul(ans,res);
    res=Mul(res,res);
    N=N>>1;
  }
}

这是简单的计算矩阵幂，在实际中的运用的话，

3.实际应用

我们以最为常见的斐波那契数列为例：
f[n]=f[n-1]+f[n-2]

f[0]=0,f[1]=1;

斐波那契数列是一个递推式求出来的

第n+1项由第n项和第n-1项递推而来

因此可以用矩阵表示：

  POJ3070/NYOJ148

  给出了矩阵相乘的定义，要你求出斐波那契的第n项对1e4取余。

#include <iostream>
#include <cstddef>
#include <cstring>
#include <vector>
using namespace std;
typedef long long ll;
const int mod=10000;
typedef vector<ll> vec;
typedef vector<vec> mat;
mat mul(mat &a,mat &b)//表示不会这样用，，，，
{
    mat c(a.size(),vec(b[0].size()));
    for(int i=0; i<2; i++)
    {
        for(int j=0; j<2; j++)
        {
            for(int k=0; k<2; k++)
            {
                c[i][j]+=a[i][k]*b[k][j];
                c[i][j]%=mod;
            }
        }
    }
    return c;
}
mat pow(mat a,ll n)
{
    mat res(a.size(),vec(a.size()));
    for(int i=0; i<a.size(); i++)
        res[i][i]=1;//单位矩阵；
    while(n)
    {
        if(n&1) res=mul(res,a);
        a=mul(a,a);
        n/=2;
    }
    return res;
}
ll solve(ll n)
{
    mat a(2,vec(2));
    a[0][0]=1;
    a[0][1]=1;
    a[1][0]=1;
    a[1][1]=0;
    a=pow(a,n);
    return a[0][1];//也可以是a[1][0];
}
int main()
{
    ll n;
    while(cin>>n&&n!=-1)
    {
        cout<<solve(n)<<endl;
    }
    return 0;
}