数论-矩阵快速幂
一道题：
求斐波那契的第i项对\(10^9+7\)求余的值（i<\(10^{10}\)）
- 算法一： 递推，会TLE
- 算法二： 你偷偷百度了一下 找到了fib的通项公式 \[fib(n)=\frac {1}{\sqrt5} \left[\left(\frac{1+\sqrt5}{2}\right)^n-\left(\frac{1-\sqrt5}{2}\right)^n\right]\] 但是这个无法求余，并且数字超过了double的有效位数
这时候只能使用矩阵快速了 考虑如下 \[\begin{bmatrix} 1 & 1\\ 1 &0 \\ \end{bmatrix} * \begin{bmatrix} 1\\ 1\\ \end{bmatrix}= \begin{bmatrix} 2\\ 1\\ \end{bmatrix} \] 也就是 \[\begin{bmatrix} 1 & 1\\ 1 &0 \\ \end{bmatrix} * \begin{bmatrix} fib(2)\\ fib(1)\\ \end{bmatrix}= \begin{bmatrix} fib(1)+fib(2)\\ 0+fib(2)\\ \end{bmatrix}= \begin{bmatrix} fib(3)\\ fib(2) \end{bmatrix} \] 我们继续往下写 \[\begin{bmatrix} 1 & 1\\ 1 & 0\\ \end{bmatrix} * \begin{bmatrix} 2\\ 1\\ \end{bmatrix}= \begin{bmatrix} 3\\ 2\\ \end{bmatrix} \] 也就是 \[\begin{bmatrix} 1 & 1\\ 1 &0\\ \end{bmatrix} * \begin{bmatrix} fib(3)\\ fib(2)\\ \end{bmatrix}= \begin{bmatrix} fib(2)+fib(3)\\ 0+fib(3)\\ \end{bmatrix}= \begin{bmatrix} fib(4)\\ fib(3) \end{bmatrix} \]
于是我们得出结论
\[ \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}^n* \begin{bmatrix} 1\\ 1 \end{bmatrix}= \begin{bmatrix} fib(n+2)\\ fib(n+1) \end{bmatrix} \]
至此，最难解决的问题就是如何进行矩阵相乘了 有一个n次方我想到使用快速幂可以求解 所以我们先要解决矩阵乘法*
矩阵的构建与赋值
//构建矩阵结构体
struct maix{
int n,m;
int mat[50][50];
};
//使用这种方法赋值 数字会按结构体定义顺序赋值n,m,mat,尤其注意数组赋值
maix a={2,2,{{},{0,1,1},{0,1,0}}};
maix b={2,1,{{},{0,1},{0,1}}};a矩阵为(从下标1开始) \[ \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix} \] ## 矩阵乘法
//乘法
maix mul(maix x,maix y){
maix res;res.n=x.n;res.m=y.m;
for(int i=1;i<=x.n;i++){
for(int k=1;k<=y.m;k++){
res.mat[i][k]=0;
for(int j=1;j<=y.n;j++)
res.mat[i][k]+=x.mat[i][j]*y.mat[j][k];
}
}
return res;
}然后就是快速幂了 1. 错误示范： 我只想求前面这个矩阵的n次方，于是写了函数
maix matrix_pow(maix x,int n){
maix res;
while(n){
if(n&1)y=mul(res,x);
x=mul(x,x);
n>>=1;
}
return y;
}有一个问题，就是我们在普通快速幂中可以做一个res用来放结果，因为我们可以随意创建一个数值1，但是在这里，我们无法创建一个同大小的单位矩阵 所以这里的解决办法是把后面那个 \[ \begin{bmatrix} fib(2)\\ fib(1) \end{bmatrix} \] 作为res，顺便也处理掉，于是我们写下正确代码
maix maix_pow(maix x,maix y,int n){
while(n){
if(n&1)y=mul(x,y); //矩阵乘法左右千万别错
x=mul(x,x);
n>>=1;
}
return y;
}输出矩阵(调试用，最后做题直接访问数组)
void print_maix(maix res){
cout<<"PR n="<<res.n<<" "<<res.m<<endl;
for(int i=1;i<=res.n;i++){
for(int j=1;j<=res.m;j++)
cout<<res.mat[i][j]<<" ";
cout<<endl;
}
cout<<"=======================\n";
}完整的代码
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
struct maix{
int n,m;
int mat[50][50];
};
void print_maix(maix res){
cout<<"PR n="<<res.n<<" "<<res.m<<endl;
for(int i=1;i<=res.n;i++){
for(int j=1;j<=res.m;j++)
cout<<res.mat[i][j]<<" ";
cout<<endl;
}
cout<<"=======================\n";
}
maix mul(maix x,maix y){
maix res;res.n=x.n;res.m=y.m;
for(int i=1;i<=x.n;i++){
for(int k=1;k<=y.m;k++){
res.mat[i][k]=0;
for(int j=1;j<=y.n;j++)
res.mat[i][k]+=x.mat[i][j]*y.mat[j][k];
}
}
return res;
}
maix maix_pow(maix x,maix y,int n){
while(n){
if(n&1)y=mul(x,y);
x=mul(x,x);
n>>=1;
}
return y;
}
int main(){
maix a={2,2,{{},{0,1,1},{0,1,0}}};
maix b={2,1,{{},{0,1},{0,1}}};
print_maix(a);
print_maix(b);
print_maix(maix_pow(a,b,10));
return 0;
}矩阵快速幂的应用
常用于求递推数列
构造方法
一般的我们会通过递推式构造矩阵，一般多项式有几项我们构造几乘几的矩阵(包括常数)
1. 数列f[n]=f[n-1]+f[n-2]+1,f[1]=f[2]=1的第n项的快速求法（不考虑高精度）
考虑3x1的矩阵 \[ \begin{bmatrix} f(n-2)\\ f(n-1)\\ 1 \end{bmatrix} \] ，希望求得某3×3的矩阵A，使得此1×3的矩阵乘以A得到矩阵: \[ \begin{bmatrix} f(n-1)\\ f(n)\\ 1 \end{bmatrix} \] 即： \[ A*\begin{bmatrix} f(n-2)\\ f(n-1)\\ 1 \end{bmatrix}= \begin{bmatrix} f(n-1)\\ f(n)\\ 1 \end{bmatrix} \] 构造矩阵 \[ \begin{bmatrix} 0&1&0\\ 1&1&1\\ 0&0&1\\ \end{bmatrix} \] 故： \[ \begin{bmatrix} 0&1&0\\ 1&1&1\\ 0&0&1\\ \end{bmatrix}^n* \begin{bmatrix} f(n-2)\\ f(n-1)\\ 1 \end{bmatrix}= \begin{bmatrix} f(n-1)\\ f(n)\\ 1 \end{bmatrix} \]
2.数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1
考虑4x1的矩阵 \[ \begin{bmatrix} f(n-2)\\ f(n-1)\\ n\\ 1 \end{bmatrix} \] 希望求得某4×4的矩阵A，使得Ac乘以此矩阵得到矩阵： \[ \begin{bmatrix} f(n-1)\\ f(n)\\ n\\ 1 \end{bmatrix} \] 即： \[ A*\begin{bmatrix} f(n-2)\\ f(n-1)\\ n\\ 1 \end{bmatrix}= \begin{bmatrix} f(n-1)\\ f(n-1)+f(n-2)\\ n+1\\ 1 \end{bmatrix} \] 容易构造出这个4×4的矩阵A，即： \[ \begin{bmatrix} 0&1&0&0\\ 1&1&1&1\\ 0& 0&1&1\\ 0&0&0&1\\ \end{bmatrix} \] 故： \[ \begin{bmatrix} 0&1&0&0\\ 1&1&1&1\\ 0& 0&1&1\\ 0&0&0&1\\ \end{bmatrix}^n*\begin{bmatrix} f(n-2)\\ f(n-1)\\ n\\ 1 \end{bmatrix}= \begin{bmatrix} f(n-1)\\ f(n-1)+f(n-2)\\ n+1\\ 1 \end{bmatrix} \]
3.数列f[n]=f[n-1]+f[n-2],f[1]=f[2]=1的前n项和s[n]=f[1]+f[2]+……+f[n]
仿照之前的思路，考虑3x1的矩阵 \[ \begin{bmatrix} f(n-2)\\ f(n-1)\\ s(n-2) \end{bmatrix} \] ，我们希望通过乘以一个3×3的矩阵A，得到1×3的矩阵: \[ \begin{bmatrix} f(n-1)\\ f(n)\\ s(n-1) \end{bmatrix} \] 即： \[ A* \begin{bmatrix} f(n-2)\\ f(n-1)\\ s(n-2) \end{bmatrix}= \begin{bmatrix} f(n-1)\\ f(n)\\ s(n-1) \end{bmatrix} \] 容易得到这个3×3的矩阵A是： \[ \begin{bmatrix} 0&1&0\\ 0&1&0\\ 0&1&1\\ \end{bmatrix} \] 这种方法的矩阵规模是\((r+1)^2\) \(f(1)=f(2)=s(1)=1\)，所以，有 \[ \begin{bmatrix} f(1)\\ f(2)\\ s(1) \end{bmatrix}= \begin{bmatrix} f(2)\\ f(3)\\ s(2) \end{bmatrix} \] 故： \[ \begin{bmatrix} 0&1&0\\ 0&1&0\\ 0&1&1\\ \end{bmatrix}^{n-1}* \begin{bmatrix} f(1)\\ f(2)\\ s(1) \end{bmatrix}= \begin{bmatrix} f(2)\\ f(3)\\ s(2) \end{bmatrix} \]
4.数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1的前n项和s[n]=f[1]+f[2]+……+f[n]
考虑5×1的矩阵 \[ \begin{bmatrix} f(n-2)\\ f(n-1)\\ s(n-2)\\ n\\ 1 \end{bmatrix} \] 我们需要找到一个5×5的矩阵A，使得它乘以A得到如下1×5的矩阵 \[ \begin{bmatrix} f(n-1)\\ f(n)\\ s(n-1)\\ n+1\\ 1 \end{bmatrix} \] 即：【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为： \[ \begin{bmatrix} 0&1&0&0&0\\ 1&1&1&0&0\\ 0&1&1&0&0\\ 0&0&0&1&1\\ 0&0&0&0&1\\ \end{bmatrix} \] 故：【f(1),f(2),s(1),3,1】* A^(n-1) = 【f(n),f(n+1),s(n),n+2,1】 \[ \begin{bmatrix} 0&1&0&0&0\\ 1&1&1&0&0\\ 0&1&1&0&0\\ 0&0&0&1&1\\ 0&0&0&0&1\\ \end{bmatrix}^{n-1}* \begin{bmatrix} f(1)\\ f(2)\\ f(1)\\ f(3)\\ 1 \end{bmatrix}= \begin{bmatrix} f(n)\\ f(n+1)\\ s(n)\\ n+2\\ 1 \end{bmatrix} \] 一般地，如果有\(f[n]=p*f(n-1)+q*f(n-2)+r*n+s\) 可以构造矩阵A为： \[ \begin{bmatrix} 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\\ \end{bmatrix} \] 更一般的，对于\(f[n]=\Sigma(a[n-i]*f[n-i])+Poly(n)\)，其中\(0<i\leq某常数c\), \(Poly (n)\)表示n的多项式，我们依然可以构造类似的矩阵\(A\)来解决问题。 设\(Degree(Poly(n))=d,\) 并规定\(Poly(n)=0\)时，\(d=-1\)，此时对应于常系数线性齐次递推关系。则本方法求前n项和的复杂度为： \[O((c+1)+(d+1))3*log(n)\]