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\n", 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;
}