数论常用内容——矩阵快速幂

在数论题以及一些非数论题中，经常出现递推的情况，如果找不到规律，强行迭代递推的话不是一个明智的选择，今天我就来给大家介绍一个适用于这种情形的方法——矩阵快速幂

矩阵快速幂是如何优化递推的呢？

首先，需要先构造合适的初始状态（第一个矩阵）然后，利用此矩阵和矩阵乘法的性质，使用快速幂的手段求出之后的状态

构造矩阵可以根据矩阵乘法的实现特点来构造，利用合适的性质可以简化运算

如何构造矩阵呢？

（一）f[n]=f[n-1]+f[n-2],f[1]=f[2]=1的第n项快速求法

考虑1×2的矩阵【f[n-2],f[n-1]】。根据Fibonacci数列的递推关系，我们可以通过乘以一个2×2的矩阵A，得到矩阵：【f[n-1],f[n]】。

即：【f[n-2],f[n-1]】*A = 【f[n-1],f[n]】=【f[n-1],f[n-1]+f[n-2]】

很容易构造出这个2×2矩阵A，即：
０ １
１ １

所以，有【f[1],f[2]】×A＝【f[2],f[3]】

又因为矩阵乘法满足结合律，故有：

【f[1],f[2]】×A ^(n-1) =【f[n],f[n+1]】

这个矩阵的第一个元素f[n]即为所求。

（二）f[n]=f[n-1]+f[n-2]+1,f[1]=f[2]=1的第n项的快速求法

仿照前例，考虑1×3的矩阵【f[n-2],f[n-1],1】，希望求得某3×3的矩阵A，使得此1×3的矩阵乘以A得到矩阵：【f[n-1],f[n],1】

即：【f[n-2],f[n-1],1】* A ＝【f[n-1],f[n],1】=【f[n-1],f[n-1]+f[n-2]+1,1】

容易构造出这个3×3的矩阵A，即：
０ １ ０
１ １ ０
０ １ １

故：【f[1],f[2],1】* A^(n-1) = 【f[n],f[n+1],1】

（三）数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1的第n项的快速求法

仿照前例，考虑1×4的矩阵【f[n-2],f[n-1],n,1】，希望求得某4×4的矩阵A，使得此1×4的矩阵乘以A得到矩阵：【f[n-1],f[n],n+1,1】

即：【f[n-2],f[n-1],n,1】* A  = 【f[n-1],f[n],n+1,1】＝【f[n-1],f[n-1]+f[n-2]+n+1,n+1,1】

容易构造出这个4×4的矩阵A，即：
０ １ ０ ０
１ １ ０ ０
０ １ １ ０
０ １ １ １

故：【f[1],f[2],3,1】* A^(n-1) = 【f[n],f[n+1],n+2,1】

（四）数列f[n]=f[n-1]+f[n-2],f[1]=f[2]=1的前n项和s[n]=f[1]+f[2]+……+f[n]求法

仿照之前的思路，考虑1×3的矩阵【f[n-2],f[n-1],s[n-2]】，我们希望通过乘以一个3×3的矩阵A，得到1×3的矩阵：【f[n-1],f[n],s[n-1]】

即：【f[n-2],f[n-1],s[n-2]】 * A  = 【f[n-1],f[n],s[n-1]】=【f[n-1],f[n-1]+f[n-2],s[n-2]+f[n-1]】

容易得到这个3×3的矩阵A是：
０ １ ０
１ １ １
０ ０ １

这种方法的矩阵规模是(r+1)*(r+1)

f(1)=f(2)=s(1)=1 ，所以，有【f(1),f(2),s(1)】* A  = 【f(2),f(3),s(2)】

故：【f(1),f(2),s(1)】* A^(n-1)  = 【f(n),f(n+1),s(n)】

（五）数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1的前n项和s[n]=f[1]+f[2]+……+f[n]求法

考虑1×5的矩阵【f[n-2],f[n-1],s[n-2],n,1】

我们需要找到一个5×5的矩阵A，使得它乘以A得到如下1×5的矩阵【f[n-1],f[n],s[n-1],n+1,1】

即：【f[n-2],f[n-1],s[n-2],n,1】* A  =【f[n-1],f[n],s[n-1],n+1,1】=【f[n-1], f[n-1]+f[n-2]+n+1,s[n-2]+f[n-1],n+1,1】

