POJ月赛题目Matrix Power Series

来源：互联网发布：idg 知乎编辑：程序博客网时间：2024/05/01 08:17

描述

Given a n × n matrix A anda positive integer k, find the sum S = A + A² + A³ +… + A^k.

输入

The input contains exactly one test case. The first line ofinput contains three positive integers n (n ≤ 30), k (k ≤ 10^9) and m (m < 10^4).Then follow n lines each containing n nonnegative integers below 32,768, givingA’s elements in row-major order.

输出

Output the elements of S modulo m in the same way as A is given.

样例输入

2 2 4

0 1

1 1

样例输出

1 2

2 3

来源

POJ Monthly

算法：
先考虑以下和X= a + a² +a³ + … + a^k，即a是一个实数，求级数和X。一般有两种想法：

（1）利用霍纳（horner法则---不会的自己翻算法导论，并自觉面墙去）；

（2）即利用级数和公式，可以在 O(log n)内计算。

本题目解法：

在一开始时就需要考虑到数据规模，如果使用类似horner的算法，则复杂度是O(k * n *n*n果断到10¹³ 的复杂度级别，果断超时，利用级数求和时，需要做除法，两矩阵的除法存在的冲要条件是矩阵为两矩阵同型，分母方阵且行列式不为0，即将除法转化为乘上分母的逆矩阵。在本题中，就算逆矩阵存在，也不一定能使用这种算法计算出结果，因为要计算行列式的值，30 * 30的行列式的值，目前并没有什么靠谱的简单算法能快速计算。

大家都知道计算两个复数(a+bi)与(c+di)的乘积，一般的算法是使用定义，即使用4次乘法就OK。不知道大家是否还记得算法导论上有个题目，即只使用3次乘法去计算两个复数的积，算法大致如下：先计算a*c、b*d、(a+b)*(c+d)，显然，这三个乘积足够推导出(a+bi)*(c+di)的值。在本题中，一样可使用类似的拆项技巧，S = A + A² + A³ +… + A^k，设E 为单位方阵，考虑S+E=A⁰+A¹+ A² + A³ +… + A^k，当k = 3时，我们只需要计算(E+A)*(E+A²)，当k = 7 时，我们只需要计算(E+A)*(E+A²)*(E+A³)，对于1、3、7、15...这种特殊的k，我们可以速度地计算出来，但是，对于一般的k，我们能怎么办呢？答案很简单，拆项...，比如k=13，拆成
S+E=A⁰+A¹+ A² + A³ +… + A¹³= (E+A)*(E+A²)*(E+A³)+ A⁸(E+A)*(E+A²)+A¹²(E+A)...

写出这个式子的时候，相信大家都明白了，求值时，我们只需要先在O(log k * n * n * n)的时间复杂度内预先计算几个矩阵的值，然后可地推地把k减小，直到0（每次至少减一半）。于是，关于递归时间复杂度有如下的递推
T(n,k) = T(k/2) + O( n * n * n)，由主定理知：递归时间O(log k * n * n * n)，再加上预处理矩阵的时间，总时间仍然是O(log k * n * n * n)，根据题目给的数据，该复杂度肯定可以接受，以下是代码：

<span style="font-size:18px;">#include <cstdlib>#include <cstring>#include <algorithm>#include <iostream>using namespace std; #ifdef _MSC_VERtypedef __int64 lld;#elsetypedef long long lld;#endif class Matrix{public:    Matrix();    Matrix(const Matrix& rhs);public:    friend istream& operator >>(istream& in, Matrix& rhs);public:    Matrix& operator *= (const Matrix&rhs);    friend const Matrix operator * (constMatrix& lhs, const Matrix& rhs);    friend Matrix MakeIndenty(const lldiScale);    Matrix& operator = (const Matrix&rhs);    Matrix& operator ^= (lld Exp);    Matrix& operator %= (lld Mod);    Matrix& operator -=(const Matrix&rhs);    Matrix& operator += (const Matrix&rhs);    friend const Matrix operator % (constMatrix& lhs, const lld Mod);    friend Matrix operator ^ (const Matrix& lhs, const lld Exp);public:    lld m_Element[30][30];    lld m_iScale;}; lld Mod, iScale, Exp; //power matrixMatrix PM[32];Matrix PME[32];Matrix PMEMul[32];//matrix explld exp[32];  Matrix&Matrix::operator-=(const Matrix& rhs){    lld iScale = m_iScale;    for(int r= 0; r < iScale; ++r){        for(int c = 0; c< iScale ;++c){            m_Element[r][c] -=rhs.m_Element[r][c];            m_Element[r][c] += Mod;        }    }   *this %= Mod;    return *this;} Matrix& Matrix::operator +=(const Matrix& rhs){    lld iScale = m_iScale;    for(int r= 0; r < iScale; ++r){        for(int c = 0; c< iScale ;++c){            m_Element[r][c] +=rhs.m_Element[r][c];        }    }    *this %= Mod;    return *this;} Matrix operator + (constMatrix& lhs, const Matrix& rhs){    Matrix result = lhs;    result += rhs;    return result;} Matrix MakeIndenty(const lldiScale){    Matrix result;    result.m_iScale = iScale;    for(int i= 0; i< iScale; ++i){        result.m_Element[i][i] = 1;    }    return result;} Matrix::Matrix(){    memset(m_Element, 0, sizeof(m_Element));    m_iScale = 0;} Matrix::Matrix(const Matrix&rhs){    memcpy(m_Element, rhs.m_Element, sizeof(m_Element));    m_iScale = rhs.m_iScale;} istream& operator >>(istream& in, Matrix& rhs){    for(int r = 0; r< rhs.m_iScale; ++r){        for(int c = 0; c< rhs.m_iScale;++c){            in>>rhs.m_Element[r][c];        }    }    rhs %= Mod;    return in;} Matrix& Matrix::operator *=(const Matrix &rhs){    lld ele[30][30];    lld iScale = m_iScale;    for(int r = 0; r< iScale; ++r){        for(int c = 0; c< iScale; ++c){            ele[r][c] = m_Element[r][c];            m_Element[r][c] = 0;        }    }    for(int r= 0; r< iScale; ++r){        for(int c= 0; c< iScale; ++c){            for(int k= 0; k< iScale; ++k){                m_Element[r][c] += ele[r][k] *rhs.m_Element[k][c];            }        }    }    *this %= Mod;    return *this;} Matrix& Matrix::operator =(const Matrix& rhs){    memcpy(m_Element, rhs.m_Element,sizeof(m_Element));    m_iScale = rhs.m_iScale;    return *this;} const Matrix operator * (constMatrix& lhs, const Matrix& rhs){    Matrix result = lhs;    result *= rhs;    return result;} Matrix& Matrix::operator ^=(lld Exp){    if(Exp == 1){        //empty    }else if(!Exp){        *this = MakeIndenty(m_iScale);    }else{        Matrix xx = *this;        xx^=(Exp>>1);        xx*=xx;        if(Exp & 1){            xx*=*this;        }        *this = xx;    }    *this %= Mod;    return *this;} Matrix operator ^ (constMatrix& lhs, const lld Exp){    Matrix xx = lhs;    xx^=Exp;    return xx;} Matrix& Matrix::operator %=(lld Mod){    lld iScale = m_iScale;    for(int r = 0; r< iScale; ++r){        for(int c = 0; c< iScale; ++c){            m_Element[r][c] %= Mod;        }    }    return *this;} const Matrix operator % (constMatrix& lhs, const lld Mod){    Matrix xx = lhs;    xx %= Mod;    return xx;} Matrix Solve(const Matrix&PreFix, const lld k); int main(){    cin>>iScale>>Exp>>Mod;    PM[0].m_iScale = iScale;    cin>>PM[0];    for(int i= 1; i<= 31; ++i){        PM[i] = PM[i-1] * PM[i-1];    }    for(int i= 0; i<= 31; ++i){        PME[i] = PM[i];        PME[i] += MakeIndenty(iScale);    }    exp[0] =1;    for(int i= 1; i< 31; ++i){        exp[i] = (exp[i-1] << 1) + 1;    }    PMEMul[0] = PME[0];    for(int i= 1; i< 31; ++i){        PMEMul[i] = PME[i] * PMEMul[i-1];    }    cout<<endl<<endl;    Matrix result = Solve(MakeIndenty(iScale),Exp);    result -= MakeIndenty(iScale);    for(int r = 0; r< iScale; ++r){        for(int c = 0; c< iScale; ++c){           cout<<result.m_Element[r][c]<<' ';        }        cout<<endl;    } } Matrix Solve(const Matrix&PreFix, const lld k){    if(k == 1){        return PreFix * PMEMul[0];    }    lld idx = (lower_bound(exp, exp + 32, k) -exp);    --idx;    lld left = k - (1<<(idx+1));    if(left == 0){        return PreFix * PMEMul[idx] + PreFix *PM[idx + 1];    }else if(left < 0){        return PreFix * PMEMul[idx];    }    return PreFix * PMEMul[idx] + Solve(PreFix* PM[idx + 1], left);}</span>

以下是OJ运行结果截图：

0 0