Description
偶然间,chnlich 发现了他小时候玩过的一个游戏“魂斗罗”,于是决定怀旧。但是这是一个奇怪的魂斗罗 MOD。
有 N 个关卡,初始有 Q 条命。
每通过一个关卡,会得到 u 分和1条命,生命上限为 Q。其中 u=min(最近一次连续通过的关数,R)。
若没有通过这个关卡,将会失去1条命,并进入下一个关卡。
当没有生命或没有未挑战过的关卡时,游戏结束,得到的分数为每关得到的分数的总和。
由于 chnlich 好久不玩这个游戏了,每条命通过每个关卡的概率均为p(0<=p<=1),原先 chnlich 的最高分纪录是 S。
现在 chnlich 想要知道,当 p 至少为多少时,chnlich 期望获得的总分数能够超过原先的最高分。
（本题题意可能有些难理解所以给了样例）
Input
输入共一行,分别表示整数 N,整数 R,整数 Q,原先的最高分整数 S。
Output
输出共一行,若不存在这样的 p,输出"Impossible."(不包含引号),否则输出 p(保留6位小数)。
Sample Input
【样例输入一】
4 2 1 5
【样例输入二】
12 3 2 12
Sample Output
【样例输出一】
0.880606
【样例输出二】
0.687201
Data Constraint
【数据说明】
对于(20%)的数据,(N<=15)
对于(50%)的数据,(N<=10000)
对于(100%)的数据,(N<=10^8,1<=R<=20,1<=Q<=5),保证(S)是一个可能出现的分数。
【补充说明】
例如,当 N=12,R=3,Q=2时
第一关:未通过 u=0 获得分数0 总分为0 剩余生命1
第二关:通过 获得分数1 总分为1 剩余生命2
第三关:通过 获得分数2 总分为3 剩余生命2
第四关:通过 获得分数3 总分为6 剩余生命2
第五关:通过 获得分数3 总分为9 剩余生命2
第六关:未通过 获得分数0 总分为9 剩余生命1
第七关:通过 获得分数1 总分为10 剩余生命2
第八关:未通过 获得分数0 总分为10 剩余生命1
第九关:未通过 获得分数0 总分为10 剩余生命0
游戏结束 总分为10
这是 chnlich 游戏的一种可能性
Solution
看完题后想到了(概率DP)，状态转移方程很容易得出，但是不好求答案。最后弃了。
正解是(二分答案+概率DP+矩乘优化)。
我们先二分一个(p)，然后对于这个(p)做(DP)。
设(f[i][j][k])表示玩到第(i)关，当前(u)为(j)，(k)连胜的概率。
那么我们很容易的到转移方程（胜和负）
而答案就是所有的(f[i][j][k]*p*min(j+1,r))的和，比较S再二分即可。
这样(O(log(1000000)*n*r*q))，50分到手。
然后我们考虑矩乘优化。
我们将(p*q)的数组编号，(1)~(len-1)
我们设f[i][j]表示从i状态转移到j状态的概率。
而len的那一列表示答案。先求出初始矩阵，然后转移即可。（详见标）
优化：我们可以发现，当(u>=q - 1)的时候，它的生命一定是满的，所以最多有大概42种状态。
那些多余状态不考虑的话，(len)将会大大减少。
记得注意精度问题！！！
Code
#include <cstdio>
#include <cstring>
#define db double
#define mem(a, x) memset(a, x, sizeof a)
#define fo(x, a, b) for (int x = a; x <= b; x++)
using namespace std;
int n, R, Q, S, len = 0, num[21][6];
db l = 0.0, r = 1.0, mid;
struct matrix
{
db a[45][45];
matrix operator *(const matrix &b)
{
matrix c;
mem(c.a, 0);
fo(k, 1, len)
fo(i, 1, len)
{
if (a[i][k] < 0.00000000001) continue;
fo(j, 1, len)
c.a[i][j] += a[i][k] * b.a[k][j];
}
return c;
}
}a, s, c;
inline int read()
{
int x = 0; char c = getchar();
while (c < '0' || c > '9') c = getchar();
while (c >= '0' && c <= '9') x = (x << 1) + (x << 3) + (c ^ 48), c = getchar();
return x;
}
inline int min(int x, int y) {return x < y ? x : y;}
bool check(db p)
{
mem(a.a, 0);
a.a[len][len] = 1.0;
fo(j, 0, R)
fo(k, 1, Q)
{
a.a[num[j][k]][len] = min(j + 1, R) * p;
a.a[num[j][k]][num[min(j + 1, R)][min(k + 1, Q)]] = p;
a.a[num[j][k]][num[0][k - 1]] = 1 - p;
}
// printf("%d
", len);
// fo(i, 1, len)
// {
// fo(j, 1, len)
// printf("%.2lf ", a.a[i][j]);
// puts("");
// }
int y = n;
mem(s.a, 0);
fo(i, 1, len) s.a[i][i] = 1;
while (y)
{
if (y & 1) s = s * a;
a = a * a; y >>= 1;
}
return s.a[num[0][Q]][len] > S;
}
int main()
{
freopen("contra.in", "r", stdin);
freopen("contra.out", "w", stdout);
n = read(), R = read(), Q = read(), S = read();
fo(j, 0, R)
fo(k, min(j + 1, Q), Q)
num[j][k] = ++len;
++len;
if (! check(1.0)) {printf("Impossible."); return 0;}
while (l + 0.000000001 < r)
{
mid = (l + r) / 2.0;
if (check(mid)) r = mid;
else l = mid;
}
printf("%.6lf", r);
return 0;
}