彩笔
2022-12-21

快速幂

在acm被称为快速幂(平方求幂) 在ctf中也被称为 $montecarlo$

先说幂运算的朴素做法，在数学中，重复连乘的运算叫做乘方，乘方的结果称为 幂
${\displaystyle n}$ n 个相同的数 ${\displaystyle b}$ b 连续相乘

$$
b^n = \overbrace {b * b * b * ···} ^ {n \ number\ of\ b}
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// c++ version
#include <iostream>
int b, n, sum; // 仅作演示, 实际应该使用python或高精度防止爆int
int main() {
    std::cin >> b >> n;
    if (n == 0) { //对0次幂特判
        std::cout << 1;
        return 0;
    }
    sum = b, n = n - 1; 
    while (n --) {
        sum *= b;
    }
    for (int sum)
    std::cout << sum; return 0;
}

让我们先来思考一个问题：7的10次方，怎样算比较快？

方法1：上述提到的朴素算法, 来计算 $7 7=49$，$49 7=343$
这样算无疑太慢了，尤其对计算机的CPU而言，每次运算只乘上一个个位数，无疑太屈才了。这时我们想到，也许可以拆分问题。
方法2：先算7的5次方，即$7 7 7 7 7$，再算它的平方，共进行了5次乘法。
方法3：先算7 7得49，则7的5次方为$49 49 * 7$，再算它的平方，共进行了4次乘法。

递归快速幂

刚刚我们用到的，无非是一个二分的思路。我们很自然地可以得到一个递归方程：(虽然我觉得并不自然)
$$
b^n=\begin{cases}b^{n-1}\cdot b,&\text{if } n \text { is odd} \ b^{\frac{n}{2}}\cdot b^{\frac{n}{2}}, &\text{if } n \text { is even but not 0}\ 1,&\text{if } n=0\end{cases}
$$

计算a的n次方，如果n是偶数（不为0），那么就先计算a的n/2次方，然后平方；如果n是奇数，那么就先计算a的n-1次方，再乘上a；

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// c++ version
#include <iostream>
int b, n, temp, sum; // 依旧避免爆int
int powF(int b, int n) {
    if (n == 0) {
        return 1;
    } else if (n % 2 == 1) { // 奇数 odd
        return powF(b, n - 1) * b;
    } else { // 偶数 even
        temp = powF(b, n / 2);
        return temp * temp;
    }
}

int main() {
    std::cin >> b >> n;
    std::cout << powF(b, n);
    return 0;
}
注意，这个temp变量是必要的，因为如果不把 $a^{\frac{n}{2}}$ 记录下来，直接写成 $qpow(a, n /2)*qpow(a, n /2)$ ，那会计算两次 $a^{\frac{n}{2}}$ ，整个算法就退化为了 $O(n)$
在实际问题中，题目常常会要求对一个大素数取模，这是因为计算结果可能会非常巨大，但是在这里考察高精度又没有必要。这时我们的快速幂也应当进行取模，此时应当注意，原则是步步取模，如果 MOD 较大，还应当开long long
大家知道，递归虽然简洁，但会产生额外的空间开销。我们可以把递归改写为循环，来避免对栈空间的大量占用，也就是非递归快速幂

非递归快速幂

我们换一个角度来引入非递归的快速幂。还是7的10次方，但这次，我们把10写成二进制的形式，也就是 $(1010)_2$
现在我们要计算

$7^{(1010)_2}$ ，可以怎么做？我们很自然地想到可以把它拆分为 $7^{(1000)_2} \cdot 7^{(10)_2}$. 实际上，对于任意的整数，我们都可以把它拆成若干个 $7^{(100…)_2}$ 的形式相乘。而这些$7^{(100…)_2}$，恰好就是 $7^1$ 、$7^2$、$7^4$……我们只需不断把底数平方就可以算出它们。

我们先看代码，再来仔细推敲这个过程：

1
2
3
4
5
6
7
8
9
10
11
//非递归快速幂
int qpow(int a, int n){
    int ans = 1;
    while(n){
        if(n&1)        //如果n的当前末位为1
            ans *= a;  //ans乘上当前的a
        a *= a;        //a自乘
        n >>= 1;       //n往右移一位
    }
    return ans;
}

