• 快速幂详解(幂运算与矩阵)


    基础知识

    矩阵乘法

    一张图说明足矣:

    代码实现(C)

    const int N = 100;
    int c[N][N];
    void multi(int a[][N], int b[][N], int n) {
        memset(c, 0, sizeof c);
        for(int i = 1; i <= n; i++) {
            for(int j = 1; j <= n; j++) {
                for(int k = 1; k <= n; k++)
                    c[i][j] += a[i][k] * b[k][j];
            }
        }
        return;
    }        

    另一种实现,将二层循环与三层循环对调,从而可以在第二层for循环时先判断 if(a[i][k] == 0) continue; ,对于矩阵有较多0的有一定效果。

    int c[N][N];
    void multi(int a[][N], int b[][N], int n) {
        memset(c, 0, sizeof c);
        for(int i = 1;i <= n; i++) {
            for(int k = 1; k <= n; k++) {
                if(a[i][k] == 0) continue; 
                for(int j = 1; j <= n; j++) {
                    c[i][j] += a[i][k] * b[k][j];
                }
            }
        }
        return;
    }   

    显然矩阵乘法的复杂度是$O(n^3)$

    矩阵乘法需要A的行数与B的列数相等的(这是A*B的前提条件,可见矩阵乘法是不满足交换律)。

    矩阵快速幂只会用到方阵(除非题目是裸的矩阵乘法),矩阵快速幂都是方阵也就避免了检验相乘的前提条件,可以放心用。

    快速幂

    经典例题:数A的幂

    我们利用分治的算法思想可以考虑如下求解一个数A的幂:

    例如:X的19次方。

    19的二进制为1 0 0 1 1。由$X^m * X^n = X^{m+n}$,则$X^{19} = X^{16} * X^2 * X^1$。

    代码实现(C++)

    double QuickPow(double x, int N) {
        if (N == 0) return 1;
        if (x == 0) return 0;
        if (N < 0) {
        x = 1 / x;
            N = -N;        
        }
        double ans = 1.0;
        while(N) {
            if(N & 1) {
                ans *= x;
            }
            x *= x;
            N >>= 1;
        }
        return ans;
    }

    过程解析

    对于$X^19$来说,初始: ans = 1; res = x; 则

    10011最后一位是1,所以是奇数:

    ans = res * ans = x; 
    res = res * res = x^2; 

    然后右移一位,1 0 0 1,则1001最后一位是1,所以是奇数:

    ans = res * ans = x * (x^2) = x^3;   
    res = res * res = x^2 * x^2 = x^4;

    然后右移一位,1 0 0,则最后一位是0,所以当前的数为偶数:

    res = res * res = x^4 * x^4 = x^8;

    然后右移一位,1 0,最后一位是0,当前数是偶数:

    res = res * res = x^8 * x^8 = x^16;

    然后右移一位,1,最后一位是1,当前数是奇数:

    ans = ans * res = (x^3) *(x^16) = x^19;
    res = res * res = x^32;

    可以看出res = X^m,m始终是与二进制位置上的权值是相对应的。当二进制位为0时,我们只让res幂指数为2,对应下一个二进制位的权值;当二进制位为1时, ans = ans * res; ,则乘上了该乘的X幂次。

    矩阵快速幂

    上面介绍了整型快速幂的实现思路,而矩阵乘法可以用来优化很多线性递推式,那么将其中的乘法换成矩阵乘法,就变成了矩阵快速幂

    代码实现(C)

    struct Matrix { //定义一个矩阵类型的结构体
        int m[maxN][maxN];
    } ans, res;
    
    //实现矩阵乘法A * B, n为方阵的阶
    Matrix Mul(Matrix A, Matrix B, int n) { 
        Matrix tmp;
        memset(tmp.m, 0, sizeof(tmp.m));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < n; k++)
                    tmp.m[i][j] = A.m[i][k] + B.m[k][j];
            }
        }
        return tmp;
    }
    
    //快速幂算法,求矩阵res的N次幂
    void QuickPower(int N, int n) {
        //整数快速幂ans初始化为1,对于矩阵乘法来说,则ans初始化为单位矩阵
        //有A * E = A
        memset(ans.m, 0, sizeof(ans.m));
        for (int i = 0; i < n; i++) ans.m[i][i] = 1;
        while (N) {
            if (N & 1) ans = Mul(ans, res);
            res = Mul(res, res);
            N = N >> 1;
        }
    }

    经典例题:斐波拉契数列

    众所周知:斐波那契数列递推公式为:$F[n] = F[n-1] + F[n-2]$, 由$F[0] = 0$,$F[1] = 1$,可以递推后面的所有数。

    使用for循环是最简单直接的解法,然而对于像POJ 3070这样的题目,求斐波那契数列,其n可高达10亿,直接递推有其局限性:

    (1)当n过大,递推次数过多,而测试时间有限时,可能会导致超时。

    (2)想要打表实现随机访问根本不可能,先把斐波那契数列求到10亿,然后想去进行随机访问。题目未给出那么多内存,数组也开不到10亿。

    因此另一个更快的解法是矩阵快速幂。

    建立矩阵递推式,找到转移矩阵

    显然,A就是转移矩阵。

    构建矩阵递推的思路

    $$A * F(n-1) = F(n)$$

    一般$F(n)$与$F(n-1)$都是按照原始递推式来构建的,当然可以先猜一个$F(n)$,主要是利用矩阵乘法凑出矩阵$A$,第一行一般就是递推式,后面的行就是不需要的项就让与其相乘的系数为0。矩阵$A$就叫做转移矩阵(一定是常数矩阵),它能把$F(n-1)$转移到$F(n)$,然后这就是个等比数列,直接写出通项:$F(n) = A^{n-1} * F(1)$,此处$F1$叫初始矩阵。所以用一下矩阵快速幂然后乘上初始矩阵就能得到$F(n)$。例如,这里的$F(n)$就两个元素(两个位置),根据自己设置的$F(n)$对应位置就是对应的值。

    示例

    1. $f(n) = a * f(n-1) + b * f(n-2) + c$(a, b, c是常数)

    $$left( egin{array}{ccc} a & b & 1 \ 1 & 0 & 0 \ 1 & 1 & 1\ end{array} ight) * left( egin{array}{c} f_{n-1} \ f_{n-2} \ c \ end{array} ight) = left( egin{array}{c} f_n \ f_{n-1} \ c \ end{array} ight)$$

    2. $f(n) = c^n - f(n-1) $(c是常数)

    $$left( egin{array}{cc} -1 & c \ 0 & c \ end{array} ight) * left( egin{array}{c} f_{n-1} \ c^{n-1} \ end{array} ight) = left( egin{array}{c} f_n \ c^n \ end{array} ight)$$

    (整理自网络)

    参考资料:

    https://www.cnblogs.com/cmmdc/p/6936196.html

    https://blog.csdn.net/wust_zzwh/article/details/52058209

    class Solution {public:    double myPow(double x, int n) {        if (n == 0) return 1;        if(x==0)  return 0;        long N = n;        if (n < 0) {            x = 1/x;            N = -N;        }        double res = 1;       //快速幂        while (N > 0) {            if ((N & 1) == 1) {                res = res * x;            }            x *= x;            N >>= 1;        }        return res;    }};

    Min是清明的茗
  • 相关阅读:
    记录美好生活:
    _2data=data.find_all(class_='cvesummarylong')#获取在srrowns里面没有的那个数据
    _1data=data.find_all(class_='srrowns')#获取所有以srrowns为标签的数据
    idea is good
    创业基础(第六章:创业资源及其管理) 来自高校:全国大学生创新创业实践联盟 分类:创新创业 学习规则:按序学习
    创业基础(第四章: 创业风险及识别与管理) 来自高校:全国大学生创新创业实践联盟 分类:创新创业 学习规则:按序学习
    hdu6447 YJJ's Salesman
    hdu6438 Buy and Resell
    论开学第二个月干了点啥
    论开学第一个月干了点啥
  • 原文地址:https://www.cnblogs.com/MinPage/p/14224286.html
Copyright © 2020-2023  润新知