BZOJ 3213 KMP Expectation

Calculate expectation with $fail$.

Description

Given a probability $p$ of throwing a coin and the $H$ surface is on the top. Accordingly, $1-p$ of $T$ surface.

Now given a target sequence consisting of $H$ and $T$, calculate how many times (in expectation) you need to make the result as such a sequence.

Solution

The target sequence must be the suffix of our result, so must a substring. We can evaluate how our result transforms when trying to fit the target sequence.

Equation is :

$f(i) = f(i-1) + 1 + (1-p_i)(f(i) - f(file_i))$.

Swap items to get expression of $f(i)$.

Need to implement high-precision number class.

#include <bits/stdc++.h>

using namespace std;

const int maxn = 1005;

typedef long long ll;

int n, a, b;
int nxt[maxn], tow[maxn];
char s[maxn];

struct bigint
{
  static const int width = 9;
  static const ll base = 1000000000;
  ll a[1005]; int len;
  bigint() {
    memset(a, 0, sizeof a); len = 1;
  }
  ll& operator [] (int x) {
    return a[x];
  }
  bigint(int x) {
    len = 0; 
    if (!x) len = 1;
    else {
      while (x) a[++len] = x % base, x /= base;
    }
  }
  bigint operator + (bigint b) {
    bigint ret; ret.len = max(len, b.len);
    for (int i = 1; i <= ret.len; ++i) {
      ret[i] = a[i] + b[i];
      //ret[i + 1] += ret[i] / base;
      //ret[i] %= base;
    }
    
    for (int i = 1; i <= ret.len; ++i)
      ret[i + 1] += ret[i] / base, ret[i] %= base;
    for (; ret[ret.len + 1]; ret.len++) ret[ret.len + 2] += ret[ret.len + 1] / base, ret[ret.len + 1] %= base;
    /*while (ret[ret.len] >= base) {
      ret[ret.len + 1] += ret[ret.len] / base;
      ret[ret.len] %= base;
      ret.len++;
    }*/
    return ret;
  }
  bigint operator * (bigint b) {
    bigint ret;
    for (int i = 1; i <= len; ++i)
      for (int j = 1; j <= b.len; ++j)
        ret[i + j - 1] += a[i] * b[j], ret[i + j] += ret[i + j - 1] / base, ret[i + j - 1] %= base;
    ret.len = len + b.len - 1;
    /*while (ret[ret.len] >= base) {
      ret[ret.len + 1] += ret[ret.len] / base;
      ret[ret.len] %= base;
      ret.len++;
    }*/
    for (; ret[ret.len + 1]; ret.len++) ret[ret.len + 2] += ret[ret.len + 1] / base, ret[ret.len + 1] %= base;
    return ret;
  }
  bigint operator / (int b) {
    bigint ret; ret.len = len; ll bit = 0;
    for (int i = len; i; --i) {
      ret[i] = (a[i] + bit * base) / b;
      bit = (bit * base + a[i]) % b;
    }
    while (!ret[ret.len]) ret.len--;
    return ret;
  }
  bigint operator - (bigint b) {
    bigint ret; ret.len = max(len, b.len);
    for (int i = 1; i <= ret.len; ++i) ret[i] = a[i] - b[i];
    
    for (int i = 1; i <= ret.len; ++i) while (ret[i] < 0) ret[i + 1]--, ret[i] += base;
    while (!ret[ret.len]) ret.len--;
    return ret;
  }
  bool operator == (bigint b) {
    if (len != b.len) return false;
    for (int i = 1; i <= len; ++i) if (a[i] != b[i]) return false;
    return true;
  }
  friend bool operator % (bigint a, int b) {
    /*if (b == 4 || b == 2 || b == 3 || b == 6) {
      (a / b * bigint(b)).print(); putchar('\n');
      a.print(); putchar('\n');
      printf("%d\n", b);
    }*/
    return a / b * bigint(b) == a;
  }
  inline void print() {
    printf("%lld", a[len]);
    for (int i = len - 1; i; --i) 
      printf("%09lld", a[i]);
  }
};
struct frac
{
  bigint up, dn;
  frac() {
    //up = bigint(0); dn = bigint(0);
  }
  frac(int u, int d) : up(u), dn(d) {}
  frac(bigint u, bigint d) : up(u), dn(d) {}
  frac operator + (frac b) {
    frac ret = frac(up * b.dn + b.up * dn, dn * b.dn);
    //ret.print();
    //ret.cancel();
    return ret;
  }
  frac operator - (frac b) {
    frac ret = frac(up * b.dn - b.up * dn, dn * b.dn);
    //ret.print();
    //ret.cancel();
    return ret;
  }
  frac operator * (frac b) {
    //if (up * b.up == bigint(0)) return frac(0, 1);
    frac ret = frac(up * b.up, dn * b.dn);
    //ret.print();
    //ret.cancel();
    return ret;
  }
  frac operator / (frac b) {
    frac ret = frac(up * b.dn, dn * b.up);
    //ret.print();
    //ret.cancel();
    return ret;
  }
  inline void cancel() {
    //this->print();
    for (int i = 2; i <= 100; ++i) {
      
      while (up % i && dn % i) {
        up = up / i, dn = dn / i;
      }
    }
  }
  void print() {
    up.print(); putchar('/'); dn.print();
    putchar('\n');
  }
};
frac p[2], f[maxn];

void getnxt() {
  int j = 0;
  nxt[0] = 0;
  tow[1] = 0;
  s[0] = '%';
  for (int i = 1; i <= n; ++i) {
    j = nxt[i - 1];
    while (j && s[j + 1] != s[i]) j = nxt[j];
    if (i > 1 && s[j + 1] == s[i]) nxt[i] = j + 1, j++;
    else nxt[i] = 0;
    if (i < n) {
      while (j && s[j + 1] == s[i + 1]) j = nxt[j];
      if (s[j + 1] != s[i + 1]) tow[i + 1] = j + 1;
      else tow[i + 1] = 0;
    }
  }
}

int main() {
  //freopen("3213.out", "w", stdout);
  scanf("%d%d", &a, &b);
  p[1] = frac(a, b); p[0] = frac(b - a, b);
  //p[1].cancel();
  //p[0].cancel();
  scanf("%s", s + 1);
  n = strlen(s + 1);
  getnxt();
  f[0] = frac(0, 1);
  //frac t1 = frac(2147483646, 2147483645);
  //frac t2 = t1;
  //frac t3 = t2 * t1;
  //t3.print();
  /*frac t1 = frac(2, 3), t2 = frac(1, 3);
  (t1 + t2).print();
  frac t3 = t1 + t2; t3.cancel(); t3.print();
  t3 = t1 / t2; t3.cancel(); t3.print();
  t3 = t1 - t2; t3.cancel(); t3.print();*/
  for (int i = 1; i <= n; ++i) {
    int typ = (s[i] == 'H');
    //printf("%d\n", tow[i]);
    //f[tow[i]].print();
    f[i] = (f[i - 1] - (p[typ ^ 1] * f[tow[i]]) + frac(1, 1)) / p[typ];
    //f[i].print();
    f[i].cancel();
    //f[i].print();
  }
  f[n].print();
  return 0;
}