矩阵快速幂

Posted 2023-04-29 xj22yangyichen

矩阵乘法

定义矩阵乘法的运算规则如下

\\[A\\left[m\\right]\\left[n\\right] * B\\left[n\\right]\\left[p\\right] = C\\left[m\\right]\\left[p\\right] \\]

其中 \\(C\\left[i\\right]\\left[j\\right]\\) 等于 \\(A\\) 的第 \\(i\\) 行乘 \\(B\\) 的第 \\(j\\) 列，举例如下

\\[\\beginbmatrix 5&6\\\\7&8 \\endbmatrix * \\beginbmatrix 1&2&3\\\\4&5&6 \\endbmatrix = \\beginbmatrix 5*1+6*4&5*2+6*5&5*3+6*6 \\\\ 7*1+8*4&7*2+8*5&7*3+8*6 \\endbmatrix \\]

矩阵乘法的性质

• 满足分配率和结合律，不满足交换率
• \\(A\\) 的列数必须等于 \\(B\\) 的行数

用斐波那契数列引入矩阵快速幂

已知斐波那契数列

\\[f_n=f_n-1+f_n-2 \\]

其中

\\[\\begincases f_0=0\\\\f_1=1 \\endcases \\]

则有

\\[\\beginbmatrix f_n+1&f_n \\endbmatrix = \\beginbmatrix f_n&f_n-1 \\endbmatrix * \\beginbmatrix 1&1\\\\1&0 \\endbmatrix = \\beginbmatrix f_n-1&f_n-2 \\endbmatrix * \\beginbmatrix 1&1\\\\1&0 \\endbmatrix^2 = \\beginbmatrix f_1&f_0 \\endbmatrix * \\beginbmatrix 1&1\\\\1&0 \\endbmatrix^n \\]

求解广义斐波那契数列

广义斐波那契数列数列的公式

\\[f_n=a*f_n-1+b*f_n-2 \\]

则有

\\[\\beginbmatrix f_n+1&f_n \\endbmatrix = \\beginbmatrix f_n&f_n-1 \\endbmatrix * \\beginbmatrix a&1\\\\b&0 \\endbmatrix = \\beginbmatrix f_1&f_0 \\endbmatrix * \\beginbmatrix a&1\\\\b&0 \\endbmatrix^n \\]

求解系数含有0的递推函数

已知公式

\\[f_n=3*f_n-1+5*f_n-3+9*f_n-4 \\]

则有

\\[\\beginbmatrix f_n+3&f_n+2&f_n+1&f_n \\endbmatrix = \\beginbmatrix f_3&f_2&f_1&f_0 \\endbmatrix * \\beginbmatrix 3&1&0&0\\\\0&0&1&0\\\\5&0&0&1\\\\9&0&0&0 \\endbmatrix^n \\]

将系数为0的项对应的矩阵中的位置置0

求解一般递推函数

已知公式

\\[f_n=a*f_n-1+b*f_n-3+c \\]

则有

\\[\\beginbmatrix f_n+2&f_n+1&f_n&c \\endbmatrix = \\beginbmatrix f_2&f_1&f_0&c \\endbmatrix * \\beginbmatrix a&1&0&0\\\\0&0&1&0\\\\b&0&0&1\\\\1&0&0&0 \\endbmatrix^n \\]

将常数项加入待求的矩阵中一并运算，并将计算矩阵相对应的位置置1以确保常数项的值不变

求解有下标的递推函数

已知公式

\\[f_n=f_n-1+2*f_n-2+n^3 \\]

这里需要用到一个知识：立方的展开

\\[n^3=(n-1)^3+3(n-1)^2+3(n-1)+1 \\]

事实上，矩阵乘法求递归问题也可以用纵向的矩阵表示，如下

