POJ 3233 Matrix Power Series(求矩陣冪的和——分塊矩陣快速冪 or 二分遞迴+矩陣快速冪)
| Time Limit: 3000MS | Memory Limit: 131072K |
| Total Submissions: 21451 | Accepted: 8993 |
Description
Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.
Input
The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k
Output
Output the elements of S modulo m in the same way as A is given.
Sample Input
2 2 4 0 1 1 1
Sample Output
1 2 2 3
題意：給定 n*n 的矩陣A和正整數 k 和 m 。 求矩陣A的冪之和。 S = A + A
題解：
方法一：構造法
因為分塊矩陣依然滿足矩陣乘法，所以可以構造成分塊矩陣冪求解，詳見 《挑戰》 P.205
程式碼如下：
#include <iostream> #include <cstdio> #include <cstring> #include <cmath> #include <vector> #include <algorithm> using namespace std; #define LL long long #define maxn 35 typedef vector<LL> vec; typedef vector<vec> mat; int m; mat mul(mat &A, mat &B)//矩陣乘法 { mat C(A.size(), vec(B[0].size())); for(int i=0; i<A.size(); ++i) { for(int k=0; k<B.size(); ++k) { for(int j=0; j<B[0].size(); ++j) C[i][j] = (C[i][j] + A[i][k]*B[k][j]) % m; } } return C; } mat quick_pow(mat A,int n)//快速冪 { mat B(A.size(), vec(A.size())); for(int i=0; i<A.size(); ++i) B[i][i]=1; while(n>0) { if(n&1) B = mul(B, A); A = mul(A, A); n >>= 1; } return B; } int main() { int n,k; LL a; while(scanf("%d%d%d",&n,&k,&m)!=EOF) { mat A(n*2, vec(n*2)); for(int i=0; i<n; ++i) { for(int j=0; j<n; ++j) { scanf("%lld",&a); A[i][j] = a; } A[n+i][i] = 1; A[n+i][n+i] = 1;//構造出分塊矩陣下方的兩個單位矩陣塊 } A = quick_pow(A, k+1);// I+A+A^2+...+A^k，構造後的 S(k) = I+A+...+A^k-1 ，所以這裡我們要求出S(k+1) for(int i=0; i<n; ++i) { for(int j=0; j<n ;++j) { a = A[n+i][j] % m; if(i==j)//減去最早加上的單位矩陣 a = (a+m-1) % m; if(j==n-1) printf("%lld\n",a); else printf("%lld ",a); } } } return 0; }
方法二： 二分+矩陣快速冪
我們知道對於矩陣 A^x 可以用矩陣快速冪 快速的求解。 我們要考慮的是怎麼減少 求解 S = A + A^2 + A^3 +...+A^k 過程中的 相加 這個操作的耗時。 下面我們把 矩陣A 看成一個普通的 數字A
對於S(10)我們有：
S(10) = A + A^2 + A^3 +A^4 + A^5 + A^5 * (A + A^2 +A^3 +A^4 + A^5 )
對於S(5)我們有：
S(5) = A + A^2 +A^3 + A^3 * (A + A^2)
對於S(2)我們有：
S(2) = A + A * (A)
題中已經給出了 A 的值， 那麼我們可以通過遞歸回溯得到S(2)的值，求出A^3的值，就可以通過上面的式子得到S(5)的值， 進而得到S(10)的值。 那麼對於每個S(k)我們都可以這樣二分遞迴求解。 注意這裡的S(n)當n為奇數需要多加上一次A^(n/2+1)的值（ 例如上面的S(5) ）。 最後加上矩陣的快速冪 與 加法操作就行了。
程式碼如下：
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define LL long long
#define maxn 40
int n,m,k;
struct Node
{
LL mat[maxn][maxn];
}A;
Node mul(Node B, Node C)//矩陣乘法
{
Node D;
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
D.mat[i][j] = 0;
for(int i=0; i<n ;++i)
{
for(int k=0; k<n; ++k)
{
for(int j=0; j<n; ++j)
D.mat[i][j] = (D.mat[i][j] + B.mat[i][k]*C.mat[k][j]) % m;
}
}
return D;
}
Node quick_pow(Node B,int x)//快速冪
{
Node C;
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)//初始化為單位矩陣
C.mat[i][j] = (i==j?1:0);
while(x>0)
{
if(x&1)
C = mul(C, B);
B = mul(B, B);
x >>= 1;
}
return C;
}
Node Add(Node B, Node C)//矩陣加法
{
Node D;
for(int i=0; i<n; ++i)
{
for(int j=0; j<n; ++j)
D.mat[i][j] = (B.mat[i][j] + C.mat[i][j]) % m;
}
return D;
}
Node dfs(int x)
{
if(x==1)
return A;
Node ans = dfs(x/2), temp;
if(x&1){
temp = quick_pow(A, x/2+1);
ans = Add(ans, mul(temp, ans));
ans = Add(temp, ans);//奇數情況要加上中間的矩陣
}
else{
temp = quick_pow(A, x/2);
ans = Add(ans, mul(temp, ans));
}
return ans;
}
int main()
{
while(scanf("%d%d%d",&n,&k,&m)!=EOF)
{
for(int i=0; i<n; ++i)
{
for(int j=0; j<n; ++j)
scanf("%lld",&A.mat[i][j]);
}
Node ans = dfs(k);
for(int i=0; i<n; ++i)
{
for(int j=0; j<n-1; ++j)
printf("%lld ",ans.mat[i][j]);
printf("%lld\n",ans.mat[i][n-1]);
}
}
return 0;
}