矩阵快速幂好题PvZ once again

题意：

植物大战僵尸算个out的游戏了，原谅被出题逼疯了的跑来挖坟了。

会玩的请无视这一段直接看题目{

游戏中僵尸向你的房子进发，吃掉沿途遇到的植物进入你的房子 你就死翘了

你在土地上种植各种植物来攻击阻挡僵尸

手推车：放置在终点，僵尸走到面前会启动推倒一整行的僵尸

大蒜：可种植的一种植物，发出恶心的气味，僵尸咬了一口就会换到邻近的另一行（如果有相邻两行，那么移动到另外两行概率是相等的）

南瓜：单纯的肉盾 被僵尸啃的

耐久度K： 植物被咬了K口后被僵尸吃掉

如有其他对游戏的不理解请clarify

}

问题是这样的：

我们的院子变成了N行M列的，而且种满了大蒜（耐久度K）（图是我盗了 我不会这么无聊的）coming的僵尸只有一只（然而这只僵尸貌似发生了变异，它每啃一口植物，同一列相同种类的植物也被啃掉一口，一口一排的样子恩恩），初始位置在第S行，因为没有放置攻击性的植物，所以僵尸就一路吃了，于是出题者很想知道僵尸死在自上而下1-N号手推车的概率各是多少

一个整数T（表示T组数据）

接下来的T组数据

每组给定四个整数 N M K S

数据范围

T<=1000

0<N<=20

0<M<=1000

0<K<=1000

1<=S<=N

思路：最初的想法是某一列转化到下一列的时候总是有相同的过程，即上一列的任意一个位置到这一列的某个位置的概率和为到这个位置的概率，忽略上一列的概率，那么上一列的任意一个位置到这一列的某个位置的概率和在每一列的转化中总是相同的。我们可以构造一个转移矩阵，M[i][j]表示从上一列的j行转移到这一列的i行的概率；有m列，即是转移了m次，用矩阵快速幂可以很快的求得；

问题到这就转化为求这个转移矩阵；一开始我用DP的方法，由1次递推到k次，考虑时间复杂度，有n个起始位置，k次转移，有n行，有T组case，总复杂度为1000*1000*20*20，超时了。。。。

最终你会发现在求这个转移矩阵同样可以用 矩阵快速幂。

考虑到他转移的特殊性，他在当前行只能转移到相邻两行，同样的你可以先构造一个转移矩阵

假设n=5，转移矩阵为：       0    0.5    0    0     0   第一行只能由第二行转移过来并且概率为0.5

                                                  1    0    0.5    0     0    第二行可由第一行转移过来，概率为1，也可由第三行转移过来，概率为0.5

                                                  0    0.5    0    0.5   0

                                                  0     0     0.5   0    1

                                                  0     0      0     0.5   0

然后初始化起点矩阵，运用k次快速幂，就可以求得M[i][j]了；

不过，还要注意当n为1的时候，直接输出概率为1.。。。。。（被卡了半小时。。）

代码奉上：

#include <cstdio>
#include <cstring>
using namespace std;
int n, m, k, s;
struct M
{
    double a[21][21];
}b, res, tmp, tmp2;
M multiply(M x, M y)
{
       M temp;
       for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)temp.a[i][j]=0;
       for(int i=1;i<=n;i++)
       {
               for(int j=1;j<=n;j++)
               {
                       for(int l=1;l<=n;l++)
                       {
                               temp.a[i][j]+=x.a[i][l]*y.a[l][j];
                       }
               }
       }
       return temp;
}
void calc(int N)
{
     while(N)
     {
             if(N&1)
                    res=multiply(res,b);
             N>>=1;
             b=multiply(b,b);
     }
}
void calc2(int N)
{
     while(N)
     {
             if(N&1)
                    tmp2=multiply(tmp2,tmp);
             N>>=1;
             tmp=multiply(tmp,tmp);
     }
}
int main()
{
    int t; scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d%d",&n,&m,&k,&s);
        if(n==1)
        {
            printf("1.0000\n");
            continue;
        }
        //第一次快速幂
        memset(tmp.a,0,sizeof(tmp.a));
        //初始化第一次的转移矩阵
        for(int i = 1; i <= n; i++)
        {
            if(i-1==0)
                tmp.a[i][i+1]=0.5;
            else if(i==n)
                tmp.a[i][i-1]=0.5;
            else
            {
                if(i==2)
                    {tmp.a[i][i+1]=0.5;tmp.a[i][i-1]=1;}
                else if(i==n-1)
                    {tmp.a[i][i+1]=1;tmp.a[i][i-1]=0.5;}
                else
                    {tmp.a[i][i+1]=0.5;tmp.a[i][i-1]=0.5;}
            }
        }
        memset(tmp2.a,0,sizeof(tmp2.a));
        for(int i =1; i <= n; i++)tmp2.a[i][i]=1.0;
        calc2(k);
        for(int i = 1; i <= n; i++)//求解第二次的转移矩阵
        {
            double h[21],h2[21];
            memset(h,0,sizeof(h));
            memset(h2,0,sizeof(h2));
            h[i]=1.0;
            for(int j = 1; j <= n; j++)
            {
                for(int l = 1; l <= n; l++)
                {
                    h2[j]+=tmp2.a[j][l]*h[l];
                }
            }
            for(int j = 1; j <= n; j++)
                b.a[j][i]=h2[j];
        }
        //第二次快速幂
        double start[21], ans[21];
        memset(start, 0, sizeof(start)); start[s]=1;
        for(int i = 1; i <= n; i++)for(int j = 1; j <= n; j++)res.a[i][j]=0;
        for(int i = 1; i <= n; i++)res.a[i][i]=1.0;
        calc(m);
        for(int i = 1; i <= n; i++)
        {
            ans[i]=0;
            for(int j = 1; j <= n; j++)
            {
                ans[i]+=res.a[i][j]*start[j];
            }
        }
        int flag=0;
        for(int i = 1; i <= n; i++)
        {
            if(!flag)
            {
                printf("%.4f",ans[i]);
                flag=1;
            }
            else
                printf(" %.4f",ans[i]);
        }
        printf("\n");
    }
    return 0;
}

顺便把矩阵快速幂模板给了：

#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <iostream>
using namespace std;
#define N 1000
int k,n;//k次，n为数组大小
struct M
{
       int a[N][N];
}origin,res;

M multiply(M x,M y)
{
       M temp;
       memset(temp.a,0,sizeof(temp.a));
       for(int i=0;i<n;i++)
       {
           for(int j=0;j<n;j++)
           {
               for(int k=0;k<n;k++)
               {
                    temp.a[i][j]+=x.a[i][k]*y.a[k][j];
               }
           }
       }
       return temp;
}
void init()
{
     memset(res.a,0,sizeof(res.a));
     for(int i = 0; i < n; i++)res.a[i][i]=1;//将res.a初始化为单位矩阵
}
void calc(int m)
{
     while(m)
     {
         if(m&1)
                res=multiply(res,origin);
         m>>=1;
         origin=multiply(origin,origin);
     }
}
int main()
{
    while(cin>>n)
    {
            init();
            calc(k);
    }
    return 0;
}