快速幂与矩阵快速幂

Posted 2021-05-20 欢迎光临一铭君一的家！

好久没有学OI了，今天突然回来打了一场比赛，T1是矩阵快速幂+高斯消元，想到快速幂算法啥的我也一直没有整理笔记，索性一起整理一下。

Part 1 什么是快速幂？

假设现在我们有两个整数：\\(n,k\\)现在要求你求出\\(n^k\\)对\\(1e9+7\\)的模是多少？

显然一个简单的算法就是\\(O(n)\\)把\\(k\\)个\\(n\\)乘起来，边乘边取模。

但是当\\(k\\)非常大，比如\\(k=1e9\\)的时候，这个算法的复杂度并不优秀。

那么怎么办呢？这个时候就要用到快速幂算法了。

另外，对于矩阵也有乘法运算，这里不再展开讲。矩阵快速幂算法就是对一个\\(n*n\\)的矩阵\\(M\\)，快速求\\(M^k\\)的算法（具体为什么是\\(n*n\\)涉及到矩阵乘法的性质，之后会提到）。

Part 2 快速幂、矩阵快速幂算法的原理

这里我们复习一下初中数学。

（1）如果将\\(n\\)自乘一次，就会变成 \\(n^2\\)，再把\\(n^2\\)自乘一次就会变成\\(n^4\\)，然后是\\(n^8\\)。以此类推，自乘\\(m\\)次的结果是\\(n^{2^{m}}\\) 。

（2）\\(a^x·a^y = a^{x+y}\\)

（3）将\\(k\\)转化为二进制观看一下：比如\\(k=11=(1011)_{2}\\)。从左到右，这些\\(1\\)分别代表十进制的\\(8,2,1\\)，可以说\\(n^{11} = n^8·n^2·n^1\\)

为什么要这样表示？因为在快速幂的过程中，我们会把\\(a\\)自乘为\\(n^2\\)，然后\\(n^2\\)自乘为\\(n^4\\)……像上面第一条说的。

相应的，如果我们的乘数\\(k\\)的二进制表示在右数第\\(x\\)位上为\\(1\\)的话，就代表\\(n^k\\)可以分解出一个因数：\\(n^{2^x}\\)，反之，在右数第\\(x\\)位上为\\(0\\)的话，就代表\\(n^k\\)不可以分解出上述因数。

所以，在求解时，我们可以创建一个\\(base\\)，初始等于\\(n\\)。然后扫描\\(k\\)的二进制表示，每扫描一位，\\(base\\)自我平方，如果\\(k\\)的这一位为\\(1\\)，那么把\\(n\\)乘上\\(base\\)（\\(base\\)相当于记录了\\(n^{2^x}\\)的值，若\\(n^k\\)分解出了这个值，那么答案\\(n\\)乘上\\(base\\)），扫描到最后一位的时候，就求出了\\(n^k\\)的值啦！

矩阵快速幂的原理同上，只不过把普通的乘法换成了矩阵乘法。

另外，关于取模运算，一般的，有公式：

\\[(a+b)\\%c=(a\\%c+b\\%c)\\%c \\]

\\[(a·b)\\%c=(a\\%c·b\\%c)\\%c \\]

这个应该都知道吧

Part 3 代码实现快速幂

显然快速幂的核心就是二进制分解\\(n^k\\)。

与运算符 a&b 的作用：

假设 a 的二进制表示为 1010110 ，b 的二进制表示是 1001101 ，那么 a&b 的结果就是 1000100 。简单来说，就是当 a、b 的某一位同时取 1 的时候，结果才为 1 ，否则为 0 。

显然，k&1的结果可以表示\\(k\\)的二进制的第一位（1 是 000...1）。

右移运算符 a>>x 的作用：

把运算符左侧的数的二进制表示右移 x 位（相当于从右边开始去掉了 x 位，其余不变）。

在\\(k\\)检测完第一位的时候，直接k>>=1去掉\\(k\\)的第一位就好。

代码实现：

