快速幂与矩阵快速幂

好久没有学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之前的状态吧。

posted @ 2021-05-20 16:44  ZTer  阅读(343)  评论(0编辑  收藏  举报