矩阵二分快速幂优化dp

思路就是先写出dp的状态转移方程，然后再把它转成形如 [f[i], f[i-1]]= A * [f[i-1], f[i-2]] 的形式，（尤其是看到要求的i值取值范围很大时，明白不能普通dp了不然容易超时）然后由“矩阵二分快速幂”来快速直接求得f[i]，而不用像普通dp一样一个一个向上做。

（2018/3/31）又默码了一遍，发现有些细节不得不提——首先是相乘是一定要注意顺序的，比如在main函数里将乘好次方的结果和初始矩阵相乘时，初始矩阵要在右边。第二点是我自己写的时候矩阵维度不像下面一样给你定了，而是作为了函数参数，那么就应该清晰一点（因为在matrix_pow函数中惯性思维是用次数n来for循环矩阵快速幂，而我粗心将传入的矩阵维度参数声明为n，就搞错了）。其次，要记住在矩阵相乘的时候把系数（对应的地方）都看好，然后要记得取模。最后，要记得用long long！因为既然你用到矩阵二分快速幂了，那么一定就是大数据才会去用的，那么记得次数n用long long！

btw，在手写推导矩阵时，注意如果遇到什么常数相加还是有办法解决的，毕竟常数*0=0，*1=本身，比如

那么A=【a b c 1 0 0 0 0 1】，init=【A(i-1)  A(i-2)  1】 

int n; // 所有矩阵都是 n * n 的矩阵
struct matrix {
   int a[100][100];    //这个维度倒无所谓，根据实际维度定也可以，不会这么大的 - - 因为是你手写推出来的矩阵。
};
matrix matrix_mul(matrix A, matrix B, int mod) {
    // 2 个矩阵相乘
    matrix C;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            C.a[i][j] = 0;
            for (int k = 0; k < n; ++k) {
                C.a[i][j] += A.a[i][k] * B.a[k][j] % mod;
                C.a[i][j] %= mod;
            }
        }
    }
    return C;
}
matrix unit() {
    // 返回一个单位矩阵
    matrix res;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (i == j) {
                res.a[i][j] = 1;
            } else {
                res.a[i][j] = 0;
            }
        }
    }
    return res;
}
matrix matrix_pow(matrix A, int n, int mod) {
    // 快速求矩阵 A 的 n 次方
    matrix res = unit(), temp = A;
    for (; n; n /= 2) {     //或者写 n>>=1 注意要加等号！
        if (n & 1) {
            res = matrix_mul(res, temp, mod);
        }
        temp = matrix_mul(temp, temp, mod);
    }
    return res;
}

应用：

我的AC代码

#include<iostream>
#include<bits/stdc++.h>
using namespace std;
const int maxn=5;
struct matrix
{
    long long m,n;
    long long a[maxn][maxn];
};
long long mod;
matrix mul(matrix A, matrix B) //矩阵乘法
{
    matrix C;
    C.m=A.m;
    C.n=B.n;
    for(int i=0;i<C.m;i++)
    {
        for(int j=0;j<C.n;j++)
        {
            C.a[i][j]=0;
            for(int k=0;k<A.n;k++)
            {
                C.a[i][j]+=A.a[i][k]*B.a[k][j]%mod;
                C.a[i][j]%mod;
            }
        }
    }
    return C;
}

matrix res(int m,int n)  //单位阵
{
    matrix E;
    E.m=m;
    E.n=n;
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<n;j++)
        {
            E.a[i][j]=0;
            if(i==j)
                E.a[i][j]=1;
        }
    }
    return E;
}

int main()
{
    long long n;
    cin>>n>>mod;
    if(n==1 || n==2)
    {
    	cout<<1;
    	return 0;
	}
    matrix E=res(2,2);
    matrix A;
    A.m=2; A.n=2;
    long long a[2][2]={1,1,1,0};
    for(int i=0;i<A.m;i++)
        for(int j=0;j<A.n;j++)
            A.a[i][j]=a[i][j];
    n=n-2; //只需要A^(n-2)
    for(; n; n>>=1)
    {
        if(n&1)
        {
            E=mul(E,A);
        }
        A=mul(A,A);
    }
    matrix start=res(2,1);
    start.m=2; start.n=1;
    long long b[2][1]={1,1};
    for(int i=0;i<start.m;i++)
        for(int j=0;j<start.n;j++)
            start.a[i][j]=b[i][j];
    matrix result=mul(E,start);
    cout<<result.a[0][0]%mod;
    return 0;
}

最主要的是找到状态转移方程：

a[i]=b[i-1]*y/100 + a[i-1](1-x/100)

b[i]=a[i-1]*x/100 + b[i-1](1-y/100)

然后把它写成[a[i]; b[i]] = A * [a[i-1];b[i-1]] ，写出这个2*2的矩阵A，就可以套模板了。

#include<iostream>
#include<bits/stdc++.h>
using namespace std;
const int maxn=5;
struct matrix
{
    double m,n;
    double a[maxn][maxn];
};
matrix mul(matrix A, matrix B) //矩阵乘法
{
    matrix C;
    C.m=A.m;
    C.n=B.n;
    for(int i=0;i<C.m;i++)
    {
        for(int j=0;j<C.n;j++)
        {
            C.a[i][j]=0;
            for(int k=0;k<A.n;k++)
            {
                C.a[i][j]+=A.a[i][k]*B.a[k][j];
            }
        }
    }
    return C;
}

matrix res(int m,int n)  //单位阵
{
    matrix E;
    E.m=m;
    E.n=n;
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<n;j++)
        {
            E.a[i][j]=0;
            if(i==j)
                E.a[i][j]=1;
        }
    }
    return E;
}

int main()
{
    long long a,b,k;
    double x,y;
    cin>>a>>b>>x>>y>>k;
    matrix E=res(2,2);
    matrix A;
    A.m=2; A.n=2;
    double aa[2][2]={1.00-x/100.00,y/100.00,x/100.00,1-y/100.00};  //关键！！！！！
    for(int i=0;i<A.m;i++)
        for(int j=0;j<A.n;j++)
            A.a[i][j]=aa[i][j];
    int n=k; //需要A^(n)
    for(; n; n>>=1)
    {
        if(n&1)
        {
            E=mul(E,A);
        }
        A=mul(A,A);
    }
    matrix start=res(2,1);
    start.m=2; start.n=1;
    double bb[2][1]={a,b};
    for(int i=0;i<start.m;i++)
        for(int j=0;j<start.n;j++)
            start.a[i][j]=bb[i][j];
    matrix result=mul(E,start);
    printf("%f %f",result.a[0][0],result.a[1][0]);
    return 0;
}