POJ-3735-Training little cats-構造矩陣+矩陣快速冪+稀疏矩陣乘法優化
http://poj.org/problem?id=3735
題意:
n只貓,三種命令:
1、第i只貓吃掉所有花生;
2、第i只貓得到一個花生;
3、交換第i,j只貓的花生;
先由k個 這些命令組成一個操作序列
然後重複操作序列m次,
n,k<=100,m<=1e9
m的次數那麼大,可以用構造矩陣,然後用快速冪的方法
引用大神http://blog.csdn.net/magicnumber/article/details/6217602的圖片
初始我們開一個1行n+1列的矩陣存放每隻貓的初始值 第n+1行的前n個元素表示每隻貓的花生數
org=【
1 0 0 0
0 1 0 0
0 0 1 0
x y z 1
】
在n+1階單位矩陣的基礎上,
給第1 只貓加1個花生的操作矩陣G為:
1 0 0 0
0 1 0 0
0 0 1 0
1 0 0 1 第n+1行的第i列加1
讓第2只貓吃完所有花生的操作矩陣E為:
1 0 0 0
0 0 0 0
0 0 1 0
0 0 0 1 第i列清空為0
交換第1 第2 只貓的操作矩陣S為:
0 1 0 0
1 0 0 0 交換第i第j列元素
0 0 1 0
0 0 0 1
********************************************************
我們可以發現,讓org初始矩陣左乘一個G矩陣,得到
【
1 0 0 0
0 1 0 0
0 0 1 0
(x+1) Y Z 1
】
恰好就是X加了1
我們可以發現,讓org初始矩陣左乘一個E矩陣,得到
【
1 0 0 0
0 1 0 0
0 0 1 0
X 0 Z 1
】
恰好就是吃光了Y
我們可以發現,讓org初始矩陣左乘一個S矩陣,得到
【
1 0 0 0
0 1 0 0
0 0 1 0
y x Z 1
】
恰好就是交換了Y和X的花生
******************************************************
那麼對於k^m次操作,我們先處理好前k次的操作
也就是用org初始矩陣不斷乘G/E/S矩陣 得到前k次的結果 (PS:由於我們可以直接知道k次每次乘了操作矩陣後的結果(交換兩列、清空行,指定位置加1),就可以不必模擬k次,直接得到k次的結果了 )
接下來就是對k次得到結果矩陣X 連續左乘m次。。。這裡就用個快速冪就OK,然而每次乘法O(N^3),log(m)≈26 所以還需要優化一下乘法部分,
因為顯然大部分元素為0,是稀疏矩陣,我們可以用係數矩陣乘法優化
一般矩陣乘法:
Matrix mul(Matrix a, Matrix b) //矩陣相乘
{
Matrix res;
for(int i = 0; i < k; i++)
for(int j = 0; j < k; j++)
{
res.mat[i][j] = 0;
for(int t = 0; t < k; t++)
{
res.mat[i][j] += a.mat[i][t] * b.mat[t][j];
res.mat[i][j] %= mod;
}
}
return res;
}
稀疏矩陣:
matrix mul(matrix a,matrix b)
{
matrix res;
memset(res.m,0,sizeof(res.m));
__int64 i,j,k;
for (i=1;i<=n+1;i++)
{
for (j=1;j<=n+1;j++)
{
if (a.m[i][j]) //稀疏矩陣乘法(非O(n^3)), 即知道答案為零直接跳過
for (k=1;k<=n+1;k++)
{
res.m[i][k]+=a.m[i][j]*b.m[j][k];
}
}
}
return res;
}
從而解決問題,
難點是稀疏矩陣優化/構造操作矩陣
程式碼:
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <iostream>
#include <queue>
#include <map>
#include <set>
#include <vector>
using namespace std;
struct matrix
{
__int64 m[105][105];
};
matrix unit_org;
matrix unit;
__int64 n,k;
int m;
void init()
{
char op;
__int64 num,num1,num2,i,j;
for (j=1;j<=k;j++)
{
scanf("%c",&op);
if (op=='g')
{
scanf("%I64d",&num);
getchar();
unit.m[n+1][num]++; //第num個+1
}
if (op=='e')
{
scanf("%I64d",&num);
getchar();
for (i=1;i<=n+1;i++)
unit.m[i][num]=0; //第num列置為零
}
if (op=='s')
{
scanf("%I64d%I64d",&num1,&num2);
getchar();
if (num1==num2) continue;
for (i=1;i<=n+1;i++) //交換兩列
{
num=unit.m[i][num1];
unit.m[i][num1]=unit.m[i][num2];
unit.m[i][num2]=num;
}
}
}
}
matrix mul(matrix a,matrix b)
{
matrix res;
memset(res.m,0,sizeof(res.m));
__int64 i,j,k;
for (i=1;i<=n+1;i++)
{
for (j=1;j<=n+1;j++)
{
if (a.m[i][j]) //稀疏矩陣乘法(非O(n^3)), 即知道答案為零直接跳過
for (k=1;k<=n+1;k++)
{
res.m[i][k]+=a.m[i][j]*b.m[j][k];
}
}
}
return res;
}
matrix pow_mi(matrix a,int b)
{
matrix res=unit_org;
while(b)
{
if (b&1)
res=mul(res,a);
a=mul(a,a);
b>>=1;
}
return res;
}
int main()
{
__int64 i,j;
while(scanf("%I64d%d%I64d",&n,&m,&k)!=EOF)
{
getchar();
if (!n&&!m&&!k) break;
memset(unit.m,0,sizeof(unit.m));
for (i=1;i<=101;i++)
unit.m[i][i]=1;
unit_org=unit;
init();
matrix ans= pow_mi(unit,m);
for (i=1;i<=n;i++)
{
if (i!=1) printf(" ");
printf("%I64d",ans.m[n+1][i]);
}
printf("\n");
}
return 0;
}