矩阵乘+概率DP bzoj2676 Contra

问题 C: Contra
时间限制: 3 Sec 内存限制: 512 MB
题目描述
　　偶然间，chnlich发现了他小时候玩过的一个游戏“魂斗罗”，于是决定怀旧。但是这是一个奇怪的魂斗罗MOD。
　　有N个关卡，初始有Q条命。
　　每通过一个关卡，会得到u分和1条命，生命上限为Q。
　　其中u=min(最近一次连续通过的关数,R)。
　　若没有通过这个关卡，将会失去1条命，并进入下一个关卡。
　　当没有生命或没有未挑战过的关卡时，游戏结束，得到的分数为每关得到的分数的总和。
　　由于chnlich好久不玩这个游戏了，每条命通过每个关卡的概率均为p(0<=p<=1)，原先chnlich的最高分纪录是S。
　　现在chnlich想要知道，当p至少为多少时，chnlich期望获得的总分数能够超过原先的最高分。
输入
　　输入共一行，分别表示整数N，整数R，整数Q，原先的最高分整数S。
输出
　　输出共一行，若不存在这样的p，输出”Impossible.”（不包含引号），否则输出p（保留6位小数）。
样例输入
样例输入一
　　4 2 1 5

样例输出一
　　0.880606

样例输入二
　　12 3 2 12

样例输出二
　　0.687201

提示
数据规模和约定
　　对于20%的数据，N<=15
　　对于50%的数据，N<=10000
　　对于100%的数据，N<=10^8,1<=R<=20,1<=Q<=5,保证S是一个可能出现的分数。

考虑二分答案。问题在于如何快速算出二分出的答案最后的期望。
f[i][j]为剩余i条命，连胜场数为j（就是得的分数为j）时的概率。考虑把两维压缩成一维，设一个编号，这个的概率转移很明显是可以写成矩阵的。
考虑倒推，这样有效避免了在中间算答案的困难。最后只需要输出初始状态。
矩阵快速幂可以快速解决问题。那么期望该怎么转移呢。
看着两个矩阵
（初始矩阵）
0 0 0 0 0 1
（转移矩阵）
0 0 0 0 0 0 //考虑把不同状态之间的概率转移写进左上角5*5的矩阵里。
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
1 2 3 4 5 1 //这一排表示当前为多少场连胜×概率，也就是达到当前状态时能得多少分的期望。
每次乘完初始矩阵里每一位就是所对应状态乘到那一位时的期望。
小优化：矩阵太大，而且有不少0，所以有0时就不用乘了。

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
#define eps 1e-15
using namespace std;
int n,r,q,s,tot,sz,id[6][24];
struct node
{
    double f[105][105];
    node(){memset(f,0,sizeof(f));}
    friend node operator *(node a,node b)
    {
        node c;
        for(int i=0;i<sz;i++)
            for(int k=0;k<sz;k++)
            {
                if(a.f[i][k]<eps)continue;
                for(int j=0;j<sz;j++)
                {
                    if(b.f[k][j]<eps)continue;
                    c.f[i][j]+=a.f[i][k]*b.f[k][j];
                }
            }
        return c;
    }
};
bool check(double x)
{
    node a,ans,b;a.f[tot][tot]=b.f[0][tot]=1;
    for(int i=1;i<=q;i++)
        for(int j=1;j<=r;j++)
        {
            a.f[tot][id[i][j]]+=j*x;
            if(i<q&&j<r)a.f[id[i+1][j+1]][id[i][j]]+=x;
            else if(i<q)a.f[id[i+1][j]][id[i][j]]+=x;
            else if(j<r)a.f[id[i][j+1]][id[i][j]]+=x;
            else a.f[id[i][j]][id[i][j]]+=x;
            if(i>1)a.f[id[i-1][1]][id[i][j]]=(1.0-x);
        }
    for(int i=0;i<=tot;i++)ans.f[i][i]=1;
    int m=n;
    while(m)
    {
        if(m&1)ans=ans*a;
        a=a*a;
        m>>=1;
    }
    b=b*ans;
    double sum=b.f[0][id[q][1]];
    return sum-s>eps;
}
int main()
{
    scanf("%d%d%d%d",&n,&r,&q,&s);sz=r*q+1;
    for(int i=1;i<=q;i++)
       for(int j=1;j<=r;j++)
            id[i][j]=tot++;
    if(!check(1.0))
    {
        printf("Impossible.\n");
        return 0;
    }
    double l=0,r=1.0,mid;
    while(r-l>eps)
    {
        mid=(l+r)/2;
        if(check(mid))r=mid;
        else l=mid;
    }
    printf("%.6lf\n",(l+r)/2);
}