基礎算法—快速冪詳解

冪運算是非常常見的一種運算，求取$a^n$,最容易想到的方法便是通過循環逐個累乘，其復雜度為$O(n)$，這在很多時候是不夠快的，所以我們需要一種算法來優化冪運算的過程。

一、快速冪——反復平方法

該怎樣去加速冪運算的過程呢？既然我們覺得將冪運算分為n步進行太慢，那我們就要想辦法減少步驟，把其中的某一部分合成一步來進行。

比如，如果$n$能被2整除，那我們可以先計算一半，得到$a^{n/2}$的值，再把這個值平方得出結果。這樣做雖然有優化，但優化的程度很小，仍是線性的復雜度。

再比如，如果我們能找到$2^k = n$，那我們就能把原來的運算優化成$((a^2)^2)^2...$，只需要$k$次運算就可以完成，效率大大提升。可惜的是，這種條件顯然太苛刻了，適用范圍很小。不過這給了我們一種思路，雖然我們很難找到$2^k = n$，但我們能夠找到$2^{k_1} + 2^{k_2} + 2^{k_3} +......+ 2^{k_m} = n$。這樣，我們可以通過遞推，在很短的時間內求出各個項的值。

我們都學習過進制與進制的轉換，知道一個$b$進制數的值可以表示為各個數位的值與權值之積的總和。比如，2進制數$1001$，它的值可以表示為10進制的$1\times2^3 + 0\times2^2 + 0\times2^1 + 1\times2^0$，即$9$。這完美地符合了上面的要求。可以通過2進制來把$n$轉化成$2^{k_m}$的序列之和，而2進制中第$i$位（從右邊開始計數，值為$1$或是$0$）則標記了對應的$2^{i - 1}$是否存在於序列之中。譬如，$13$為二進制的$1101$，他可以表示為$2^3 + 2^2 + 2^0$，其中由於第二位為$0$，$2^1$項被舍去。

如此一來，我們只需要計算$a、a^2、a^4、a^8......a^{2^{k_m}}$的值（這個序列中的項不一定都存在，由$n$的二進制決定）並把它們乘起來即可完成整個冪運算。借助位運算的操作，可以很方便地實現這一算法，其復雜度為$O(\log n)$。

typedef long long ll;
ll mod;
ll qpow(ll a, ll n)//計算a^n % mod
{
    ll re = 1;
    while(n)
    {
        if(n & 1)//判斷n的最后一位是否為1
            re = (re * a) % mod;
        n >>= 1;//舍去n的最后一位
        a = (a * a) % mod;//將a平方
    }
    return re % mod;
}

取模運算一般情況下是需要的，當然也可以省去。

二、矩陣快速冪

快速冪只是通過二進制拆分$n$來加速冪運算的手段，當然並不只適用於求取數字的冪次，對於矩陣的$n$次方，也可以用同樣的手段求取。除了乘法的規則與上面的快速冪不同之外，其他方面並沒有太大的差別。

不過這有什么意義呢？利用矩陣的冪次，我們可以快速地完成遞推。

比如，在POJ3070 Fibonacci中，就需要我們快速地求取斐波那契數列的第$n$項（取模），對於$n=10^{18}$，一步步推過去顯然太慢了，那么我們可以考慮構造矩陣來幫我們完成遞推。

首先復習一下矩陣的乘法:

$$\left[\begin{array}{cccc}{a_{11}} & {a_{12}} & {\dots} & {a_{1 n}} \\ {a_{21}} & {a_{22}} & {\dots} & {a_{2 n}} \\ {\dots} & {\dots} & {\dots} & {\dots} \\ {a_{n 1}} & {a_{n 2}} & {\dots} & {a_{n n}}\end{array}\right] \left[\begin{array}{cccc}{b_{11}} & {b_{12}} & {\dots} & {b_{1 n}} \\ {b_{21}} & {b_{22}} & {\dots} & {b_{2 n}} \\ {\dots} & {\dots} & {\dots} & {\dots} \\ {b_{n 1}} & {b_{n 2}} & {\dots} & {b_{n n}}\end{array}\right]=\left[\begin{array}{cccc}{c_{11}} & {c_{12}} & {\dots} & {c_{1 n}} \\ {c_{21}} & {c_{22}} & {\dots} & {c_{2 n}} \\ {\dots} & {\dots} & {\dots} & {\dots} \\ {c_{n 1}} & {c_{n 2}} & {\dots} & {c_{n n}}\end{array}\right]$$

 

其中$C_{i j}=\sum_{k=1}^{n} a_{i k} * b_{k j}$

對於$i \ge 3$

有$Fib_i = Fib_{i-1} + Fib_{i-2}$

