整数快速幂(取模)矩阵快速幂及其应用

Posted wenzhixin

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了整数快速幂(取模)矩阵快速幂及其应用相关的知识,希望对你有一定的参考价值。

摘要:

  本文主要介绍了整数快速幂、矩阵快速幂及其应用,以题为例重点展示了使用细节。


 

  我们要计算一个整数x的n次方,即x^n,普通的方法是连乘,这里介绍一种效率非常高的计算幂运算的算法——反复平方法。

  首先考虑加速幂运算的方法,如果n=2^k,则可以将x^n = ((x2)2)..,即只要做k次平方运算就可以求得x^n。然后由此我们可以想到,先将n表示为2的幂次之和,即x^n = 2k1 + 2k2 + 2k3... ,那么 x^n = x2^k1 * x2^k2  * x2^k1 ...,只需在求x2^i 的同时进行计算就好了。最终得到O(logn)的计算幂运算的算法。

  比如计算x^22 = x^16 * x^4 * x^2,其中22的二进制数是10110,也就是需要反复平方3次。代码如下:

 1 typedef long long ll;
 2 ll qpow(ll x, ll n) {
 3     ll res = 1;
 4     while(n) {
 5         if(n&1)
 6             res = res * x;    //如果二进制最低位为1,则乘上x^(2^i) 
 7         x = x * x;            //将x平方 
 8         n >>= 1;             //n/2
 9     }
10     return res;
11 }

  在实际应用中有时还需要求解x^n%mod。代码如下:

 1 typedef long long ll;
 2 ll qpow(ll x, ll n, ll mod) {
 3     ll res = 1;
 4     while(n) {
 5         if(n&1)
 6             res = res * x % mod;    //如果二进制最低位为1,则乘上x^(2^i) 
 7         x = x * x % mod;            //将x平方 
 8         n >>= 1;                    //n/2
 9     }
10     return res;
11 }

  看一道例题:UVA 10006 Carmichael Numbers

  判断是否是C数,需要满足以下两个条件

  1.不是素数.

  2.对任意的1<x<n都有x^n和x同余模n.

  代码如下:

 1 #include <cstdio>
 2 #include <cmath>
 3 typedef long long ll;
 4 
 5 ll qpow(ll x, ll n, ll mod) {
 6     ll res = 1;
 7     while(n) {
 8         if(n&1)
 9             res = res * x % mod;
10         x = x * x % mod;
11         n >>= 1; 
12     }
13     return res;
14 }
15 bool isprime(ll x) {
16     if(x == 0 || x == 1)
17         return 0;
18     ll k = (ll)sqrt(x);
19     for(ll i = 2; i < k; i++) {
20         if(x % i == 0)
21             return 0;
22     }    
23     return 1;
24 }
25 int main()
26 {
27     ll n;
28     while(scanf("%lld", &n) == 1 && n != 0) {
29         if(isprime(n)) {
30             printf("%lld is normal.
", n);
31             continue;
32         }
33         ll i;
34         for(i = 2; i < n; i++) {
35             if(qpow(i, n, n) != i % n)
36                 break;
37         }
38         if(i == n) 
39             printf("The number %lld is a Carmichael number.
", n);
40         else
41             printf("%lld is normal.
", n);
42     }
43     return 0;
44 }

  现在要求一个矩阵A的m次幂,也就是A^m,首先应该会两个矩阵的乘法,然后知道A^m的结果一定是一个同型矩阵,最后需要理解上面的整数快速幂。剩下的就是将整数换成矩阵操作。代码如下:

 1 const int Matrix_size = 2
 2 struct Matrix {//定义矩阵 
 3     int x[Matrix_size][Matrix_size]; 
 4 }ans, res;
 5 /*
 6 定义矩阵乘法,A * B, 它们的都是n阶方阵 
 7 */ 
 8 Matrix Mmul(Matrix A, Matrix B, int n) {
 9     Matrix tmp;
10     for(int i = 1; i <= n; i++) {
11         for(int j = 1; j <= n; j++) {
12             tmp.x[i][j] = 0;
13         }
14     } 
15     
16     for(int i = 1 ; i <= n; i++) {
17         for(int j = 1; j <= n; j++) {
18             for(int k = 1; k <= n; k++) {
19                 tmp.x[i][j] += A.[i][k] * B[k][j];
20             }
21         }
22     }
23     return tmp;
24 }
25 /*
26 求res的N次方,n是res的阶数 
27 */
28 void qmpow(int N, int n) {
29     for(int i = 1; i <= n; i++) {
30         for(int j = 1; j <= n; j++) {
31             res.x[i][j] = i == j ? 1 : 0;
32         }
33     }
34     
35     while(N) {
36         if(N&1)
37             ans = Mmul(ans, res);
38         res = Mmul(res, res);
39         N >>= 1;
40     }
41 }

   利用上面的矩阵快速幂算法可以快速的求解一个矩阵的n次幂,那么求一个矩阵的n次幂有什么用呢?

  1.求第n项斐波那契数

  根据斐波那契数的定义 F0 = 0,F1 = 1;

                    F= Fn - 1 + Fn - 2(n>=2).

  可以用矩阵表示为:

技术分享图片

  进一步递推得到:

技术分享图片

  这里需要求的是右边系数矩阵的n-2次幂,然后代入前两项即可求得f(n),也就是第n项斐波那契数。

下面看一道例题:HDU 6198 number number number

题意

先给出了斐波那契数列的定义

 F0 = 0, F1 = 1;

 Fn = F n - 1 + F n - 2.

再给出坏数的定义,给一个整数k,如果一个整数n不能有k个任意的(可重复)的斐波那契数组成,就成为是一个坏数。现给出k,问最小的坏数是多少,答案对998244353取模。

解题思路

应暴力的方法是不行了,因为给出的k很大。先观察前几项可得:

当k = 1时,4 = 5 - 1 = F(2 * 1 + 3) - 1;

当k = 2时,12 = 13 - 1  = F(2 * 2 + 3) - 1;

问题变成了求解第2 * k + 3项斐波那契数,又因为k很大,就需要使用矩阵快速幂解决了。

根据定义我们定义前两项是1和1,系数矩阵是1,1,1,0,求第n项需要计算系数矩阵的n-2次幂,然后乘上前两项,得到F(n)和F(n - 1)。

代码如下:

 1 #include <cstdio>
 2 #include <cstring>
 3 typedef long long ll;
 4 const int mod = 998244353;
 5 struct Matrix {
 6     ll x[2][2];
 7 };
 8 Matrix Mmul(Matrix a, Matrix b) {
 9     Matrix tmp;
10     memset(tmp.x, 0, sizeof(tmp.x));
11     for(int i = 0; i < 2; i++) {
12         for(int j = 0; j < 2; j++) {
13             for(int k = 0; k < 2; k++) {
14                 tmp.x[i][j] = (tmp.x[i][j] + a.x[i][k] * b.x[k][j] % mod) % mod;
15             }
16         }
17     }
18     return tmp;
19 }
20 Matrix Mqpow(Matrix a, ll n) {
21     Matrix tmp;
22     for(int i = 0; i < 2; i++) {
23         for(int j = 0; j < 2; j++) {
24             tmp.x[i][j] = i == j ? 1 : 0;
25         }
26     }
27     
28     while(n) {
29         if(n&1)
30             tmp = Mmul(tmp, a);
31         a = Mmul(a, a);
32         n >>= 1;
33     }
34     return tmp;
35 }
36  int main()
37 {
38     ll k;
39     while(scanf("%lld", &k) != EOF) {
40         Matrix st;
41         st.x[0][0] = 1; st.x[0][1] = 1;
42         st.x[1][0] = 1; st.x[1][1] = 0;
43         
44         Matrix init;
45         init.x[0][0] = 1; init.x[0][1] = 0;
46         init.x[1][0] = 1; init.x[1][1] = 0;
47         
48         st = Mqpow(st, 2 * k + 1);
49         st = Mmul(st, init);
50         printf("%lld
", (st.x[0][0] - 1 + mod) % mod);
51     }    
52     return 0;    
53 }

  关于矩阵快速幂还有其他一些重要的应用,时间有限,之后再做补充。

  关于矩阵快速幂的介绍和应用举例就到这里,主要运用线性代数的知识,做题的时候要找到合适的递推式,然后利用矩阵快速幂优化。

以上是关于整数快速幂(取模)矩阵快速幂及其应用的主要内容,如果未能解决你的问题,请参考以下文章

递推求值快速幂矩阵

codevs1281 矩阵乘法 快速幂 !!!手写乘法取模!!! 练习struct的构造函数和成员函数

Happy Necklace (矩阵快速幂 + 递推 + 取模)

快速幂介绍及其模板

Jzzhu and Sequences (矩阵快速幂 + 取模)

整数快速乘法/快速幂+矩阵快速幂+Strassen算法 (转)