POJ3233 矩阵的N次幂求和 二分

k 为偶数：A^k = A ^(k/2) * A ^(k/2);

k 为奇数： A^k = A ^ (k/2) * A^(k/2) * A

  

n  = 2k 为偶数

A^1 + A^2 + A^3 +....... A^(2k) = (A + A^2 + A^3 + .....A ^ k) + A^k * (A + A^2 + A^3 + .....A ^ k);

n = 2k + 1为奇数

A^1 + A^2 + A^3 +....... A^(2k + 1) = (A + A^2 + A^3 + .....A ^ k) +A ^(k + 1) + A^(k + 1) * (A + A^2 + A^3 + .....A ^ k);

#include <iostream>#include <cstdio>#include <cstring>#include <cmath>using namespace std;struct MATRIX{int m[35][35];};int m;MATRIX MOD(MATRIX num1,int n){int i,j;MATRIX res;for(i = 0; i < n; i++){for(j = 0; j < n; j++){res.m[i][j] = res.m[i][j]  %  m;}}return res;}MATRIX  mul(MATRIX num1, MATRIX num2,int n){int x;int i,j;MATRIX res;memset(&res,0,sizeof(res));for(i = 0; i <n; i++){for(j = 0; j < n; j++){for(x = 0; x < n; x++){res.m[i][j] += (num1.m[i][x] * num2.m[x][j]) % m;res.m[i][j] %= m;}}}return res;}MATRIX sum(MATRIX num1,MATRIX  num2, int n){int i,j;MATRIX result;for(i = 0; i < n; i++){for(j = 0; j < n; j++){result.m[i][j] = (num1.m[i][j] + num2.m[i][j]) % m;}}return result;//Why not res}MATRIX POW(MATRIX mat,int n,int t){if(t == 1){return  mat;}MATRIX res;if(t % 2 == 0){res = POW(mat,n, t/2);MATRIX current;memcpy(¤t, &res, sizeof(res));res = mul(current,res,n);}else{res = POW(mat,n,t/2);res = mul(res,res,n);res = mul(res,mat,n);}return res;}MATRIX total(MATRIX mat,int k,int n){MATRIX current;MATRIX res;if(k == 1){return mat;}if(k % 2 == 0){current = total(mat,k/2,n);res = sum(current,mul(current , POW(mat,n,k/2),n),n);return res;}else{current = total(mat,k/2,n);  MATRIX flag = POW(mat,n,k/2 + 1);  MATRIX flag1;  flag1 = sum(current, flag,n);  res = sum(flag1,mul(flag,current,n),n);  return res;}}int main(){int n,k;while(scanf("%d %d %d",&n,&k,&m) != EOF){MATRIX num;MATRIX res;memset(&num,0,sizeof(num));memset(&res,0,sizeof(res));int i,j;for(i = 0; i < n; i++){for(j = 0; j < n; j++){scanf("%d",&num.m[i][j]);}}MOD(num, n);res = total(num, k,n);MOD(num, n);for(i = 0; i < n; i++){for(j = 0; j < n; j++){printf("%d%c",res.m[i][j],j != n-1? ' ':'\n');}}}return 0;}