整數快速冪（取模）、矩陣快速冪及其應用

摘要：

　　本文主要介紹了整數快速冪、矩陣快速冪及其應用，以題爲例重點展現了使用細節。

---

 

　　咱們要計算一個整數x的n次方，即x^n，普通的方法是連乘，這裏介紹一種效率很是高的計算冪運算的算法——反覆平方法。

　　首先考慮加速冪運算的方法，若是n=2^k，則能夠將x^n = ((x²)²)..，即只要作k次平方運算就能夠求得x^n。而後由此咱們能夠想到，先將n表示爲2的冪次之和，即x^n = 2^k1 + 2^k2 + 2^k3... ，那麼 x^n = x^{2^k1} * x^{2^k2}  * x^{2^k1} ...，只需在求x^{2^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", 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", n); 40         else
41             printf("%lld is normal.\n", 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項斐波那契數

　　根據斐波那契數的定義 F₀ = 0，F₁ = 1;

　　　　　　　　　           F_n = F_{n - 1} + F_{n - 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\n", (st.x[0][0] - 1 + mod) % mod); 51  } 52     return 0; 53 }

　　關於矩陣快速冪還有其餘一些重要的應用，時間有限，以後再作補充。

　　關於矩陣快速冪的介紹和應用舉例就到這裏，主要運用線性代數的知識，作題的時候要找到合適的遞推式，而後利用矩陣快速冪優化。