bzoj[2510] 弱题（矩阵快速幂+dp）

做一道矩阵优化dp的练习题：弱题

2510: 弱题

Time Limit: 10 Sec Memory Limit: 128 MB
Submit: 474 Solved: 257
[Submit][Status][Discuss]

Description

有M个球，一开始每个球均有一个初始标号，标号范围为1～N且为整数，标号为i的球有$a_i$个，并保证$Σa_i = M$。
每次操作等概率取出一个球（即取出每个球的概率均为$\frac 1 M$），若这个球标号为$k（k < N）$，则将它重新标号为$k + 1$；若这个球标号为$N$，则将其重标号为1。（取出球后并不将其丢弃）
现在你需要求出，经过$K$次这样的操作后，每个标号的球的期望个数。

Input

第1行包含三个正整数$N，M，K$，表示了标号与球的个数以及操作次数。
第2行包含$N$个非负整数$a_i$，表示初始标号为$i$的球有$a_i$个。

Output

应包含N行，第i行为标号为i的球的期望个数，四舍五入保留3位小数。

Sample Input

2 3 2
3 0

Sample Output

1.667
1.333

样例说明

第1次操作后，由于标号为2球个数为0，所以必然是一个标号为1的球变为标号为2的球。所以有2个标号为1的球，有1个标号为2的球。

第2次操作后，有$\frac 1 3$的概率标号为2的球变为标号为1的球（此时标号为1的球有3个），有$\frac 2 3$的概率标号为1的球变为标号为2的球（此时标号为1的球有1个），所以标号为1的球的期望个数为$\frac 1 3*3+\frac 2 3 = \frac 5 3$。同理可求出标号为2的球期望个数为$\frac 4 3$。

【数据规模与约定】

对于10%的数据，$N ≤ 5, M ≤ 5, K ≤ 10$；

对于20%的数据，$N ≤ 20, M ≤ 50, K ≤ 20$；

对于30%的数据，$N ≤ 100, M ≤ 100, K ≤ 100$；

对于40%的数据，$M ≤ 1000, K ≤ 1000$；

所有数据，$N ≤ 1000, M ≤ 100,000,000, K ≤ 2,147,483,647$

$\\$
$\\$
$\\$
$\\$
$\\$
$\\$
$\\$

解：

看到k的范围，很容易就会想到矩阵快速幂优化dp。因为我们只能把$k$来$log$掉。dp中做到这一点的就是矩阵。但是有一个问题，这里的转移矩阵边长为1，可能$n^3$不是很过的了。那我们要考虑优化矩阵乘法：先看一种$n=3$的情况：

有没有发现这个矩形的特殊之处？是的，这个矩形每一行的值都相同，且位置正好错开1。我们把它叫做循环矩阵，很显然如果我们知道第一行就可以直接推出整个矩阵。但这并不是循环矩阵最好的性质。

这道题可以利用循环矩阵的性质其实是：如果$A$是循环矩阵，$B$也是循环矩阵，那么$A+B,A-B,A*B$都是循环矩阵。加减是很显然的，我们证明一下乘法：

$C(i,j)=\sum_{k=1}^n A(i,k)*B(k,j)$

$C(i+1,j+1)$
$=\sum_{k=1}^nA(i+1,k+1)*B(k+1,j+1)$
$=\sum_{k=1}^nA(i,k)*B(k,j)$

由图可知$A(i,0)=A(i,n)$，同理$B(i,0)=B(i,n)$。

得证。

所以这个转移矩阵无论多少次方都是一个循环矩阵，那么我们算乘法时只用计算第一行，然后$n^2$推其余行就行。我们的矩乘做到了$n^2$。但是这样虽然时间过了，但是如果把1000*1000的矩阵封装在结构体里，还存double，亲测上来就RE了，连输入的机会都不给。可能要用一个函数来写矩乘。

后来才发现自己傻逼了，既然我们只需要第一行就可以推整体，那么我们为什么要推整体，何不存就只存第一行？于是我们使用一维数组存一个矩阵，又可以愉快地用结构体写乘法了。而且最后答案也只需要一行，我们其实都不需要推出矩阵。

下面附上代码：

#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;
int n,m,k;
struct lxy{
    double m[1005];
    lxy operator * (const lxy &q) const {
      lxy ret;
      for(int i=1;i<=n;i++)
        ret.m[i]=0;
      for(int i=1;i<=n;i++)
        for(int k=1;k<=n;k++)
          ret.m[i]+=m[k]*q.m[(i-k+n)%n+1];//这里经过公式推导发现竖行可以使用第一排的某个顺序替代
      return ret;
    }
}A,key;

lxy quick(lxy B,int p)
{
    lxy C;
    for(int i=1;i<=n;i++)
      C.m[i]=0;
    C.m[1]=1;
    while(p!=0)
    {
        if((p&1)==1) C=C*B;
        B=B*B;
        p=p>>1;
    }
    return C;
}

int main()
{
    scanf("%d%d%d",&n,&m,&k);
    for(int i=1;i<=n;i++)
      scanf("%lf",&A.m[i]);
    for(int i=1;i<=n;i++)
    {
        if(i==1) key.m[i]=(m-1.0)/m;
        else if(i==2) key.m[i]=1.0/m;
        else key.m[i];
    }
    key=quick(key,k);
    A=A*key;
    for(int i=1;i<=n;i++)
      printf("%.3f\n",A.m[i]);
    return 0;
}