[BZOJ 3456]城市规划

aji多项式求逆毁我青春,,,

设\(f_i\)表示\(i\)个点的有标号无向联通图,考虑所有可能的图(记\(F_i\)为\(i\)个点的有标号无向图,显然\(F_i = 2^{\binom{i}{2}}\))和\(f\)的关系(使用图计数的经典套路:枚举1所在的联通块大小):

\[F_n = \sum_{i = 1}^n \binom{n-1}{i-1} f_i F_{n - i}\]

看起来事卷积?但是这个卷积没有办法直接用FFT/NTT求(当然分离一下项啊,移下项就可以分治NTT力)。

考虑进一步化简柿子。完全展开后会发现右边有个非常碍眼的\((n-1)!\),所以两边除一下:

\[\frac{F_n}{(n-1)!} =\sum_{i = 1}^n \frac{f_i}{(i-1)!}\cdot \frac{F_{n - i}}{(n - i)!}\]

然后这个问题就很毒瘤了:我们要求的答案多项式(除上那个阶乘)和一个已知多项式做卷积,可以得到另一个已知多项式……这样就需要多项式除法了,于是乎多项式逆元派上了用场。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const int maxn = 1020010;
const ll ha = 1004535809LL;
const ll bs = 3LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % ha;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x,  ha - 2LL);
}
 
int flip(int bi, int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if((1 << i) & x) {
      ans += (1 << (bi - i - 1));
    }
  }
  return ans;
}
void NTT(ll *A, int bi, bool flag = false) {
  int n = 1 << bi;
  for(int i = 0; i < n; i ++) {
    int v = flip(bi, i);
    if(v < i) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    ll xi = pow_mod(3LL, (ha - 1LL) / (ll(L << 1)));
    if(flag) xi = inv(xi);
    for(int i = 0; i < n; i += (L << 1)) {
      ll w = 1LL;
      for(int j = i; j < i + L; j ++) {
        ll v1 = A[j], v2 = A[j + L];
        A[j] = (v1 + (w * v2) % ha) % ha;
        A[j + L] = (v1 - (w * v2) % ha + ha) % ha;
        w = (w * xi) % ha;
      }
    }
  }
}
void poly_mul(ll *A, ll *B, int bi, ll *C) {
  static ll T1[maxn], T2[maxn];
  int n = (1 << bi);
  std::copy(A, A + n, T1);
  std::copy(B, B + n, T2);
  NTT(T1, bi); NTT(T2, bi);
#ifdef LOCAL
  puts("poly_mul :");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", A[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", B[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T2[i]);
  }
  puts("");
#endif 
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * T2[i]) % ha;
  }
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
#endif 
  NTT(T1, bi, true);
  ll inv_n = inv(n);
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * inv_n) % ha;
  }
  std::copy(T1, T1 + n, C);
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", C[i]);
  }
  puts("");
#endif 
}
 
void poly_inv(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = inv(B[0]);
  } else {
    poly_inv((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) + 1)) {
      bi ++; sz <<= 1;
    }
    ll inv_sz = inv(sz);
    static ll tmp[maxn];
    std::copy(B, B + mod, tmp);
    std::fill(tmp + mod, tmp + sz, 0LL);
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * BB[i]) % ha;
      tmp[i] = (tmp[i] * (ha - 1LL)) % ha;
      tmp[i] = (tmp[i] + 2LL) % ha;
      tmp[i] = (tmp[i] * BB[i]) % ha;
    }
    NTT(tmp, bi, true);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * inv_sz) % ha;
    }
    std::copy(tmp, tmp + mod, BB);
    std::fill(BB + mod, BB + sz, 0LL);
  }
}
 
