數論篇2——快速冪全解（模平方根算法）

模平方根算法

求a的b次方有庫函數 pow(a, b)，可是它返回值是double類型，而且在不同開發環境下，數據有精度誤差（比如某DEV，詳見），如果自己寫for循環，當b特別大時，超范圍、超時都妥妥的。所以，就有了模平方根算法，也就是通常說的快速冪。

原理：

//遞歸寫法
int pow_power(int a, int b,int MOD){//a的b次方
    if(b == 0) return 1;
    int res = pow_power(a, b/2);
    res = res * res % MOD;
    if(b&1) res = res * a % MOD;
    return res;
}
//迭代寫法
int quickPower(int x, int n, int mod) {
    int t = x, res = 1;
    while (n) {
        if (n & 1) 
            res = ((res % mod) * (t % mod)) % mod;
        t = ((t % mod) * (t % mod)) % mod;
        n >>= 1;
    }

    return res % mod;
}

根據原理，還可以寫出來快速乘（龜速乘），在乘法會爆long long范圍時，只能這樣做了。

int mul(int a, int b, int p){//快速乘，計算a*b%p 
    int res = 0;
    while(b){
        if(b & 1) res = (ret + a) % p;
        a = (a + a) % p;
        b >>= 1;
    }
    return res;
}

高精度快速冪

上述算法，如果數據級別較大，改用long long，但當需要處理大於10^10的數據時，依然需要實現高精度快速冪。

一個簡單的板子，沒有做取模，位數限制在了500位，根據題目要求可靈活修改。

#include <iostream>
#include <string.h>
#include <math.h>

using namespace std;

int a[1001],b[1001], res[1001], temp[1001];

void carry(int* arr) {
    for (int i = 0; i < 500; i++) {
        arr[i + 1] += arr[i] / 10;
        arr[i] %= 10;
    }
}

void multi(int* a, int* b, int* res) {
    memset(temp, 0, sizeof(temp));
    for (int i = 0; i < 500; i++) {
        for (int j = 0; j < 500; j++) {
            temp[i + j] += a[j] * b[i];
        }
    }
    carry(temp);
    memcpy(res, temp, sizeof(temp));
}

void quick_power(int *a,int k) {
    res[0] = 1;
    while (k) {
        if (k & 1)
            multi(res, a, res);
        k >>= 1;
        if (k)
            multi(a, a, a);
    }
}

int main() {
    int p, k = 1;
    string s;
    cin >> s;
    for (int i = 0; i < s.length(); i++) {
        a[i] = s[s.length() - 1 - i] - '0';
    }
    cin >> p;
    quick_power(a, p);
    for (int i = 499; i >= 0; i--, k++) {
        cout << res[i];
        if (k % 50 == 0)cout << endl;
    }
    cout << endl;
    return 0;
}

矩陣快速冪

首先要了解：矩陣乘法

//Amn；Bnm
int** MatMulti(int **A, int **B, int m, int n) {
    int **C = new int*[n + 1];
    for (int i = 1; i <= n; i++) {
        C[i] = new int[n + 1];
        memset(C[i], 0, sizeof(int)*(n + 1));
    }
    
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= m; j++) {
            for (int k = 1; k <= n; k++) {
                if (A[i][k] == 0 || B[k][j] == 0)continue;
                C[i][j] += A[i][k] * B[k][j];
            }    
        }
    }
    return C;
}

上面只是簡單的計算矩陣的乘積，會感覺很抽象，因為上述矩陣並沒有具體的含義。

以常見的斐波那契數列為例：

F[n] = F[n-1] + F[n-2]. 由 F[0] = 0, F[1] = 1，可以遞推后面的所有數。

求第i項的復雜度是 O(n)；效率太低了。對於 n 的規模較大的題目無法在規定時間內求出。

把斐波那契數列的遞推式表示成矩陣就是：

 將記作矩陣A

於是有

 所以只要求出來A^n，就可以求出F_n了。

 題目鏈接：http://poj.org/problem?id=3070

題目給了一個更簡便的遞推式：

#include <iostream>
using namespace std;

void initialize(int m[][2]) {
    m[0][0] = 1; m[0][1] = 1; m[1][0] = 1; m[1][1] = 0;
}
void MatrixMulti(int m1[][2],int m2[][2],int m[][2]) {
    int t[2][2] = {0};
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            for (int k = 0; k < 2; k++) {
                t[i][j] = (t[i][j] + m1[i][k]%10000 * m2[k][j]%10000) % 10000;
            }
        }
    }
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            m[i][j] = t[i][j];
        }
    }
}
int MatrixQuickPower(int m[][2],int n) {
    int res[2][2] = { {1,0},{0,1} };
    while (n) {
        if (n & 1)
            MatrixMulti(res, m, res);
        MatrixMulti(m, m, m);
        n >>= 1;
    }
    return res[0][1];
}
int main() {
    int n, m[2][2];
    while (cin >> n) {
        if (n == -1)break;
        initialize(m);
        cout << MatrixQuickPower(m, n) << endl;
    }
    return 0;
}

矩陣快速冪是用來求解遞推式的，所以第一步先要列出遞推式:

 比如：F(n)=F(n-1)+F(n-2)

第二步是建立矩陣遞推式，找到轉移矩陣:

接下來就可以求解了

一些小性質

（1）計算 指數運算 結果的位數