You have an unbalanced coin, the state of which (head or tail) after tossing depends on the previous state. If the head side is up now, the probability that the coin is still head up after tossing is $p$, and the probability that the tail side comes up is $1-p$. Similarly, if the tail side is up now, the probability that the coin is still tail up after tossing is $p$, and the probability that the head side comes up is $1-p$.
Assume the coin is head up initially. You will then toss it for $n$ times, during these $n$ times, every time the head side is up, you will get one score. What is the probability that you will get exactly $k$ score for $k=0,1,2,\cdots,n$?
intinit(int n){ int len = 1, p = 0; for (; len < n; len <<= 1, ++p) ; for (int i = 1; i < len; ++i) { rev[i] = rev[i >> 1] >> 1 | (i & 1) << p - 1; } for (int i = 0; i < len >> 1; ++i) { wn[i] = Complex(cos(2 * PI / len * i), sin(2 * PI / len * i)); } return len; }
voidfft(Complex *a, int len, bool inv = 0){ for (int i = 0; i < len; ++i) { if (i < rev[i]) { std::swap(a[i], a[rev[i]]); } } for (int i = 1; i < len; i <<= 1) { for (int j = 0; j < len; j += i << 1) { for (int k = 0; k < i; ++k) { Complex x = a[j + k], y = wn[len / (i << 1) * k] * a[j + i + k]; a[j + k] = x + y; a[j + i + k] = x - y; } } } if (inv) { std::reverse(a + 1, a + len); for (int i = 0; i < len; ++i) { a[i] = a[i] / len; } } }
structMatrix { int n[2][2]; double *a[2][2];
friend Matrix operator*(Matrix const &x, Matrix const &y) { int len = init(x.n[0][0] + y.n[0][0] - 1); for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { for (int k = 0; k < len; ++k) { A[i][j][k] = k < x.n[i][j] ? x.a[i][j][k] : 0; B[i][j][k] = k < y.n[i][j] ? y.a[i][j][k] : 0; C[i][j][k] = 0; } fft(A[i][j], len); fft(B[i][j], len); } } for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { for (int k = 0; k < 2; ++k) { for (int l = 0; l < len; ++l) { C[i][j][l] = C[i][j][l] + A[i][k][l] * B[k][j][l]; } } fft(C[i][j], len, 1); } } Matrix ret; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { ret.n[i][j] = 0; for (int k = 0; k < 2; ++k) { ret.n[i][j] = std::max(ret.n[i][j], x.n[i][k] + y.n[k][j] - 1); } ret.a[i][j] = newdouble[ret.n[i][j]]; for (int k = 0; k < ret.n[i][j]; ++k) { ret.a[i][j][k] = C[i][j][k].a; } } } return ret; } } base, ans;
intmain(){ int n; double p; scanf("%d%lf", &n, &p); base.n[0][0] = base.n[0][1] = 2; base.n[1][0] = base.n[1][1] = 1; base.a[0][0] = newdouble[2]; base.a[0][0][0] = 0; base.a[0][0][1] = p; base.a[0][1] = newdouble[2]; base.a[0][1][0] = 0; base.a[0][1][1] = 1 - p; base.a[1][0] = newdouble; base.a[1][0][0] = 1 - p; base.a[1][1] = newdouble; base.a[1][1][0] = p; ans.n[0][0] = ans.n[1][1] = 1; ans.a[0][0] = newdouble; ans.a[0][0][0] = 1; ans.a[1][1] = newdouble; ans.a[1][1][0] = 1; int b = n; for (; b;) { if (b & 1) { ans = ans * base; } if (b >>= 1) { base = base * base; } } for (int i = 0; i <= n; ++i) { printf("%.10f\n", ans.a[0][0][i] + ans.a[1][0][i]); } return0; }