可以構造出矩陣遞推式$$\left[\begin{array}{ll}{1} & {1} \\ {1} & {0}\end{array}\right]\left[\begin{array}{c}{F i b_{i}} \\ {F i b_{i-1}}\end{array}\right]=\left[\begin{array}{c}{F i b_{i+1}} \\ {F i b_{i}}\end{array}\right]$$

那么$$\left[\begin{array}{ll}{1} & {1} \\ {1} & {0}\end{array}\right]^{n-2}\left[\begin{array}{l}{1} \\ {1}\end{array}\right]=\left[\begin{array}{c}{F i b_{n}} \\ {F i b_{n-1}}\end{array}\right]$$

我們就可以利用矩陣快速冪以$O(\log n)$求取斐波那契數列的第$n$項了。

const ll mod = 10000;
const int maxv = 2;

struct Matrix {
    ll a[maxv][maxv]; //矩陣

    Matrix operator*(const Matrix &b) const& { 
        //矩陣乘法，復雜度O(maxv^3)，也可看作常數，但maxv較大（大於5）時會使運算時間提高好幾個數量級
        Matrix ans;
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                ans.a[i][j] = 0;
                for (int k = 0; k < maxv; ++k) {
                    ans.a[i][j] += a[i][k] * b.a[k][j] % mod;
                    ans.a[i][j] %= mod;
                }
            }
        }
        return ans;
    }
    static Matrix qpow(Matrix x, ll n) {//矩陣快速冪，將乘法復雜度看作常數則復雜度為O(log n)
        Matrix ans;
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                if (i == j)
                    ans.a[i][j] = 1;
                else
                    ans.a[i][j] = 0;
            }
        }//初始化為單位矩陣，參考普通數字的快速冪這里初始化為1
        while (n) {//其余細節基本相同
            if (n & 1)
                ans = ans * x;
            x = x * x;
            n >>= 1;
        }
        return ans;
    }
    Matrix(ll temp[maxv][maxv]) {//構造方法
        for (int i = 0; i < maxv; ++i) {
            for (int j = 0; j < maxv; ++j) {
                a[i][j] = temp[i][j];
            }
        }
    }
    Matrix() { }
};

ll fib(ll n) {//求取斐波那契數列第n項（本題取模）
    if (n == 0)
        return 0;
    if (n <= 2)
        return 1;
    
    ll temp[maxv][maxv] = {
        1, 1,
        1, 0
    };
    Matrix m(temp);
    m = Matrix::qpow(m, n - 2);
    return (m.a[0][1] + m.a[0][0]) % mod;
}

 當然了，構造矩陣的方法是多種多樣的，例如上面的題目中就給出了另一種構造出斐波那契數列的方法。

矩陣快速冪能解決的遞推也遠不止斐波那契數列，下面根據我個人的習慣列舉幾個比較常見、比較“套路”的情況：

（1）對於$a_n = x\cdot a_{n-1} + y\cdot a_{n-2} + c^n$，可以得到

$$
\left[\begin{array}{ccc}{x} & {y} & {c} \\ {1} & {0} & {0} \\ {0} & {0} & {c}\end{array}\right]\left[\begin{array}{c}{a_{i}} \\ {a_{i-1}} \\ {c^{i}}\end{array}\right]=\left[\begin{array}{c}{a_{i+1}} \\ {a_{i}} \\ {c^{i+1}}\end{array}\right]
$$

（2）對於$a_n = a_{n-1} + a_{n-2} + n^2$，可以得到

$$
\left[\begin{array}{ccccc}{1} & {1} & {1} & {2} & {1} \\ {1} & {0} & {0} & {0} & {0} \\ {0} & {0} & {1} & {2} & {1} \\ {0} & {0} & {0} & {1} & {1} \\ {0} & {0} & {0} & {0} & {1}\end{array}\right]\left[\begin{array}{c}{a_{i}} \\ {a_{i-1}} \\ {i^{2}} \\ {i} \\ {1}\end{array}\right]=\left[\begin{array}{c}{a_{i+1}} \\ {a_{i}} \\ {(i+1)^{2}} \\ {i+1} \\ {1}\end{array}\right]
$$

（3）對於$a_n=x\cdot a_{n-1} + y\cdot a_{n-2}$，求$a_n$的前綴和$s_n$

$$
\left[\begin{array}{lll}{1} & {x} & {y} \\ {0} & {x} & {y} \\ {0} & {1} & {0}\end{array}\right]\left[\begin{array}{c}{s_{i}} \\ {a_{i}} \\ {a_{i-1}}\end{array}\right]=\left[\begin{array}{c}{s_{i+1}} \\ {a_{i+1}} \\ {a_{i}}\end{array}\right]
$$