#define ll long long
ll Mod=998244353;
inline ll FastPow(ll a,ll k){//a是底数，k是指数
    if(k==0)//特判一下
        return a%Mod;
    ll ans=1;//答案，这里我把读入进函数的a当作上文中提到的base
    while(k){
        if(k&1)//最后一位是1
            ans=(ans%Mod*a%Mod)%Mod;//答案乘以base
        a=(a%Mod*a%Mod)%Mod;//base自平方
        k>>=1;//去掉二进制第一位
    }
    return ans;//返回答案
}

Part 4 关于矩阵乘法

想要做矩阵快速幂，先要了解矩阵乘法。（考虑到我数学水平有限，这里不展开讲）

先定义矩阵\\(A_{n\\times k}\\)、\\(B_{k\\times m}\\)。

\\(A= \\left\\{ \\begin{matrix} a_{11} & a_{12} & \\cdots & a_{1k}\\\\ a_{21} & a_{22} & \\cdots & a_{2k} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ a_{n1} & a_{n2} & \\cdots & a_{nk} \\end{matrix} \\right\\}\\)，\\(B= \\left\\{ \\begin{matrix} b_{11} & b_{12} & \\cdots & b_{1m}\\\\ b_{21} & b_{22} & \\cdots & b_{2m} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ b_{k1} & b_{k2} & \\cdots & b_{km} \\end{matrix} \\right\\}\\)

矩阵乘法的定义如下：

（1）\\(A\\)（左矩阵）的列数必须与 \\(B\\)（右矩阵）的行数相同，这也是矩阵乘方要求是正方形矩阵的原因，即矩阵自乘要求行和列数相等。

（2） \\(C\\) 的行数与 \\(A\\)（左矩阵）相同，列数与\\(B\\)（右矩阵）相同，即\\(C_{n\\times m}\\)

（3） \\(C\\) 的第 \\(i\\) 行第 \\(j\\) 列的元素 \\(C_{ij}\\) 由 \\(A\\) 的第 \\(i\\) 行元素与 \\(B\\) 的第 \\(j\\) 列元素对应相乘，再取乘积之和。

用一个公式来说就是：

\\[C_{ij}=\\sum_{i=1}^{k}a_{ik}·b_{kj} \\]

关于单位矩阵：

单位矩阵是一个左上到右下对角线的元素都为\\(1\\)，其余元素为\\(0\\)的矩阵，一个矩阵乘以单位矩阵（在行列数允许的情况下），根据矩阵乘法的定义，得到的结果是这个矩阵本身。

OK，了解到这里已经足够我们做矩阵快速幂了。

Part 5 矩阵快速幂代码实现

对指数\\(k\\)的操作同上，这里不再复述。

重要的是，重载乘法运算符，使得我们在上面的式子对于矩阵也适用。

代码实现：

const ll Mod=1e9+7;
ll n,k;
ll M[maxn][maxn];//这里的M是初始矩阵

struct matrix{
    ll M[maxn][maxn];//定义一个结构体
}A,res;

matrix operator * (matrix &a,matrix &b){//重载这个结构体的*运算符
    matrix c;//c是运算结果
    memset(c.M,0,sizeof(c.M));
    for(int i=0;i<n;i++)//改为矩阵乘法的原理
        for(int j=0;j<n;j++)
            for(int k=0;k<n;k++)
                c.M[i][j]=(c.M[i][j]%Mod+((a.M[i][k]%Mod)*(b.M[k][j]%Mod))%Mod)%Mod;//乘法取模
    return c;//返回结果矩阵c
}

inline void FastMatrixPow(matrix a,ll k){
    memset(res.M,0,sizeof(res.M));
    for(int i=0;i<n;i++)
        res.M[i][i]=1;//创建单位矩阵（相当于1）作为初始答案
    while(k){
        if(k&1)
            res=res*a;
        a=a*a;//同上操作，a作为base自乘
        k>>=1;
    }
}

Part 6 结语

其实快速幂算法也是分治算法的一种，体现的是分而治之的思想。核心就是把指数分解成几个2的幂的乘积的形式，然后分别求解。

今天是重回OI的第5天，希望可以快速恢复到去年CSP之前的状态吧。