[LibreOJ 6436][PKUSC2018]神仙的游戏

考场上写炸了FFT的zzs事野蛮人(确信)

兄啊这题怎么还卡常啊!

不妨设\(n = |S|\)。那么如果串有一个长度为$i$的border,就意味着串有一个长度为\(n - i\)的循环节。

那么考虑有两个非?的位置\(i\)和\(j\)(不妨设\(i < j\)),且\(S_i\neq S_j\)。那么我们设\(l = j - i\),则\(l\)的约数都不能成为循环节的长度,否则会出现矛盾(一个位置又要是0又要是1)。

那么考虑对于所有可能的\(l\)判定是否存在这种不合法的\(i\)与\(j\)。我们利用编辑距离函数:

\[f_l = \sum_{i} A_i A_{i + l}(A_i - A_{i + l})^2\]

这里\(A_i\)相当于\(S_i\)的一个重编码,\(S_i\)为?时应当为0,其他情况一种为1一种为2。这个函数不方便我们卷积,所以我们考虑把\(A\)翻转一下(新的序列称为\(B\)):

\[
\begin{aligned}
f_l &= \sum_{i} B_{n - i - 1} A_{i + l}(B_{n - i - 1} A_{i + l})^2\\
&= \sum_{i} B_{n - i - 1}^3 A_{i + l} - 2B_{n - i - 1}^2 A_{i + l}^2 + B_{n - i - 1}A_{i + l}^3
\end{aligned}\]

然后问题变成了三个卷积。FFT一下就行了(请注意IDFT只需要一次……否则会被卡常)。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cmath>
#include <complex>
#include <vector>
using R = double;
struct C {
  R x, y;
  C(R xx = 0.00, R yy = 0.00) {
    x = xx; y = yy;
  }
  R real() {
    return x;
  }
};
C operator +(const C &a, const C &b) {
  return C(a.x + b.x, a.y + b.y);
}
C operator -(const C &a, const C &b) {
  return C(a.x - b.x, a.y - b.y);
}
C operator *(const C &a, const C &b) {
  return C(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

using ll = long long;
constexpr R pi = acos(-1);
int flip(int bi, int x) {
  int ret = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ret += (1 << (bi - i - 1));
    }
  }
  return ret;
}
void FFT(C *A, int bi, R flag = 1.00) {
  int n = 1 << bi;
  for(int i = 0; i < n; i ++) {
    int v = flip(bi, i);
    if(i < v) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    R rd = pi / (R(L));
    C xi_n(cos(rd), sin(flag * rd));
    for(int j = 0; j < n; j += (L << 1)) {
      C xi(1.00, 0.00);
      for(int i = j; i < j + L; i ++) {
        C x = A[i], y = A[i + L];
        A[i] = x + xi * y;
        A[i + L] = x - xi * y;
        xi = xi * xi_n;
      }
    }
  }
}

const int maxn = 500005;
R tot[maxn];
C T[maxn << 2], TT[maxn << 2], D[maxn << 2];

char S[maxn];
inline int conv(char c) {
  if(c == '1') return 2;
  if(c == '0') return 1;
  return 0;
}
int A[maxn], B[maxn];
bool boom[maxn];
int main() {
  scanf("%s", S); int n = strlen(S);
  for(int i = 0; i < n; i ++) {
    A[i] = B[n - 1 - i] = conv(S[i]);
  }
  
  int len = 1, bi = 0;
  while(len < (2 * n)) {
    len <<= 1; bi ++;
  }
  
  C z(0.00, 0.00);
  // std::fill(T, T + len, z);
  // std::fill(TT, TT + len, z);
  for(int i = 0; i < n; i ++) {
    T[i] = A[i] * A[i] * A[i]; TT[i] = B[i];
  }
  FFT(T, bi); FFT(TT, bi);
  for(int i = 0; i < len; i ++) {
    D[i] = D[i] + (T[i] * TT[i]);
  }
  
  std::fill(T, T + len, z);
  std::fill(TT, TT + len, z);
  for(int i = 0; i < n; i ++) {
    T[i] = A[i] * A[i]; TT[i] = B[i] * B[i];
  }
  FFT(T, bi); FFT(TT, bi);
  for(int i = 0; i < len; i ++) {
    D[i] = D[i] + C(-2.00, 0.00) * (T[i] * TT[i]);
  }
  
  std::fill(T, T + len, z);
  std::fill(TT, TT + len, z);
  for(int i = 0; i < n; i ++) {
    T[i] = A[i]; TT[i] = B[i] * B[i] * B[i];
  }
  FFT(T, bi); FFT(TT, bi);
  for(int i = 0; i < len; i ++) {
    D[i] = D[i] + (T[i] * TT[i]);
  }
  FFT(D, bi, -1.00);
  
  for(int i = 1; i < n; i ++) {
    if(fabs(D[n - 1 + i].real() / (R(len))) > 0.99) {
      boom[i] = true;
    }
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = i * 2; j <= n; j += i) {
      boom[i] = (boom[i] || boom[j]);
    }
  }
  ll ret = 0LL;
  for(int i = 1; i <= n; i ++) {
    if(!boom[n - i]) {
      ret ^= (ll(i)) * (ll(i));
    }
  }
  printf("%lld\n", ret);
  return 0;
}

