POJ 3233 二分求等比数列 矩阵快速幂

    xiaoxiao2022-05-22  34

    题目链接: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; }
    转载请注明原文地址: https://ju.6miu.com/read-1109400.html

    最新回复(0)