Unbalanced Coin

题目描述

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

题意概述

有一枚初始正面朝上的硬币。每次投掷时,硬币有$p$的概率保持不变,$1-p$的概率翻面。求$n$次投掷后,正面朝上的总次数分别为$0,1,2,\cdots,n$的概率。

数据范围:$1 \le n \le 10^5$。

算法分析

令$f_{i,0/1,j}$表示$i$次投掷后,硬币正面/反面朝上,且正面朝上的总次数为$j$概率。容易想到转移方程:

$$\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}$$

这样直接转移的复杂度是$O(n^2)$,考虑如何优化。

令$f_{i,0/1}$的普通型生成函数为$F_{i,0/1}$,即$F_{i,0/1}=\sum_{j=0}^i f_{i,0/1,j}x^j$。则上述转移方程可表示为:

$$\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}$$

容易想到其矩阵形式:

$$\begin{equation} \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]\end{equation}$$

如此便可以用矩阵快速幂来优化了。其中多项式乘法需要用FFT计算。

代码

#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;
}

DNA Sequence

题目描述

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

题意概述

给定$m$个长度不超过$10$的DNA片段,问所有长度为$n$的DNA序列中有几个不包含所有给定的DNA片段。
数据范围:$0 \le m \le 10, \; 1 \le n \le 2 \times 10^9$。

算法分析

对$m$个DNA片段建立AC自动机。令$f_{i, j}$表示长度为$i$的DNA序列在自动机上跑到节点$j$的方案数。转移可以用一个矩阵来表示,因此可以用矩阵快速幂来加速转移。

代码

#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;
}

Tower

题意概述

有一个数列,$a_1=1, \; a_n=2a_2a_{n-1}-a_{n-2}$。给定$a_2$,求这个数列前$N$项的平方和模$M$。有$T$组数据。
数据范围:$1 \le T \le 10^5, \; 1 \le a_2, M \le 10^9, \; 2 \le N \le 10^9$。

算法分析

设$p=2a_2$,$S_n$表示前$n$项的平方和(不包括前两项)。
那么
$$
\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}
$$
用快速幂来加速递推,可以在$O(\log N)$的时间复杂度内完成。