从快速幂到dp 优化:矩阵快速幂

幂运算

na 相加我们当然不会写成一个循环,na 相乘我们自然要用幂运算。

幂运算裸题

题目链接

PAT L1-012: 计算指数

题意

输入数字n,求出2n 的值,n[1,10]

解法

用上cmath 头文件中的pow() 函数即可,由于本题数据范围极小,所以根本不会出现什么精度问题。

过题代码

#include <iostream>
#include <cmath>
using namespace std;

int main() {
    int n;
    cin >> n;
    cout << 2 << "^" << n << " = " << pow(2, n) << endl;
    return 0;
}

快速幂运算

这只是相对于比较小的数据而言,但是如果要计算的数据范围较大,大到连double 都存不下时,用pow() 函数显然会输出错误的结果,而如果写一个循环,显然不可能在规定时间内得到结果,如果题目没有要求结果对某个数取模,算一算结果的位数,大概就要用到java 的BigInteger 类,不过以下仅讨论结果对某个数取模时的写法,即使要用到BigInteger 类,也是相同的写法,只是换个语言罢了。
先在这里提以下两个取模公式,详细请点传送门

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

证明略。

算法思想

对于幂运算我们知道,

abc=(ab)cab+c=ab×ac

有了上面两个公式,我们就可以将幂运算进行下面的转换:
xn={(xn12)2×x(n%2=1)(xn2)2(n%2=0)

而右边小括号里的xn12xn2 又是一个幂运算,于是我们就可以递归转换下去,直到n=0。现在来看一下这样转换之后的计算次数有多少次:
如果是要计算xn,就要先计算xn2xn12,要计算这个值,就要先计算xn4xn14 或……以这样的方式递减下去,我们可以看到,将在log2n 的次数内下降至x0,所以要计算xn 的时间复杂度就为O(logn)

幂运算题改

题意

n=18654203254124875 时求出2n 的结果,由于该结果非常大,所以只要求输出mod(2n,109+7) 的值。

题解:快速幂

代码

#include <iostream>
using namespace std;

#define LL long long

const LL MOD = 1000000007;

LL mi(LL res, LL n) {
    LL ans;
    for(ans = 1; n != 0; n >>= 1) {
        if(n % 2 == 1) {
            ans = (ans * res) % MOD;
        }
        res = (res * res) % MOD;
//        cout << "ans = " << res << "^" << n / 2 << " * " << ans << endl << endl;
//        cout << "ans = " << ans << "    res = " << res << "    n = " << n << endl << endl;
    }
    return ans;
}

int main() {
    LL n;
    cin >> n;
    cout << "2^" << n << " = " << mi(2, n) << endl;

    return 0;
}

对代码快速幂函数部分有疑问的,可以去掉循环中两个cout 语句的注释,来看看每次循环的值的变化是什么样的。

快速幂裸题

题目链接

POJ 1995: Raising Modulo Numbers

题意

给定HAB 的值以及M,求出

(A1B1+A2B2+A3B3+...+AHBH)modM

的值,其中M,H[1,45000]AB 不同时为0。

题解:快速幂

过题代码

#include <iostream>
using namespace std;

#define LL long long

LL mi(LL res, LL n, LL MOD) {
    LL ans;
    for(ans = 1; n != 0; n >>= 1) {
        if(n % 2 == 1) {
            ans = (ans * res) % MOD;
        }
        res = (res * res) % MOD;
    }
    return ans;
}

int main() {
    ios::sync_with_stdio(false);
    LL Z, M, N;
    LL ans, a, b;

    cin >> Z;
    while(Z--) {
        ans = 0;
        cin >> M >> N;
        for(int i = 0; i < N; ++i) {
            cin >> a >> b;
            ans = (ans + mi(a, b, M)) % M;
        }
        cout << ans << endl;
    }

    return 0;
}

简单dp

关于动态规划

dp 是动态规划的简称,动态规划大家第一次听大概是在上课的时候,老师教到递归求斐波那契数列时提到的动态规划记录状态的解法吧。动态规划什么原理什么性的,很多博客已经说得很清楚了,大多是复制粘贴,这里不再赘述,想要详细了解的,请点传送门
如果觉得传送门里的例子太难理解,这里来一个简单的:从楼下到楼上有n 层台阶,一只青蛙要上楼,它有两种跳法,一种是一步一台阶,一种是一步二台阶,问给定台阶的数量n,这只青蛙有多少种跳法。
不太熟悉dp 或者经常做蓝桥杯题目的,可能一看这题就开始写搜索的代码了,确实,很多dp 也可以说是搜索的优化,因为搜索的状态太多,于是就将重复的搜索状态记录下来,这样就可以减少大量无用的搜索,这大概就是动态规划的作用吧。
从题中我们可以知道,除了第一层和第二层,每到一层台阶,都有两种到达这层台阶的方法,一种是从第n2 层台阶跳两层上来,一种是从第n1 层台阶跳一步上来,也就是说,到达第n 层台阶的方法数量就是到达第n1 层的方法数量加上到达第n2 层的方法数量,设dp(n) 为到达第n 层台阶的跳法,于是有以下递推公式 :

dp(n)={1(n=1)2(n=2)dp(n1)+dp(n2)(n>2)

咦,这不就是一个斐波那契吗,然后写个循环就出来了。
个人觉得斐波那契的动态规划解法只是一种循环写法,是个人都能想到,还不能和动态规划扯上什么关系,可能写完了也不知道动态规划是什么,这个例子应该能比较形象地体现出动态规划在算法的优化上的作用。
当然,动态规划可不全是斐波那契数列,这里只是讲个简单的例子。

斐波那契数列

题目链接

九度OJ 1387: 斐波那契数列

题意

输出斐波那契数列第n 项,其中n[1,70]

题解

没什么好解的,一个循环预处理,注意long long 就行了。

过题代码

#include <cstdio>
using namespace std;

#define LL long long

int main() {
    int n;
    LL fib[100];
    fib[0] = 0;
    fib[1] = 1;

    for(int i = 2; i <= 70; ++i) {
        fib[i] = fib[i - 1] + fib[i - 2];
    }
    while(scanf("%d", &n) != EOF) {
        printf("%lld\n", fib[n]);
    }

    return 0;
}

矩阵快速幂

dp 已经这么快了,还可以优化?看标题↑。
看到这里,还记得前面提过的快速幂吗?n 个相同的数相乘可以用快速幂,n 个相同的矩阵相乘,自然也可以用快速幂。可这和dp 有什么关系呢?主要是因为:部分dp 是可以推导出某个能用矩阵相乘的形式表示的递推公式,斐波那契数列就是一个很好的例子。

算法思想

快速幂的背景是大数和取模上面已经说过,不再重复,这里来看如何把斐波那契数列的递推公式表示成矩阵相乘的形式:

fib(n)=fib(n1)+fib(n2)

(fib(n)fib(n1))=(1110)(fib(n1)fib(n2))

(fib(n)fib(n1))=(1110)n1(fib(1)fib(0))(n>0)

对于这种形式的递推公式的矩阵构造应该能够很轻松看得出来吧,其实快速幂的运用不仅仅是在这样的递推公式上,对于某些图将其邻接矩阵表示出来,也可以构造出一个矩阵快速幂,其算法时间复杂度为O(k3logn)k

斐波那契矩阵快速幂裸题

题目链接

51Nod 1242: 斐波那契数列的第N

题意

求斐波那契数列第n 项值,如果太大则对109+9 取模。

题解:矩阵快速幂

过题代码

#include <cstdio>
#include <cstring>
using namespace std;

#define LL long long
const LL MOD = 1000000009;
const int SIZE = 2;

struct Matrix {
    LL num[SIZE][SIZE];

    Matrix() {
        memset(num, 0, sizeof(num));
        for(int i = 0; i < SIZE; ++i) {
            num[i][i] = 1;
        }
    }
    void Set() {
        num[0][0] = num[0][1] = num[1][0] = 1;
        num[1][1] = 0;
    }
    void operator*=(const Matrix &b) {
        LL ans[SIZE][SIZE];
        for(int i = 0; i < SIZE; ++i) {
            for(int j = 0; j < SIZE; ++j) {
                ans[i][j] = 0;
                for(int k = 0; k < SIZE; ++k) {
                    ans[i][j] = (ans[i][j] + num[i][k] * b.num[k][j] % MOD) % MOD;
                }
            }
        }
        memcpy(num, ans, sizeof(ans));
    }
};

void mi(Matrix &res, LL n) {
    Matrix ans;
    for(; n != 0; n >>= 1) {
        if((n & 1) == 1) {
            ans *= res;
        }
        res *= res;
    }
    memcpy(res.num, ans.num, sizeof(ans.num));
}

int main() {
    LL n;
    Matrix matrix;
    scanf("%lld", &n);

    matrix.Set();
    mi(matrix, n - 1);

    printf("%lld", matrix.num[0][0] % MOD);

    return 0;
}

可以看到,这里的mi 函数与上面快速幂的写法几乎完全相同,只是定义了矩阵的结构体和矩阵乘法而已,当然矩阵乘法的写法肯定是要像能够默下来那样熟练,最好快速幂的写法也能默得下来,毕竟应该算是一种基础的算法,题目才不会这样告诉你:请用矩阵快速幂解决这道动态规划问题。

矩阵快速幂应用

说了这么多,就来检查一下是否理解了矩阵快速幂吧,最好看完下面的题目之后能够自己先想出怎么运用矩阵快速幂,再来看题解。

Okabe and El Psy Kongroo

题目链接

Codeforces 821E: Okabe and El Psy Kongroo

题意

有一个人要出去散步,但是他只能在安全的区域内散步,将他的散步区域限定在平面直角坐标系中的第一象限,从坐标(0,0) 开始,往右走到(k,0) 的位置,问安全的走法有多少种。
其中安全的区域有n 段,每段用ai,bi,ci(i[1,n]) 三个数表示:从x=aix=bi 之间,只能在y[0,ci] 之间散步。已知他只能从坐标(x,y) 走到(x+1,y1),(x+1,y),(x+1,y+1) 三个坐标上。
n[1,100],a,b,k[1,1018],c[0,15] 数据保证a1=0,ankbn,且当i[2,n] 时,ai=bi1。由于结果可能非常大,将求得的结果对109+7 取模。

题解

这题是不是与上面的青蛙跳台阶很像?如果熟练的话可以看得出来,这是一个可以用dp 优化的搜索的题目,从走路的方式很容易看出,如果设走到坐标(x,y) 的走法为dp(x,y),则:

dp(x,y)=dp(x1,y)+dp(x1,y1)+dp(x1,y+1)

其中任何一项的纵坐标超过c,其值都为0。
递推公式有了,不过这是一个二维的dp,并不像前面那个一维的容易写出矩阵(当然要构造矩阵啦,k 的范围可是1018 次方,横坐标一个一个推过去肯定超时),这题可以将每一个横坐标对应的所有点看作是一个状态,那么对于每一段(注意这里每一段对应的都是不同的矩阵)的状态转移方程就是:
dp(x,0)=dp(x1,0)+dp(x1,1)dp(x,1)=dp(x1,0)+dp(x1,1)+dp(x1,2)dp(x,2)=dp(x1,1)+dp(x1,2)+dp(x1,3)dp(x,c)=dp(x1,c1)+dp(x1,c)

这样这个矩阵就容易推导多了:
(dp(x,0)dp(x,1)dp(x,2)dp(x,c))=(1100111001100001)(dp(x1,0)dp(x1,1)dp(x1,2)dp(x1,c))

最后两个细节上的问题,一个是在两段的分界线上,c 值不同应当将两个c 值中最小的那个c 值以上的所有状态都设为0。另一个是算法的时间复杂度,显然被分的段数越多,c 的取值范围越大,算法的时间复杂度也越高,整体的时间复杂度大概为O(nc3logk),庆幸nc3 才337500,logk 也只在64 左右。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <sstream>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long

const int maxn = 16;
const LL MOD = 1000000007;
struct Matrix {
    LL num[maxn][maxn];

    Matrix() {
        memset(num, 0, sizeof(num));
        for(int i = 0; i < maxn; ++i) {
            num[i][i] = 1;
        }
    }
    void Set(int Size) {
        memset(num, 0, sizeof(num));
        for(int i = 0; i < Size; ++i) {
            for(int j = i - 1; j <= i + 1; ++j) {
                if(j >= 0 && j < Size) {
                    num[i][j] = 1;
                }
            }
        }
    }
    void mult(const Matrix &b, const int &Size) {
        LL ans[maxn][maxn];
        for(int i = 0; i < Size; ++i) {
            for(int j = 0; j < Size; ++j) {
                ans[i][j] = 0;
                for(int k = 0; k < Size; ++k) {
                    ans[i][j] = (ans[i][j] + num[i][k] * b.num[k][j] % MOD) % MOD;
                }
            }
        }
        memcpy(num, ans, sizeof(ans));
    }
};
LL N, K, a, b, c, Ans[2][maxn], now;
Matrix tmp;

void mi(Matrix &res, LL n, int Size) {
    Matrix ans;
    for(; n != 0; n >>= 1) {
        if((n & 1) == 1) {
            ans.mult(res, Size);
        }
        res.mult(res, Size);
    }
    memcpy(res.num, ans.num, sizeof(ans.num));
}

void setZero(LL *num, int n) {
    for(int i = n + 1; i < maxn; ++i) {
        num[i] = 0;
    }
}

int main() {

    #ifdef LOCAL
    freopen("test.txt", "r", stdin);
    #endif // LOCAL

    ios::sync_with_stdio(false);

    Ans[now][0] = 1;

    scanf("%I64d%I64d", &N, &K);

    while(N--) {
        scanf("%I64d%I64d%I64d", &a, &b, &c);
        if(a < K) {
            if(b > K) {
                b = K;
            }
            setZero(Ans[now], c);
            tmp.Set(c + 1);
            mi(tmp, b - a, c + 1);
            for(int i = 0; i <= c; ++i) {
                Ans[!now][i] = 0;
                for(int j = 0; j <= c; ++j) {
                    Ans[!now][i] = (Ans[!now][i] + tmp.num[i][j] * Ans[now][j] % MOD) % MOD;
                }
            }
            now = !now;
            setZero(Ans[now], c);
        }
    }

    printf("%I64d\n", Ans[now][0]);

    return 0;
}
阅读更多
版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/CSDNjiangshan/article/details/75136730
个人分类: 博文
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭
关闭