bzoj2676

二分概率+矩乘+dp

也是二分概率，然后dp[i][j][k]表示当前到了i，有j条命，下一次的收益是k，然后矩乘转移，但是我自己的似乎wa了，抄了liu_runda的才行，具体不知道为什么

注释的是我自己写的，谁能告诉我哪里错了？

#include<bits/stdc++.h>
using namespace std;
const int N = 122;
int n, r, q, tot;
int id[N][N];
double S;
struct matrix {
    double a[N][N];
    matrix() { for(int i = 0; i <= tot; ++i) for(int j = 0; j <= tot; ++j) a[i][j] = 0.0; }
    matrix friend operator * (const matrix &a, const matrix &b) {
        matrix ret;
        for(int k = 0; k <= tot; ++k)
            for(int i = 0; i <= tot; ++i) if(a.a[i][k] >= 1e-15)
                for(int j = 0; j <= tot; ++j) 
                    ret.a[i][j] += a.a[i][k] * b.a[k][j];
        return ret;     
    }
    void set() {
        for(int i = 0; i <= tot; ++i) a[i][i] = 1.0;
    }
};
double calc(double p)
{
    matrix A, B;
    A.set();
    /*
    第i轮j条命这一次得k分 
    dp[i][j][k]=(dp[i+1][j + 1][k + 1] + k) * p
    dp[i][j][k] += dp[i + 1][j - 1][1] * (1.0 - p)
    矩阵乘法 
    */
    B.a[tot][tot] = 1.0;
    for(int i = 1; i <= q; ++i)
        for(int j = 1; j <= r; ++j)
        {
            if(i > 1) B.a[id[i - 1][1]][id[i][j]] = 1.0 - p;
            B.a[tot][id[i][j]] = p * (double)j;
            if(i < q && j < r) B.a[id[i + 1][j + 1]][id[i][j]] = p;
            else if(i < q) B.a[id[i + 1][j]][id[i][j]] = p;
            else if(j < r) B.a[id[i][j + 1]][id[i][j]] = p;
            else B.a[id[i][j]][id[i][j]] = p;
        } 
//  for(int j = 0; j <= q; ++j) 
//      for(int k = 1; k <= r; ++k) B.a[tot][id[j][k]] = p * (double)k;
//  for(int j = 0; j <= q; ++j)
//      for(int k = 1; k <= r; ++k) B.a[id[min(j + 1, q)][min(k + 1, r)]][id[j][k]] = p;
//  for(int j = 2; j <= q; ++j)
//      for(int k = 1; k <= r; ++k) B.a[id[j - 1][1]][id[j][k]] = 1.0 - p;           
    for(int i = n; i; i >>= 1, B = B * B) if(i & 1) A = A * B;
//  double ret = 0.0;
//  for(int i = 0; i <= tot; ++i)
//  {
//      for(int j = 0; j <= tot; ++j) printf("%.6f ", A.a[i][j]);        
//      puts("");   
//  } 
//  for(int i = 0; i < tot; ++i) ret += A.a[tot][i];
//  printf("p = %.6f ret = %.6f\n", p, ret);
    return A.a[tot][id[q][1]];   
}
int main()
{
    scanf("%d%d%d%lf", &n, &r, &q, &S);
    for(int i = 1; i <= q; ++i)
        for(int j = 1; j <= r; ++j) id[i][j] = tot++;
//  printf("tot = %d\n", tot);
    double l = 0.0, r = 1.0, ans = -1.0;
    while(r - l > 1e-10) 
    {
        double mid = (l + r) / 2.0;
        if(calc(mid) > S) r = ans = mid;
        else l = mid;
    }
    if(ans == -1.0) puts("Impossible.");
    else printf("%.6f\n", ans);
    return 0;
}