容易构造出A为：
０ １ ０ ０ ０
１ １ １ ０ ０
０ ０ １ ０ ０
０ １ ０ １ ０
０ １ ０ １ １

故：【f(1),f(2),s(1),3,1】* A^(n-1)  = 【f(n),f(n+1),s(n),n+2,1】

（六）f[n]=p*f[n-1]+q*f[n-2]+r*n+s

可以构造矩阵A为：
0  q  0  0  0
1  p  1  0  0
0  0  1  0  0
0  r  0  1  0
0  s  0  1  1

（七）f[n]=Sigma(a[n-i]*f[n-i])+Poly(n)

其中0<i<=某常数c, Poly (n)表示n的多项式，我们依然可以构造类似的矩阵A来解决问题。

设Degree(Poly(n))=d, 并规定Poly(n)=0时，d=-1，此时对应于常系数线性齐次递推关系。则本方法求前n项和的复杂度为：((c+1)+(d+1))3*logns

例如：A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2)；

给定三个值N，X，Y求S(N):S(N) = A(0)2 +A(1)2+……+A(n)2。

解：
考虑1*4 的矩阵【s[n-2],a[n-1]^2,a[n-2]^2,a[n-1]*a[n-2]】

我们需要找到一个4×4的矩阵A，使得它乘以A得到1×4的矩阵
【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】

即：【s[n-2],a[n-1]^2,a[n-2]^2,a[n-1]*a[n-2]】* A = 【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】= 【s[n-2]+a[n-1]^2 , x^2 * a[n-1]^2 + y^2 * a[n-2]^2 + 2*x*y*a[n-1]*a[n-2]，a[n-1]^2 , x*a[n-1]^2 + y*a[n-2]a[n-1]】

可以构造矩阵A为：
1     0    0    0
1    x^2   1    x
0    y^2   0    0
0    2xy   0    y

故：【S[0],a[1]^2,a[0]^2,a[1]*a[0]】 * A^(n-1) = 【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】

所以：【S[0],a[1]^2,a[0]^2,a[1]*a[0]】 * A^(n) = 【s[n],a[n+1]^2,a[n]^2,a[n+1]*a[n]】

若A = (B * C ) 则AT = ( B * C )T = CT * BT

故可得如下图所示的矩阵

构造完矩阵如何得到答案呢？

这里，矩阵快速幂，和数字的快速幂算法是一样的，只是把数字换成了矩阵而已，矩阵的乘法线性代数应该都学过，这里不再赘述，下面是我写的矩阵快速幂代码实现

const int MAXN = 15;//最大294
const int mod = 10000007;

//矩阵定义
struct  Matrix{
    __int64 a[MAXN][MAXN];
    int r;//行数
    int c;//列数
};

//输出测试函数
void print(Matrix mat){
    for(int i=0;i<mat.r;i++){
        for(int j=0;j<mat.c;j++){
            if(j) cout<<" ";
            cout<<mat.a[i][j];
        }
        cout<<endl;
    }
    cout<<endl;
}

//矩阵初始化
void init(Matrix &mat, __int64 n){
    for (int i = 0;i < mat.r;i++){
        for (int j = 0;j < mat.c;j++){
            mat.a[i][j] = n;
        }
    }
}

//矩阵乘法
Matrix mul(Matrix mat1, Matrix mat2){
    Matrix res;
    res.c = mat2.c;
    res.r = mat1.r;
    init(res, 0);
    for (int i = 0;i < mat1.r;i++){
        for (int j = 0;j < mat2.c;j++){
            for (int k = 0;k < mat1.c;k++){
                res.a[i][j] += mat1.a[i][k] * mat2.a[k][j];
                res.a[i][j] %= mod;
            }
        }
    }
    return res;
}

//矩阵快速幂
Matrix fmp(Matrix ans,__int64 n){
    Matrix base;
    base.c=base.r=ans.r;
    init(base,0);
    for(int i=0;i<base.r;i++){
        base.a[i][i]=1;
    }
    while (n){
        if (n & 1) base = mul(base, ans);
        ans = mul(ans, ans);
        n >>= 1;
    }
    return base;
}