【矩阵快速幂】——转载伟佳

转自http://njzwj.github.io/2016/05/30/matrix-fast-power

矩阵快速幂

在处理递推问题时，通常我们使用的是O(n)的算法。比如说计算Fibonacci数列的第k个元素的值，一般利用递推式${\mathrm{Fn=Fn-1+Fn-2Fn=Fn-1+Fn-2\; 从n=1,2n=1,2\; 的条件递推。递推过程可以表示为如下形式：\; \{Fn+1=Fn+Fn-1Fn=Fn\{Fn+1=Fn+Fn-1Fn=Fn\; 等效于下面的矩阵形式：\; (Fn+1Fn)=(FnFn-1)(1110)(Fn+1Fn)=(FnFn-1)(1110)\; 很显然下式也成立：\; (Fn+1Fn)=(F2F1)(1110)n-1(Fn+1Fn)=(F2F1)(1110)n-1}}^{}$

这样一来我们就把一个递推的问题转化成了一个矩阵乘法的问题。只要计算出构造出的变换矩阵的若干次幂，那么用它去乘初始状态的向量，就可以算出结果了。当然计算矩阵的乘法可以迭代计算，也就是说不断右乘变换矩阵的方法，不过这就和递推没什么区别了，对于$\mathrm{n=109n=109\; 或者以上的数据来说，会面临超时的危险。由于方阵乘法有优秀的结合律的性质所以我们自然想到了用快速幂来解决这类问题。}$

假设$\mathrm{XX\; 是一个方阵，要计算XyXy\; ，首先将yy\; 展开为二进制的形式：\; y=a0\ast 20+a1\ast 21+a2\ast 22+\dots +ai\ast 2iy=a0\ast 20+a1\ast 21+a2\ast 22+\dots +ai\ast 2i\; 其中aiai\; 表示yy\; 二进制第i位的数值，比如5可以表示为：\; 5=1\ast 20+0\ast 21+1\ast 225=1\ast 20+0\ast 21+1\ast 22\; 那么X14X14\; 可以表示为：\; X14=X2X4X8X14=X2X4X8\; 换言之任意整数次幂都能被表示为若干项X2nX2n\; 的乘积，而X2nX2n\; 可以通过自身迭代得到(X2n+1=(X2n)2X2n+1=(X2n)2\; )，在矩阵阶数固定的情况加计算矩阵n次幂的时间复杂度是O(logn)。}$

 

HDOJ 4686 - Arc of Dream

题目链接

An Arc of Dream is a curve defined by following function: $$\mathrm{AoD(n)=n-1\sum i=0aibiAoD(n)=\sum i=0n-1aibi\; where\; a0=A0a0=A0\; ai=ai-1\ast AX+AYai=ai-1\ast AX+AY\; b0=B0b0=B0\; bi=bi-1\ast BX+BYbi=bi-1\ast BX+BY\; What\; is\; the\; value\; of\; AoD(N)\; modulo\; 1,000,000,007?}$$

Input

There are multiple test cases. Process to the End of File. Each test case contains 7 nonnegative integers as follows:

N
A0 AX AY
B0 BX BY

N is no more than ${\mathrm{10181018\; ,\; and\; all\; the\; other\; integers\; are\; no\; more\; than\; 2\times 1092\times 109\; .}}^{}$

解题思路

根据题意构造如下变换：

求变换矩阵的n-1次幂即可。

代码

留意mod运算为负的情况。

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
const int maxn=5;
struct Matrix {
    ll m[maxn][maxn];
    Matrix () {
        clear();
    }
    void diag1() {
        for (int i=0;i<maxn;++i)
            for (int j=0;j<maxn;++j)
                m[i][j]=(i==j);
    }
    void clear() {
        memset(m, 0, sizeof(m));
    }
    void out() {
        for (int i=0;i<maxn;++i) {
            for (int j=0;j<maxn;++j) printf("%lld ", m[i][j]);
            printf("\n");
        }
    }
    Matrix operator*(const Matrix &rhs) const {
        Matrix ans;
        for (int i=0;i<maxn;++i) {
            for (int j=0;j<maxn;++j) {
                ans.m[i][j]=0;
                for (int k=0;k<maxn;++k) {
                    ans.m[i][j]+=m[i][k]*rhs.m[k][j];
                    ans.m[i][j]%=mod;
                }
            }
        }
        return ans;
    }
};

Matrix operator^(const Matrix& m, ll n) {
    Matrix base=m, res;
    res.diag1();
    while (n) {
        if (n&1) res=res*base;
        base=base*base;
        n>>=1;
    }
    return res;
}

int main() {
    ll a0, b0, ax, ay, bx, by, n;
    Matrix m;
    while (scanf("%lld%lld%lld%lld%lld%lld%lld", &n, &a0, &ax, &ay, &b0, &bx, &by)!=EOF) {
        m.clear();
        if (n==0) {
            printf("0\n");
            continue;
        }
        m.m[0][0]=ax%mod;
        m.m[4][0]=ay%mod;
        m.m[1][1]=bx%mod;
        m.m[4][1]=by%mod;
        m.m[0][2]=m.m[0][3]=(ax*by)%mod;
        m.m[1][2]=m.m[1][3]=(ay*bx)%mod;
        m.m[2][2]=m.m[2][3]=(ax*bx)%mod;
        m.m[4][2]=m.m[4][3]=(ay*by)%mod;
        m.m[3][3]=m.m[4][4]=1;
        m=m^(n-1);
        a0%=mod;
        b0%=mod;
        printf("%lld\n", (mod*5+(a0*m.m[0][3])%mod+(b0*m.m[1][3])%mod+
        ((a0*b0)%mod*m.m[2][3])%mod+((a0*b0)%mod*m.m[3][3])%mod+(m.m[4][3])%mod)%mod);
    }
    return 0;
}