POJ 3233 二分求等比数列矩阵快速幂_二分幂求等比数列-优快云博客

本文提供了一种使用矩阵快速幂解决POJ 3233问题的方法，通过递归分解求解S=(A+A^2+A^3+…+A^k)modm，介绍了矩阵乘法、加法及快速幂运算的实现。

题目链接：http://poj.org/problem?id=3233
题意：求 S =( A + A^2 + A^3 + … + A^k) mod m

二分原理(转自AcDreamer大佬)：
（1）当这里写图片描述时，
（2）当时，那么有

（3）当这里写图片描述时，那么有

代码：

#include <cstdio>
#include <iostream>
#include <cstring>
#define sf scanf
#define pf printf
using namespace std;
const int maxn = 35;
int n;
__int64 mod;
struct Mat{
    __int64 C[maxn][maxn];
    Mat(){memset(C,0,sizeof C);for(int i = 0;i < n;++i) C[i][i] = 1;}
    void out(){
        for(int i = 0;i < n;++i){
            pf("%I64d",C[i][0]);
            for(int j = 1;j < n;++j) pf(" %I64d",C[i][j]);
            pf("\n");
        }
    }
};
Mat operator*(const Mat& A,const Mat& B){
    Mat ret;
    for(int i = 0;i < n;++i){
        for(int j = 0;j < n;++j){
            ret.C[i][j] = 0;
            for(int k = 0;k < n;++k){
                ret.C[i][j] += A.C[i][k]*B.C[k][j];
                ret.C[i][j] %= mod;
            }
        }
    }
    return ret;
}
Mat operator+(const Mat& A,const Mat& B){
    Mat ret;
    for(int i = 0;i < n;++i)
        for(int j = 0;j < n;++j)
            {ret.C[i][j] = A.C[i][j] + B.C[i][j];ret.C[i][j] %= mod;}
    return ret;
}
Mat operator^(Mat A,int k){
    Mat ret;
    while(k){
        if(k & 1) ret = ret * A;
        A = A * A;
        k = k >> 1;
    }
    return ret;
}

Mat getSum(const Mat A,int k){
    if(k == 1){
        return A;
    }
    else{
        Mat ret = getSum(A,k / 2);
        if(k & 1){
            Mat tmp = A ^ (k / 2 + 1);
            return ret + (ret * tmp) + tmp;
        }
        else{
            return ret + ret * (A ^ (k / 2));
        }
    }
}

int main(){
    int k;
    while( ~sf("%d%d%I64d",&n,&k,&mod)){
        Mat st;
        for(int i = 0;i < n;++i)
            for(int j = 0;j < n;++j) {sf("%I64d",&st.C[i][j]);st.C[i][j] %= mod;}
        getSum(st,k).out();
    }
    return 0;

}