int main() {
  static ll fac[maxn], A[maxn], B[maxn], BB[maxn];
  int n; scanf("%d", &n);
  int bi = 0, sz = 1;
  while(sz <= n + 1) {
    bi ++; sz <<= 1;
  }
  fac[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    fac[i] = (fac[i - 1] * (ll(i))) % ha;
  }
  B[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    B[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    B[i] = (B[i] * inv(fac[i])) % ha;
  }
  for(int i = 1; i <= n; i ++) {
    A[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    A[i] = (A[i] * inv(fac[i - 1])) % ha;
  }
  poly_inv(n + 1, B, BB);
  poly_mul(A, BB, bi, A);
  printf("%lld\n", (A[n] * fac[n - 1]) % ha);
  return 0;
}

[BZOJ 3328]PYXFIB

又学了个新套路呢qwq

题面要求的其实是:

\[\sum_{i = 0}^n [k|i]\binom{n}{i} F_i\]

不考虑那个\([k|i]\),式子是非常像一个二项式展开的……

考虑构造矩阵的二项式展开,可以发现(其中\(A\)是斐波那契数列的二项式展开):

\[(I + A)^n = \sum_{i = 0}^n \binom{n}{i} A^i\]

然后\(k=1\)的情况我们就会做辣!接下来的问题是,如何取出一个生成函数中次数为\(k\)倍数的项求和?

然后我们发现单位根有个非常优良的性质:

\[\frac{1}{k} \sum_{i=1}^{k} \xi_k^{ni} = [k|n]\]

这里的\(n\)是一个常数,但是其实可以对应到原式里的某一项的次数。至于证明……满足\(k|n\)的情况左边显然是一堆1取平均值,其他情况可以通过等比数列求和证明左边为0。

于是可以构造多项式\((I+xA)^n\),把所有\(k\)次单位根做为\(x\)(求这个的话,可以利用原根)带入,对结果取平均值即可。

代码:

/**************************************************************
    Problem: 3328
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9080 ms
    Memory:832 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cmath>
#include <vector>
typedef long long ll;
typedef ll Mat[2][2];
ll p; int k;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % p;
  while(b) {
    if(1LL & b) ans = (ans * res) % p;
    res = (res * res) % p;
    b /= 2LL;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, p - 2LL);
}
 
void factor(int x, std::vector<int> &V) {
  int bd = sqrt(x + 0.5);
  for(int i = 2; i <= bd; i ++) {
    if(x % i == 0) {
      V.push_back(i);
      while(x % i == 0) x /= i;
    }
  }
  if(x > 1) V.push_back(x);
}
int get_phi() {
  if(p == 2LL) return 1LL;
  std::vector<int> V;
  factor(p - 1LL, V);
  for(int i = 2; i <= (p - 1); i ++) {
    bool flag = true;
    for(int j = 0; j < V.size(); j ++) {
      int up = (p - 1) / V[j];
      if(pow_mod(i, up) == 1LL) {
        flag = false; break;
      }
    }
#ifdef LOCAL
    if(flag) printf("xi : %d\n", i);
#endif
    if(flag) return i;
  }
}
 
void mat_mul(const Mat &A, const Mat &B, int n, int m, int t, Mat &C) {
  Mat D; memset(D, 0, sizeof(D));
  for(int i = 0; i < n; i ++) {
    for(int j = 0; j < t; j ++) {
      for(int k = 0; k < m; k ++) {
        D[i][j] = (D[i][j] + (A[i][k] * B[k][j]) % p) % p;
      }
    }
  }
  memcpy(C, D, sizeof(D));
}
void mat_pow(const Mat &A, int n, ll b, Mat &ret) {
  Mat C, res;
  memset(C, 0, sizeof(C));
  for(int i = 0; i < n; i ++) C[i][i] = 1LL;
  memset(res, 0, sizeof(res));
  memcpy(res, A, sizeof(res));
  while(b) {
    if(1LL & b) mat_mul(C, res, n, n, n, C);
    mat_mul(res, res, n, n, n, res);
    b /= 2LL;
  }
  memcpy(ret, C, sizeof(C));
}
 
int main() {
  int T; scanf("%d", &T);
  while(T --) {
    ll n;
    scanf("%lld%d%lld", &n, &k, &p);
    ll xi = pow_mod(get_phi(), (p - 1) / k);
    ll w = 1LL;
    ll ans = 0;
    for(int i = 0; i < k; i ++) {
      Mat X;
      memset(X, 0, sizeof(X));
      X[0][0] = X[1][0] = X[0][1] = w;
      X[0][0] = (X[0][0] + 1LL) % p;
      X[1][1] = (X[1][1] + 1LL) % p;
      mat_pow(X, 2, n, X);
#ifdef LOCAL
      puts("X : ");
      for(int i = 0; i < 2; i ++) {
        for(int j = 0; j < 2; j ++) {
          printf("%lld ", X[i][j]);
        }
        puts("");
      }
#endif
      ans = (ans + X[0][0]) % p;
      w = (w * xi) % p;
    }
    ans = (ans * inv(k)) % p;
    printf("%lld\n", ans);
  }
  return 0;
}

[BZOJ 4919][Lydsy1706月赛]大根堆

ao神啊这题,,,

很显然的思路是设计状态,然后线段树合并或者线段树启发式合并转移来优化一下,但是真的不好写……

考虑求LIS的那个二分+单调数组来做,那个东西的本质其实就是维护了一堆各种长度的断点一样的东西(考虑答案为\(n\),\(n - 1\)……的情况,在值的选取上各个情况之间会有断点)。那么我们就用一个multiset来维护断点,然后启发式合并即可。

考虑一棵子树,把根的子树进行合并,根加进来之后对断点会有何影响?那么考虑现在那个multiset里那个东西的后继(要lower_bound,因为如果有和根权值一样的点的话拐点不会发生变化。这里讨论有后继的情况),如果要取到那个点的长度的LIS的话,现在只需要小于等于根的权值的点就行了,取到这个点也不会使LIS变大(它不能和根共存),所以那个点就不是拐点了。并且无论是否出现这种情况,根那个点都会成为一个新的拐点。

然后DFS一遍就行啦。

代码:

/**************************************************************
    Problem: 4919
    User: danihao123
    Language: C++
    Result: Accepted
    Time:992 ms
    Memory:25648 kb
****************************************************************/
 
#include <cstdio>
#include <set>
#include <vector>
const int maxn = 200005;
std::vector<int> G[maxn];
 
std::multiset<int> *st[maxn];
void merge(int x, int y) {
  if(st[x] -> size() < st[y] -> size()) std::swap(st[x], st[y]);
  while(st[y] -> size() > 0) {
    std::multiset<int>::iterator it = st[y] -> begin();
    int v = *it; st[y] -> erase(it);
    st[x] -> insert(v);
  }
}
int d[maxn];
void dfs(int x) {
  for(int i = 0; i < G[x].size(); i ++) {
    int v = G[x][i];
    dfs(v);
    merge(x, v);
  }
  std::multiset<int>::iterator it = st[x] -> lower_bound(d[x]);
  if(it != (st[x] -> end())) st[x] -> erase(it);
  st[x] -> insert(d[x]);
}
 
int main() {
  int n; scanf("%d", &n);
  for(int i = 1; i <= n; i ++) {
    int f; scanf("%d%d", &d[i], &f);
    G[f].push_back(i);
  }
  for(int i = 1; i <= n; i ++) {
    st[i] = new std::multiset<int>();
  }
  dfs(1);
  printf("%d\n", st[1] -> size());
  for(int i = 1; i <= n; i ++) {
    delete st[i];
  }
  return 0;
}

[BZOJ 5093][Lydsy1711月赛]图的价值

统计的时候,考虑每个店对每个图的贡献,可以看出答案是:

\[n2^{\tfrac{(n - 1)(n - 2)}{2}}\,\sum_{i = 0}^{n - 1} \binom{n - 1}{i} i^k\]

求和号前面还好说,主要看求和号后面的咋搞。

后面的那一部分是个经典题吧……但我还是来推导一番吧:

\[
\begin{aligned}
\sum_{i = 0}^{m} \binom{m}{i} i^k&= \sum_{i = 0}^{m} \binom{m}{i}\sum_{j = 0}^k S(k, j) j!\binom{i}{j}\\
&= \sum_{j = 0}^k S(k, j) j! \sum_{i = 0}^m \binom{m}{i}\binom{i}{j} \\
&= \sum_{j = 0}^k S(k, j) m^{\underline{j}} \sum_{i = 0}^{m} \binom{n - j}{i - j} \\
&= \sum_{j = 0}^k S(k, j) m^{\underline{j}} 2^{m - j}
\end{aligned}
\]

然后如果我们能高效的求出某一行的第二类斯特林数,那么问题就迎刃而解了。

考虑斯特林反演(这本质上是一个容斥,用二项式反演也可以推导):

\[S(k, j) = \sum_{i = 0}^j \frac{(-1)^{j - i}}{(j - i)!} \cdot \frac{i^k}{i!}\]

这个东西很显然是一个卷积……而且模数还很友好(UOJ模数),用NTT求出一行的第二类斯特林数即可。

代码:

/**************************************************************
    Problem: 5093
    User: danihao123
    Language: C++
    Result: Accepted
    Time:22632 ms
    Memory:25824 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 998244353LL;
const ll bs = 3LL;
inline ll pow_mod(const ll &a, ll b) {
  ll ans = 1LL, res = a;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
inline ll inv(const ll &x) {
  return pow_mod(x, ha - 2LL);
}
 
inline int flip(const int &bi, const int &x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if((1 << i) & x) {
      ans |= (1 << (bi - i - 1));
    }
  }
  return ans;
}
inline void ntt(ll *A, const int &bi, const int &len, bool flag = false) {
  for(int i = 0; i < len; i ++) {
    int v = flip(bi, i);
    if(v < i) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < len; L <<= 1) {
    ll xi = pow_mod(bs, (ha - 1LL) / (ll(L << 1)));
    if(flag) xi = inv(xi);
    for(int i = 0; i < len; i += (L << 1)) {
      ll w = 1LL;
      for(int j = i; j < i + L; j ++) {
        ll x = A[j], y = A[j + L];
        A[j] = (x + (w * y) % ha) % ha;
        A[j + L] = (x - (w * y) % ha + ha) % ha;
        w = (w * xi) % ha;
      }
    }
  }
}
 
const int maxn = 800005;
ll A[maxn], B[maxn];
ll fac[maxn], down[maxn];
int main() {
  ll n; int k; scanf("%lld%d", &n, &k);
  if(n == 1) {
    puts("0"); return 0;
  }
  int len = 1, bi = 0;
  while(len <= (2 * k)) {
    len <<= 1; bi ++;
  }
  fac[0] = 1LL;
  for(int i = 1; i <= k; i ++) {
    fac[i] = (fac[i - 1] * (ll(i))) % ha;
  }
  for(int i = 0; i <= k; i ++) {
    A[i] = (inv(fac[i]) * pow_mod(ha - 1LL, i)) % ha;
  }
  for(int i = 0; i <= k; i ++) {
    B[i] = inv(fac[i]);
    B[i] = (B[i] * pow_mod(i, k)) % ha;
  }
  ntt(A, bi, len); ntt(B, bi, len);
  for(int i = 0; i < len; i ++) {
    A[i] = (A[i] * B[i]) % ha;
  }
  ntt(A, bi, len, true);
  ll inv_l = inv(len);
  for(int i = 0; i < len; i ++) {
    A[i] = (A[i] * inv_l) % ha;
  }
#ifdef LOCAL
  for(int i = 0; i <= k; i ++) {
    printf("%lld ", A[i]);
  }
  puts("");
#endif
  down[0] = 1LL;
  for(int i = 1; i <= std::min(ll(k), n - 1LL); i ++) {
    down[i] = (down[i - 1] * (n - (ll(i)))) % ha;
  }
  ll ans = 0LL;
  for(int i = 1; i <= std::min(ll(k), n - 1LL); i ++) {
    ll delta = (A[i] * down[i]) % ha;
    delta = (delta * pow_mod(2, n - 1LL - (ll(i)))) % ha;
    ans = (ans + delta) % ha;
  }
  ans = (ans * pow_mod(2LL, (n - 1LL) * (n - 2LL) / 2LL)) % ha;
  ans = (ans * n) % ha;
  printf("%lld\n", ans);
  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 3622]已经没有什么好害怕的了

已经……没有什么好害怕的了

OI给了我足够的欢乐了,OI让我看到了新的世界,让我有了人生中最好的几位朋友,我没有什么可以遗憾的了……就算可能要面临最残酷的结局吧……

有点跑题了呢……还是说这个题吧

显然可以看出来这个题要求你的从大到小的匹配(姑且称之为好的匹配)恰好要有\(\frac{n + k}{2}\)个,剩下的全都不是好的匹配。

首先把糖果(记为\(A\))和药片(记为\(B\))分别排序,对于每一个\(A\)中元素就很容易得到\(B\)中有多少比它小的。考虑设计状态\(f[i][j]\)表示对(排序后的)\(A\)的前\(i\)项,已经确定了的好的匹配有至少\(j\)个。这个DP转移显然。

然后发现\(f[n][j]\times (n - j)!\)就是对于所有元素的完美匹配中,有至少\(j\)个是好的匹配的方案数。这正是广义容斥原理的用武之地。考虑记\(t_j = f[n][j]\times (n - j)!\),那么答案就是(令\(\frac{n + k}{2} = l\)):

\[\sum_{i = l}^n (-1)^{n - i}\binom{i}{l} t_i\]

代码:

/**************************************************************
    Problem: 3622
    User: danihao123
    Language: C++
    Result: Accepted
    Time:1844 ms
    Memory:63680 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000000009LL;
const int maxn = 2005;
ll C[maxn][maxn], F[maxn];
inline void process() {
  int N = 2000;
  C[0][0] = 1LL;
  for(int i = 1; i <= N; i ++) {
    C[i][0] = C[i][i] = 1LL;
    for(int j = 1; j < i; j ++) {
      C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % ha;
    }
  }
  F[0] = 1LL;
  for(int i = 1; i <= N; i ++) {
    F[i] = (F[i - 1] * (ll(i))) % ha;
  }
}
 
int A[maxn], B[maxn];
ll cp[maxn]; int n;
void gen_cp() {
  std::sort(A + 1, A + 1 + n); std::sort(B + 1, B + 1 + n);
  for(int i = 1; i <= n; i ++) {
    cp[i] = n - (B + n + 1 - std::upper_bound(B + 1, B + 1 + n, A[i]));
#ifdef LOCAL
    printf("cp[%d] : %lld\n", i, cp[i]);
#endif
  }
}
 
ll f[maxn][maxn];
ll dp(int dec) {
  process(); gen_cp();
  f[0][0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    for(int j = 0; j <= n; j ++) {
      f[i][j] = f[i - 1][j];
      if(j > 0) {
        f[i][j] += (f[i - 1][j - 1] * std::max(0LL, cp[i] - ll(j - 1))) % ha;
        if(f[i][j] >= ha) f[i][j] -= ha;
      }
#ifdef LOCAL
      printf("f[%d][%d] : %lld\n", i, j, f[i][j]);
#endif
    }
  }
  for(int j = 0; j <= n; j ++) {
    f[n][j] = (f[n][j] * F[n - j]) % ha;
  }
  ll ans = 0;
  for(int j = dec; j <= n; j ++) {
    ll a = ((j - dec) % 2 == 0) ? 1LL : (ha - 1LL);
    ll delta = (C[j][dec] * f[n][j]) % ha;
    ans += (a * delta) % ha;
    if(ans >= ha) ans -= ha;
  }
  return ans;
}
 
int main() {
  int k; scanf("%d%d", &n, &k);
  if((n - k) & 1) {
    puts("0"); return 0;
  }
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &A[i]);
  }
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &B[i]);
  }
  printf("%lld\n", dp((n + k) / 2));
  return 0;
}

[BZOJ 5091][Lydsy0711月赛]摘苹果

秒,秒啊.jpg

首先根据期望线性性,每个点的期望可以分开算。不妨设\(f_{i, j}\)表示\(j\)在某个第\(i\)步被走到的概率。那么每个点\(j\)的期望访问次数就是\(\sum_{i = 0}^k f_{i, j}\)。

然后考虑去求那个\(f_{i, j}\)。显然\(f_{0, j} = \frac{d_j}{2m}\)。但是考虑\(f_{1, j}\)的转移方程:

\[f_{1, j} = \sum \frac{f_{0, u}}{d_u}\]

然后展开之后发现每个\(\frac{f_{0, u}}{d_u}\)都是\(\frac{1}{2m}\),于是乎\(f_{1, j} = \frac{d_j}{2m}\)。

如此一来,对于任意\(i\),都有\(d_{i, j} = \frac{d_j}{2m}\),然后就很好做了……

代码:

/**************************************************************
    Problem: 5091
    User: danihao123
    Language: C++
    Result: Accepted
    Time:532 ms
    Memory:2384 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000000007LL;
ll pow_mod(ll a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b >> 1);
  ans = (ans * ans) % ha;
  if(1LL & b) ans = (ans * a) % ha;
  return ans;
}
ll inv(ll v) {
  return pow_mod(v, ha - 2LL);
}
 
const int maxn = 100005;
ll a[maxn], d[maxn];
int main() {
  int n, m, k; scanf("%d%d%d", &n, &m, &k);
  for(int i = 1; i <= n; i ++) {
    scanf("%lld", &a[i]);
  }
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v);
    d[u] ++; d[v] ++;
  }
  ll ans = 0;
  ll inv_2m = inv(2 * m);
  for(int i = 1; i <= n; i ++) {
    ans += (((d[i] * inv_2m) % ha) * a[i]) % ha;
    if(ans >= ha) ans -= ha;
  }
  ans = (ans * (ll(k))) % ha;
  printf("%lld\n", ans);
  return 0;
}

[BZOJ 3156]防御准备

又做了一个简单的斜率优化题TAT

首先,设\(f_i\)表示最后一个放塔的点事\(i\)时的最优解,那么将原方程化简得:

\[f_i = a_i + \frac{i^2 - i}{2} + max(-ij + \frac{j^2 + j}{2} + f_j)\]

然后求直线形式,得到:

\[ij + d_i = f_j + \frac{j^2 + j}{2}\]

用类似于玩具装箱(上一篇题解)的方式搞一搞即可。

代码:

/**************************************************************
	Problem: 3156
	User: danihao123
	Language: C++
	Result: Accepted
	Time:2496 ms
	Memory:16480 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <deque>
#include <cmath>
#include <climits>
typedef long long ll;
typedef ll T;
struct Point {
  T x, y;
  Point(T qx = 0LL, T qy = 0LL) {
    x = qx; y = qy;
  }
};
typedef Point Vector;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
Vector operator *(const Vector &a, T lam) {
  return Vector(a.x * lam, a.y * lam);
}
Vector operator *(T lam, const Vector &a) {
  return Vector(a.x * lam, a.y * lam);
}
inline T dot(const Vector &a, const Vector &b) {
  return (a.x * b.x + a.y * b.y);
}
inline T times(const Vector &a, const Vector &b) {
  return (a.x * b.y - a.y * b.x);
}

const int maxn = 1000005;
T f[maxn], a[maxn];
int n;
void dp() {
  f[1] = a[1];
  std::deque<Point> Q;
  Q.push_back(Point(1LL, f[1] + 1LL));
  for(T i = 2; i <= n; i ++) {
    T k = i;
    Vector st(1LL, k);
    while(Q.size() > 1 && times(Q[1] - Q[0], st) > 0LL) {
      Q.pop_front();
    }
    f[i] = a[i] + ((i * i - i) / 2LL);
    f[i] += Q.front().y - Q.front().x * i;
    Point ins(i, f[i] + ((i * i + i) / 2LL));
    while(Q.size() > 1 && times(ins - Q.back(), Q.back() - Q[Q.size() - 2]) > 0LL) {
      Q.pop_back();
    }
    Q.push_back(ins);
  }
}

int main() {
  scanf("%d", &n);
  for(int i = n; i >= 1; i --) {
    scanf("%lld", &a[i]);
  }
  dp();
  ll ans = f[n];
#ifdef LOCAL
  printf("f[n] : %lld\n", ans);
#endif
  for(T i = 1; i < n; i ++) {
#ifdef LOCAL
    printf("f[%d] : %lld\n", i, f[i]);
#endif
    ans = std::min(ans, f[i] + (n - i + 1LL) * (n - i) / 2LL);
  }
  printf("%lld\n", ans);
  return 0;
}

[BZOJ 1010][HNOI2008]玩具装箱toy

很久之前是学过并写过斜率优化的……但是很快就忘了。现在感觉自己理解了,感觉是真的懂了……抽空写篇文章解释一下吧……

先单独说这一个题。将DP方程完全展开,并且设\(P_i = S_i + i\),\(c = L + 1\),可得:

\[f_i = c^2 + P_i^2 - 2P_i c + max(P_j^2 + 2P_j c + f_j - 2P_i P_j)\]

然后\(c^2 + P_i^2 - 2P_i c\)这部分是常数项不需要管了,我们就想想max里面那些(姑且设之为\(d_i\))咋整好了。

设\(d_i = P_j^2 + 2P_j c + f_j - 2P_i P_j\),稍作移项,得:

\[2P_i P_j + d_i = P_j^2 + 2P_j c + f_j\]

于是乎,\(d_i\)可以看做斜率为\(2P_i\)的直线过点\((P_j, P_j^2 + 2P_j c + f_j)\)得到的截距。而那些点我们之前都知道了,问题就变成了已知斜率,求过某点集中的点的最大截距。

想象一个固定斜率的直线从下往上扫,那么碰到的第一个点就是最优解。首先这个点一定在下凸壳上,其次下凸壳上这点两侧的线段的斜率肯定一个比\(2P_i\)大另一个比它小。并且最好的一点是这个斜率还是单调的,那么分界点一定是单调递增的。

代码:

/**************************************************************
    Problem: 1010
    User: danihao123
    Language: C++
    Result: Accepted
    Time:132 ms
    Memory:2416 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <deque>
#include <cmath>
typedef long long ll;
typedef ll T;
struct Point {
  T x, y;
  Point(T qx = 0LL, T qy = 0LL) {
    x = qx; y = qy;
  }
};
typedef Point Vector;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
Vector operator *(const Vector &a, T lam) {
  return Vector(a.x * lam, a.y * lam);
}
Vector operator *(T lam, const Vector &a) {
  return Vector(a.x * lam, a.y * lam);
}
inline T dot(const Vector &a, const Vector &b) {
  return (a.x * b.x + a.y * b.y);
}
inline T times(const Vector &a, const Vector &b) {
  return (a.x * b.y - a.y * b.x);
}
 
const int maxn = 50005;
T C[maxn], S[maxn], P[maxn];
T f[maxn];
int n; ll c;
void process() {
  for(int i = 1; i <= n; i ++) {
    S[i] = S[i - 1] + C[i];
    P[i] = S[i] + (ll(i));
  }
}
void dp() {
  std::deque<Point> Q;
  Q.push_back(Point(0LL, 0LL));
  for(int i = 1; i <= n; i ++) {
    ll k = 2 * P[i];
    Vector st(1, k);
    while(Q.size() > 1 && times(Q[1] - Q[0], st) > 0LL) {
      Q.pop_front();
    }
    f[i] = c * c + P[i] * P[i] - 2LL * P[i] * c;
    f[i] += Q.front().y - k * Q.front().x;
#ifdef LOCAL
    printf("f[%d] : %lld\n", i, f[i]);
#endif
    Vector ins(P[i], f[i] + P[i] * P[i] + 2LL * P[i] * c);
    while(Q.size() > 1 && times(ins - Q.back(), Q.back() - Q[Q.size() - 2]) > 0LL) {
#ifdef LOCAL
      printf("Deleting (%lld, %lld)...\n", Q.back().x, Q.back().y);
#endif
      Q.pop_back();
    }
    Q.push_back(ins);
#ifdef LOCAL
    printf("Inserting (%lld, %lld)...\n", ins.x, ins.y);
#endif
  }
}
 
int main() {
  scanf("%d%lld", &n, &c); c ++;
  for(int i = 1; i <= n; i ++) {
    scanf("%lld", &C[i]);
  }
  process(); dp();
  printf("%lld\n", f[n]);
  return 0;
}

[BZOJ 2561]最小生成树

窝只会做水体了qwq

因为这条边既存在在最大生成树里,又存在在最小生成树里。那就说明\(u\)到\(v\)的路径上,在最大生成树上这条边是最小的,在最小生成树上这条边是最大的。

所以说我们不能用大于\(L\)的边来联通\(u\)和\(v\),也不能用小于\(L\)的边,于是乎……跑两次最小割即可。

代码:

/**************************************************************
	Problem: 2561
	User: danihao123
	Language: C++
	Result: Accepted
	Time:1528 ms
	Memory:15952 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <utility>
#include <algorithm>
#include <queue>
const int maxn = 20005;
const int maxm = 200005;
const int maxe = maxm << 2;
int n, m;
int first[maxn];
int next[maxe], to[maxe], flow[maxe], cap[maxe];
int gcnt = 0;
inline void init_graph() {
  gcnt = 0;
  std::fill(first + 1, first + 1 + n, 0);
}
inline void add_edge(int u, int v, int f) {
  gcnt ++;
  next[gcnt] = first[u];
  first[u] = gcnt;
  to[gcnt] = v;
  flow[gcnt] = 0; cap[gcnt] = f;
}
inline int rev(int i) {
  return (((i - 1) ^ 1) + 1);
}
inline void ins_edge(int u, int v, int f) {
  add_edge(u, v, f);
  add_edge(v, u, 0);
}

int line[maxm][3];
int d[maxn];
int s, t;
inline bool bfs() {
  static bool vis[maxn];
  std::fill(vis + 1, vis + 1 + n, false);
  std::fill(d + 1, d + 1 + n, 0);
  std::queue<int> Q;
  Q.push(s); vis[s] = true; d[s] = 1;
  while(!Q.empty()) {
    int u = Q.front(); Q.pop();
    for(int i = first[u]; i; i = next[i]) {
      int v = to[i];
      if(!vis[v] && cap[i] > flow[i]) {
        d[v] = d[u] + 1; vis[v] = true;
        Q.push(v);
      }
    }
  }
  return vis[t];
}
int cur[maxn];
int dfs(int x, int a) {
  if(a == 0 || x == t) return a;
  int ans = 0;
  for(int &i = cur[x]; i; i = next[i]) {
    int v = to[i];
    int f;
    if(d[v] == d[x] + 1 && (f = dfs(v, std::min(a, cap[i] - flow[i]))) > 0) {
      ans += f; a -= f;
      flow[i] += f; flow[rev(i)] -= f;
      if(a == 0) break;
    }
  }
  if(a > 0) d[x] = -1;
  return ans;
}
int dinic() {
  int ans = 0;
  while(bfs()) {
    for(int i = 1; i <= n; i ++) cur[i] = first[i];
    ans += dfs(s, 0x7fffffff);
  }
  return ans;
}

int main() {
  scanf("%d%d", &n, &m);
  for(int i = 1; i <= m; i ++) {
    int *l = line[i];
    scanf("%d%d%d", &l[0], &l[1], &l[2]);
  }
  int L; scanf("%d%d%d", &s, &t, &L);
  for(int i = 1; i <= m; i ++) {
    int *l = line[i];
    if(l[2] < L) {
      ins_edge(l[0], l[1], 1);
      ins_edge(l[1], l[0], 1);
    }
  }
  int ans = dinic();
  init_graph();
  for(int i = 1; i <= m; i ++) {
    int *l = line[i];
    if(l[2] > L) {
      ins_edge(l[0], l[1], 1);
      ins_edge(l[1], l[0], 1);
    }
  }
  ans += dinic();
  printf("%d\n", ans);
  return 0;
}