[BZOJ 4259]残缺的字符串

终于肝掉了这道提……

考虑用编辑距离函数处理这个问题。考虑每个可能的结尾\(i\),定义\(f_i\)(这里把\(*\)看做\(0\)把):

\[f_i = \sum_{j = 0}^{m - 1} A_j B_{i - m + 1 + j}\,(A_j - B_{i - m + 1 + j})^2\]

考虑翻转\(A\),式子变成:

\[f_i = \sum_{j = 0}^{m - 1} A_{m - 1 - j}\,B_{i - m + 1 + j}\,(A_j - B_{i - m + 1 + j})^2\]

式子给人一种卷积的即视感……再化简一下,得到:

\[f_i = \sum_{j = 0}^{m - 1} A_{m - 1 - j}^3\,B_{i - m + 1 + j}\,- 2A_{m - 1 - j}^2\,B_{i - m + 1 + j}^2\,+ A_{m - 1 - j}\,B_{i - m + 1 + j}^3\]

这是三个卷积……跑三遍FFT即可得到\(f_i\),当且仅当\(f_i = 0\)时,以\(i\)结尾的长度为\(m\)的子串才能匹配。

顺便提个几个坑……这个题用long double + std::complex会因为常数巨大而T掉。并且注意eps不要调那么细啊……这个东西精度误差极大,我直接调成1了。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <algorithm>
#include <utility>
#include <queue>
#include <complex>
constexpr int maxn = 300005;
using R = double;
using C = std::complex<R>;
constexpr R pi = acos(-1);
inline int flip(int bi, int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ans += (1 << (bi - i - 1));
    }
  }
  return ans;
}
inline void FFT(C *A, int bi, int n, R flag = 1.00) {
  for(int i = 0; i < n; i ++) {
    int fl = flip(bi, i);
    if(i < fl) std::swap(A[fl], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    R rd = pi / (R(L));
    C xi(cos(rd), sin(flag * rd));
    for(int i = 0; i < n; i += (L << 1)) {
      C w(1.00, 0.00);
      for(int j = i; j < i + L; j ++) {
        C v1 = A[j], v2 = A[j + L];
        A[j] = v1 + v2 * w;
        A[j + L] = v1 - v2 * w;
        w = w * xi;
      }
    }
  }
}

constexpr R eps = 1.00;
inline int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0.00) {
      return -1;
    } else {
      return 1;
    }
  }
}
int main() {
  static char A[maxn], B[maxn];
  static int V1[maxn << 2], V2[maxn << 2];
  static C P1[maxn << 2], P2[maxn << 2], P3[maxn << 2];
  static R F[maxn << 2];
  int m, n; scanf("%d%d", &m, &n);
  if(m > n) {
    puts("0"); return 0;
  }
  scanf("%s%s", A, B);
  std::reverse(A, A + m);
  for(int i = 0; i < m; i ++) {
    V1[i] = (A[i] == '*') ? 0 : (A[i] - 'a' + 1);
  }
  for(int i = 0; i < n; i ++) {
    V2[i] = (B[i] == '*') ? 0 : (B[i] - 'a' + 1);
  }
  int sz = 1, bi = 0;
  while(sz <= (m + n)) {
    sz *= 2; bi ++;
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i] * V1[i] * V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += P3[i].real() / (R(sz));
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i] * V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i] * V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += (R(-2)) * P3[i].real() / (R(sz));
  }
  std::fill(P1, P1 + sz, C(0, 0));
  std::fill(P2, P2 + sz, C(0, 0));
  for(int i = 0; i < m; i ++) {
    P1[i] = V1[i];
  }
  for(int i = 0; i < n; i ++) {
    P2[i] = V2[i] * V2[i] * V2[i];
  }
  FFT(P1, bi, sz); FFT(P2, bi, sz);
  for(int i = 0; i < sz; i ++) {
    P3[i] = P1[i] * P2[i];
  }
  FFT(P3, bi, sz, -1.00);
  for(int i = 0; i < sz; i ++) {
    F[i] += P3[i].real() / (R(sz));
  }
  static int ans[maxn]; int cnt = 0;
  for(int i = 0; i < n; i ++) {
#ifdef LOCAL
    printf("F[%d] : %.10Lf\n", i, F[i]);
#endif
    if(sign(F[i]) == 0 && i >= m - 1) {
      ans[cnt ++] = i - m + 1 + 1;
    }
  }
  printf("%d\n", cnt);
  for(int i = 0; i < cnt; i ++) {
    if(i) putchar(' ');
    printf("%d", ans[i]);
  }
  putchar('\n');
  return 0;
}

