## 题目描述

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$?

## 算法分析

\begin{align} f_{i,0,j} &= f_{i-1,0,j-1} \times p+f_{i-1,1,j-1} \times (1-p) \\ f_{i,1,j} &= f_{i-1,0,j} \times (1-p)+f_{i-1,1,j} \times p \end{align}

\begin{align} F_{i,0} &= F_{i-1,0} \times px+F_{i-1,1} \times (1-p)x \\ F_{i,1} &= F_{i-1,0} \times (1-p)+F_{i-1,1} \times p \end{align}

$$$$\left[\begin{matrix} px & (1-p)x \\ 1-p & p\end{matrix}\right] \left[\begin{matrix} F_{i-1,0} \\ F_{i-1,1} \end{matrix}\right]= \left[\begin{matrix} F_{i,0} \\ F_{i,1} \end{matrix}\right]$$$$

## 代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

int const N = 200005;
double const PI = acos(-1);

int rev[N];

struct Complex {
double a, b;
Complex(double _a = 0, double _b = 0) : a(_a), b(_b) {}
friend Complex operator+(Complex const &x, Complex const &y) {
return Complex(x.a + y.a, x.b + y.b);
}
friend Complex operator-(Complex const &x, Complex const &y) {
return Complex(x.a - y.a, x.b - y.b);
}
friend Complex operator*(Complex const &x, Complex const &y) {
return Complex(x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a);
}
friend Complex operator/(Complex const &x, double const &y) {
return Complex(x.a / y, x.b / y);
}
} wn[N], A[2][2][N], B[2][2][N], C[2][2][N];

int init(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;
}

void fft(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;
}
}
}

struct Matrix {
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] = new double[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;

int main() {
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] = new double[2];
base.a[0][0][0] = 0;
base.a[0][0][1] = p;
base.a[0][1] = new double[2];
base.a[0][1][0] = 0;
base.a[0][1][1] = 1 - p;
base.a[1][0] = new double;
base.a[1][0][0] = 1 - p;
base.a[1][1] = new double;
base.a[1][1][0] = p;
ans.n[0][0] = ans.n[1][1] = 1;
ans.a[0][0] = new double;
ans.a[0][0][0] = 1;
ans.a[1][1] = new double;
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]);
}
return 0;
}

## 题目描述

It’s well known that DNA Sequence is a sequence only contains A, C, T and G, and it’s very useful to analyze a segment of DNA Sequence，For example, if a animal’s DNA sequence contains segment ATC then it may mean that the animal may have a genetic disease. Until now scientists have found several those segments, the problem is how many kinds of DNA sequences of a species don’t contain those segments.
Suppose that DNA sequences of a species is a sequence that consist of A, C, T and G, and the length of sequences is a given integer $n$.

## 代码

#include <cstdio>
#include <cstring>

static const int N = 11;
static const int POOL = 105;
char s[N];

int get_id(char c) {
return c == 'A' ? 0 : c == 'C' ? 1 : c == 'T' ? 2 : 3;
}

typedef int Array[POOL][POOL];

int numn = 1, que[POOL];
Array mp, a, tmp;
struct ACAutomaton {
int nxt[4], fail;
bool match;
} node[POOL];

void insert(char s[]) {
int root = 1, len = strlen(s);
for (int i = 0; i < len; ++ i) {
int p = get_id(s[i]);
if (! node[root].nxt[p]) node[root].nxt[p] = ++ numn;
root = node[root].nxt[p];
}
node[root].match = true;
}

void build() {
int qb = 0, qe = 1; que[0] = 1;
while (qb < qe) {
int u = que[qb ++];
for (int i = 0; i < 4; ++ i)
if (node[u].nxt[i]) {
node[node[u].nxt[i]].fail = 1;
if (u > 1) {
int root = node[u].fail;
while (root > 1 && ! node[root].nxt[i]) root = node[root].fail;
if (node[root].nxt[i]) {
node[node[u].nxt[i]].fail = node[root].nxt[i];
if (node[node[root].nxt[i]].match) node[node[u].nxt[i]].match = true;
}
}
que[qe ++] = node[u].nxt[i];
}
}
}

void mul(Array a, Array b) {
memset(tmp, 0, sizeof tmp);
for (int i = 1; i <= numn; ++ i)
for (int j = 1; j <= numn; ++ j)
if (a[i][j]) for (int k = 1; k <= numn; ++ k) (tmp[i][k] += 1ll * a[i][j] * b[j][k] % 100000) %= 100000;
memcpy(a, tmp, sizeof tmp);
}

void power(Array a, int b) {
while (b) {
if (b & 1) mul(a, mp);
mul(mp, mp), b >>= 1;
}
}

int main() {
int m, n; scanf("%d%d", &m, &n);
for (int i = 0; i < m; ++ i) scanf(" %s", s), insert(s);
build();
for (int i = 1; i <= numn; ++ i) if (! node[i].match)
for (int j = 0; j < 4; ++ j) {
int root = i;
while (root > 1 && ! node[root].nxt[j]) root = node[root].fail;
if (node[root].nxt[j]) ++ mp[i][node[root].nxt[j]]; else ++ mp[i][root];
}
for (int i = 1; i <= numn; ++ i) a[i][i] = 1;
power(a, n);
int ans = 0;
for (int i = 1; i <= numn; ++ i) if (! node[i].match) (ans += a[1][i]) %= 100000;
printf("%d\n", ans);
return 0;
}


## 算法分析

\begin{align} & a_n^2=(pa_{n-1}-a_{n-2})^2=p^2a_{n-1}^2-2pa_{n-1}a_{n-2}+a_{n-2}^2 \\ & S_n=S_{n-1}+a_n^2 \\ & a_na_{n-1}=a_{n-1}(pa_{n-1}-a_{n-2})=pa_{n-1}^2-a_{n-1}a_{n-2} \end{align}

$$\begin{bmatrix} p^2 & 1 & -2p & 0 \\ 1 & 0 & 0 & 0 \\ p & 0 & -1 & 0 \\ 1 & 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} a_n^2 \\ a_{n-1}^2 \\ a_na_{n-1} \\ S_{n-1} \end{bmatrix}= \begin{bmatrix} a_{n+1}^2 \\ a_n^2 \\ a_{n+1}a_n \\ S_n \end{bmatrix}$$