最初ans为1，然后我们一位一位算：

1010的最后一位是0，所以a^1这一位不要。然后1010变为101，a变为a^2。

101的最后一位是1，所以a^2这一位是需要的，乘入ans。101变为10，a再自乘。

10的最后一位是0，跳过，右移，自乘。

然后1的最后一位是1，ans再乘上a^8。循环结束，返回结果。

这里的位运算符，>>是右移，表示把二进制数往右移一位，相当于/2；&是按位与，&1可以理解为取出二进制数的最后一位，相当于 $\mod2==1$。这么一等价，是不是看出了递归和非递归的快速幂的关系了？虽然非递归快速幂因为牵扯到二进制理解起来稍微复杂一点，但基本思路其实和递归快速幂没有太大的出入。

快速幂的拓展

上面所述的都是整数的快速幂，但其实，在算 a^n 时，只要a的数据类型支持乘法且满足结合律，快速幂的算法都是有效的。矩阵、高精度整数，都可以照搬这个思路。下面给出一个模板：

注意，较复杂类型的快速幂的时间复杂度不再是简单的 O(\log n) ，它与底数的乘法的时间复杂度有关。

例如，矩阵快速幂的一个经典应用是求斐波那契数列：

（洛谷P1962） 斐波那契数列

题目背景

大家都知道，斐波那契数列是满足如下性质的一个数列：
$$
F_n = \begin{cases}1& (n \le 2) \ F_{n-1}+F_{n-2}& (n\ge 3) \end{cases}
$$
题目描述

请你求出$F_n \bmod 10^9 + 7$的值。

（以下内容涉及到基本的线性代数知识）

设矩阵

$$
A=\begin{bmatrix}0 &1\ 1 & 1\end{bmatrix}
$$ ，我们有

$$
A\begin{bmatrix}F_n\ F_{n+1}\end{bmatrix} = \begin{bmatrix}F_{n+1}\ F_n+F_{n+1}\end{bmatrix}=\begin{bmatrix}F_{n+1}\ F_{n+2}\end{bmatrix}
$$

，于是 :

$$
\begin{aligned} \begin{bmatrix}F_n\ F_{n+1}\end{bmatrix} &= A\begin{bmatrix}F_{n-1}\ F_n\end{bmatrix}\&=A^2\begin{bmatrix}F_{n-2}\ F_{n-1}\end{bmatrix}\&=…\&=A^{n-1}\begin{bmatrix}F_1\ F_2\end{bmatrix}\&=A^{n-1}\begin{bmatrix}1\ 1\end{bmatrix} \end{aligned}
$$

这样，我们把原来较为复杂的问题转化成了求某个矩阵的幂的问题，这就可以应用快速幂求解了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include <cstdio>
#define MOD 1000000007
typedef long long ll;

struct matrix
{
    ll a1, a2, b1, b2;
    matrix(ll a1, ll a2, ll b1, ll b2) : a1(a1), a2(a2), b1(b1), b2(b2) {}
    matrix operator*(const matrix &y)
    {
        matrix ans((a1 * y.a1 + a2 * y.b1) % MOD,
                   (a1 * y.a2 + a2 * y.b2) % MOD,
                   (b1 * y.a1 + b2 * y.b1) % MOD,
                   (b1 * y.a2 + b2 * y.b2) % MOD);
        return ans;
    }
};

matrix qpow(matrix a, ll n)
{
    matrix ans(1, 0, 0, 1); //单位矩阵
    while (n)
    {
        if (n & 1)
            ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

int main()
{
    ll x;
    matrix M(0, 1, 1, 1);
    scanf("%lld", &x);
    matrix ans = qpow(M, x - 1);
    printf("%lld\n", (ans.a1 + ans.a2) % MOD);
    return 0;
}

![[Pasted image 20221221012506.png]]

引用几乎是复制了这篇文章: https://zhuanlan.zhihu.com/p/95902286

彩笔

彩笔

熟练拼写:XCPC,Crypto