\\[\\beginbmatrix f_n\\\\f_n-1\\\\n^3\\\\n_2\\\\n\\\\1 \\endbmatrix = \\beginbmatrix f_2\\\\f_1\\\\8\\\\4\\\\2\\\\1 \\endbmatrix * \\beginbmatrix 1&2&1&3&3&1 \\\\ 1&0&0&0&0&0 \\\\ 0&0&1&3&3&1 \\\\ 0&0&0&1&2&1 \\\\ 0&0&0&0&1&1 \\\\ 0&0&0&0&0&1 \\endbmatrix^n-2 \\]

快速幂计算（整数快速幂/矩阵快速幂）

 

//*************快速幂计算****************************************

 

朴素算法实现：

 1 ll get_pow(ll x, ll n)  //** (这里的n要求不小于0，如果n小于0则令n=-n，并且最终返回1.0/ans即可)
 2 {
 3     ll ans=1;
 4     while(n--)
 5     {
 6         ans*=x%MAX;
 7         ans%=MAX;
 8     }
 9     return ans;
10 }

 

 

快速幂算法：

原理：

二分：

假设我们现在要计算pow(x,n），那么有当n为偶数时pow(x, n)==pow(x*x, n/2)，当n为奇数时，pow(x, n)==pow(x, n-1)*x, 此时n-1为偶数，可按前面的公式继续迭代；

循环往复，即可计算出答案；

 

二进制：

计算pow(x,n)，先将n转化为二进制形式，n=2^a+2^b....

例如：计算pow(x,21)，19=2^4+2^2+2^1；其中2^2可以由(2^1)*(2^1)得到，同理, 2^4可以由(2^2)*(2^2)得到；

所有有pow(x,21)==x^(2^4+2^2+2^1)，其时间复杂度为o(long2(n));

 

以上两种思路的出发点不同，不过其本质一致，代码也相同；

 

代码：

 1 ll get_pow(ll x, ll n)  //** (这里的n要求不小于0，如果n小于0则令n=-n，并且最终返回1.0/ans即可)
 2 {
 3     int ans=1;
 4     while(n)
 5     {
 6         if(n&1)
 7         {
 8             ans=(ans*x)%MAX;
 9         }
10         x=(x*x)%MAX;
11         n>>=1;
12     }
13     return ans;
14 }

 

 

//*********************矩阵快速幂计算**********************************

矩阵快速幂和快速幂算法原理一样，只是操作对象换成了矩阵；

http://acm.nyist.net/JudgeOnline/problem.php?pid=148（题目链接）

 

ac代码：

 1 #include <bits/stdc++.h>
 2 #define MAXN 2
 3 #define mod 10000
 4 #define ll long long
 5 using namespace std;
 6 
 7 struct Matrix
 8 {
 9     ll x[MAXN][MAXN];
10 };
11 
12 Matrix temp={1, 1, 1, 0};
13 
14 ll n;  //***n为幂数，结果对mod取模
15 
16 Matrix multi(Matrix a, Matrix b) //***矩阵乘法
17 {
18     Matrix c;
19     memset(c.x, 0, sizeof(c.x));
20     for(int i=0; i<MAXN; i++)
21     {
22         for(int j=0; j<MAXN; j++)
23         {
24             for(int k=0; k<MAXN; k++)
25             {
26                 c.x[i][j]+=(a.x[i][k]*b.x[k][j])%mod;
27             }
28         }
29     }
30     return c;
31 }
32 
33 ll Pow(ll n)
34 {
35     Matrix ans;
36     temp.x[0][0]=temp.x[0][1]=temp.x[1][0]=1;
37     temp.x[1][1]=0;
38     memset(ans.x, 0, sizeof(ans.x));
39     ans.x[0][0]=ans.x[1][1]=1;
40     while(n)
41     {
42         if(n&1) ans=multi(ans, temp);
43         temp=multi(temp, temp);
44         n>>=1;
45     }
46     return ans.x[0][1];
47 }
48 
49 int main(void)
50 {
51     std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
52     while(cin >> n && n!=-1)
53     {
54         cout << Pow(n)%mod << endl;
55     }
56     return 0;
57 }