[BZOJ 3527][ZJOI2014]力

现在才明白卷积真的是the deep, dark fantasies……(逃

首先约去\(q_i\),得到:

\[E_j = \sum_{i < j}\frac{q_i}{(j - i)^2} - \sum_{j < i}\frac{q_i}{(j - i)^2}\]

注意到如果很容易得到等式前半部分的高效求法,后半部分把序列反过来就能做了。

那么我们会发现,设\(g_x = \frac{1}{x^2}\),然后式子前半部分(姑且称作\(p_j\))可表示为:

\[p_j = \sum_{i < j} q_i g_{j - i}\]

这不就是个卷积吗?FFT一波即可。

代码:

/**************************************************************
    Problem: 3527
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9984 ms
    Memory:28952 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <utility>
#include <queue>
#include <cassert>
#include <complex>
typedef long double R;
typedef std::complex<R> C;
const R eps = 1e-3;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0) {
      return -1;
    } else {
      return 1;
    }
  }
}
const int maxn = 400005;
const double pi = acos(-1);
int bi, len;
int flip(int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ans += 1 << (bi - i - 1);
    }
  }
  return ans;
}
void fft(C *A, double g = 1) {
  for(int i = 0; i < len; i ++) {
    int R = flip(i);
#ifdef LOCAL
    // printf("The flipping of %d is %d\n", i, R);
#endif
    if(i < R) {
      std::swap(A[R], A[i]);
    }
  }
  for(int L = 1; L < len; L <<= 1) {
    C xi_n(cos(pi / (double(L))), sin(g * pi / (double(L))));
    for(int i = 0; i < len; i += L << 1) {
      C xi(1, 0);
      for(int j = i; j < i + L; j ++) {
        C a = A[j], b = A[j + L];
        A[j] = a + xi * b;
        A[j + L] = a - xi * b;
        xi = xi * xi_n;
      }
    }
  }
}
 
int main() {
  int n;
  static C A[maxn], rd[maxn]; static R ans[maxn];
  static R q[maxn];
  scanf("%d", &n);
  bi = 0; len = 1;
  while(len <= 2 * n) {
    len <<= 1;
    bi ++;
  }
  for(int i = 0; i < n; i ++) {
    scanf("%Lf", &q[i]); A[i] = q[i];
  }
  assert(sign(rd[0].real()) == 0);
  for(int i = 1; i < n; i ++) {
    rd[i] = 1.00 / ((R(i)) * (R(i)));
#ifdef LOCAL
    printf("rd[%d] : %.5Lf\n", i, rd[i].real());
#endif
  }
  fft(A); fft(rd);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] += A[i].real() / len;
#ifdef LOCAL
    printf("delta_v of %d : %.5Lf\n", i, A[i].real() / len);
#endif
  }
  std::reverse(q, q + n);
  for(int i = 0; i < len; i ++) A[i] = 0;
  for(int i = 0; i < n; i ++) {
    A[i] = q[i];
  }
  fft(A);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] -= A[n - 1 - i].real() / len;
    printf("%.3Lf\n", ans[i]);
  }
  return 0;
}