Hexo
0%

4-快速幂-矩阵快速幂

发表于 更新于

快速幂-矩阵快速幂

快速幂就是快速算底数的n次幂。其时间复杂度为 O(log₂N)， 与朴素的O(N)相比效率有了极大的提高

(0)同余定理

(a×b)%c=(a%c×b%c)%c

(a+b)%c=(a%c + b%c)%c

(a-b)%c=(a%c - b%c + c)%c

快速幂-矩阵快速幂

(1)快速幂

  1. 法一

$a^{2k}=(a^{k})^{2}$ 幂次为偶数情况下

$a^{2k+1}=(a^{k})^{2}*a$ 幂次为奇数情况下

代码以递归形式实现$a^{n}$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
long long mod=1e9+7;
long long quick_pow1(long long a,long long n)
{
    //a^n
    if (n<0)
    {
        cerr << "error quick pow" << endl;
        return -1;
    }
    if (a==0)
        return 0;
    if (n==0)
        return 1;
    a%=mod;
    if ((n&1)==1)
        return a*quick_pow1(a,n-1)%mod;
    else
    {
        long long half=quick_pow1(a,n/2);
        return half*half%mod;
    }
}

  1. 法二

求$a^{k}$则将幂次k按照2的幂次展开，即方法：先将k转换为二进制数，把二进制数首先写成加权系数展开式，计算展开式中的每个加数，即展开成：以a为底，幂次为2的n次幂的数，相乘，前提是k>=0

例子：求，则29D=11101B=+++=1+4+8+16=29D
即= *

代码以迭代形式实现

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
long long mod=1e9+7;
long long quick_pow2(long long a,long long n)
{
    if (n<0)
    {
        cerr << "error quick pow" << endl;
        return -1;
    }
    if (a==0)
        return 0;
    a%=mod;
    long long ans=1;
    while (n!=0)
    {
        if ((n&1)==1)
        {
            ans*=a;
            ans%=mod;
        }
        a*=a;
        a%=mod;
        n>>=1;
    }
    return ans;
}

(3)矩阵快速幂

例子：求斐波那契数列

常规解法是直接dp模拟一遍，时间复杂度O(n)，空间复杂度用滚动数组可做到O(1)

解：显然 =
这里构建了一个矩阵A，使得 A =
又可以推导得 =
即最终公式1: A = ，A =
且矩阵无交换律，但是有结合律，即A = 且 A =
可以推导得 A = A (A ) = =
即最终通项公式2: * = ， 这里斐波那契额数列A = , 将要求的f(k)代入公式中的f(n)求即可

所以转化成新问题：求，此时A是一个矩阵。此时只需结合快速幂的方法二即可。

tips：为了设计方便，通常将通项公式中的所有矩阵都拓展到A的维度
即此题求斐波那契数列，通项公式优化为 * =

求斐波那契数列，用矩阵快速幂方法

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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
/*
* f(0)=f(1)=1
* f(i)=f(i-1)+f(i-2) 求f(n)，通过矩阵快速幂
* 则矩阵通项公式
A=[0 1]   即 A^n * [f(0)] = [f(n)  ]
  [1 1]            [f(1)] = [f(n+1)]
*/
//矩阵A的维度是M*M，这里是2*2
const int M=2;
//矩阵A的类
class Ma
{
public:
    int a[M][M];
    Ma()
    {
        memset(a,0,sizeof(a));
    }
    void become_unit_matrix()
    {
        //设置成单位矩阵
        for (int i=0;i<M;++i)
            for (int j=0;j<M;++j)
                a[i][j]= i==j ? 1 : 0;
    }
    void show()
    {
        cout << "----" << endl;
        cout << "Matrix is " << M << '*' << M << endl;
        for (int i=0;i<M;++i)
        {
            for (int j=0;j<M;++j)
                cout << a[i][j] << ' ';
            cout << endl;
        }
        cout << "----" << endl;
    }
    Ma operator*(const Ma &B) const
    {
        //重载两个矩阵相乘
        Ma ans;
        for (int i=0;i<M;++i)
            for (int j=0;j<M;++j)
                for (int k=0;k<M;++k)
                    ans.a[i][j]+=(this->a[i][k])*(B.a[k][j]);
        return ans;
    }
    Ma operator^(int n) const
    {
        //重载幂次运算，并结合快速幂
        Ma ans;
        //单位矩阵乘任何矩阵还是任何矩阵本身
        ans.become_unit_matrix();
        Ma A=*this;
        while(n!=0)
        {
            if ((n&1)==1)
                ans=ans*A;
            A=A*A;
            n>>=1;
        }
        return ans;
    }
};

int get_fibo(int n)
{
    // 0 1 2 3 4 5  6  7  8
    // 1 1 2 3 5 8 13 21 34
    Ma A;
    A.a[0][0]=0;
    A.a[0][1]=A.a[1][0]=A.a[1][1]=1;
    //A.show();
    Ma F;
    F.a[0][0]=F.a[1][0]=1;
    //F.show();
    Ma ans=(A^n)*F;
    //ans.show();
    return ans.a[0][0];
}