2020 Apr Records | fa

菜鸡 fa_555 会把一些动态简要地记在这里。主要还是留给自己以后看的。

所有代码不保证包括无关紧要的部分

[SCOI2007]k短路

solution

~~这是一道卡 A* 的 A* 模板题~~

简单说一下 A* 算法。对于一个搜索状态，设计函数 $f, g, h, h^\ast$。其中 $g$ 为搜索到当前状态的代价；$h$ 为估计剩余的代价；$h^\ast$ 为实际剩余的代价，$h < h^\ast$，且差距越小效率越高；$f = g + h$，对 $f$ 维护一个优先队列即可。

code

int N, K;

struct Edge { int nxt, to, w; };

namespace G1 {
  std::bitset<51> vis;
  int head[51], dis[51];
  Edge e[2451];

  inline void addEdge(int u, int v, int w);

  void spfa(int s);
} // G1

namespace G2 {
  int *h = G1::dis, head[51];
  Edge e[2451];

  struct Node {
    std::bitset<51> vis;
    int s, u, len, f;
    std::vector<int> z;

    bool operator<(const Node &rhs) const {
      if (f != rhs.f) return f > rhs.f;
      for (int i = 0, mx = std::min(z.size(), rhs.z.size()); i < mx; ++i)
        if (z[i] != rhs.z[i]) return z[i] > rhs.z[i];
      return z.size() > rhs.z.size();
    }
  };

  inline void addEdge(int u, int v, int w);

  void Astar(int s, int t) {
    std::priority_queue<Node> q;
    q.push({std::bitset<51>(1ll << s), s, s, 0, h[s], std::vector<int>()});
    int cnt = 0;
    for (Node now, nxt; q.size(); ) {
      now = q.top(), q.pop();
      if (now.u == t && ++cnt == K) {
        cout << now.s;
        for (int u : now.z)
          cout << -u;
        cout.flush();
        exit(0);
      }
      for (int i = head[now.u], v; i; i = e[i].nxt) {
        if (now.vis[v = e[i].to]) continue;
        nxt = {now.vis, now.s, v, now.len + e[i].w, h[v], now.z};
        nxt.vis.set(v), nxt.z.push_back(v), nxt.f += nxt.len;
        q.push(nxt);
      }
    }
  }
} // G2

int main() {
  int M, S, T;

  cin >> N >> M >> K >> S >> T;
  for (int u, v, w; M; --M) {
    cin >> u >> v >> w;
    G1::addEdge(v, u, w), G2::addEdge(u, v, w);
  }
  G1::spfa(T), G2::Astar(S, T);
  puts("No");
  return 0;
}

[NOI2005]智慧珠游戏

DLX 是之前学的，这里记一记怎么建出一个精确覆盖模型。

对整个拼盘按照一定的方法编号，每个格子对应一个号码（下面的代码采用先从左到右，再从上到下依次编号）。
把所有零件的形状存下来。这里采取的是确定一个关键点，其余点用相对坐标表示的方法。
设一个有若干行、每行 67 列的 01 矩阵，每行表示某个零件的一种摆法。其中左 12 列表示 12 个零件种类（第 $i$ 列为 $1$ 表示该行代表第 $i$ 种零件的某种摆法）；右 55 列表示拼盘上的 55 个空位，第 $i$ 列为 $1$ 表示这种摆法占据了第 $i$ 个格子。
对于 12 种零件分别枚举摆法。
- 首先枚举上文提到的“关键点”的坐标，一共 55 种情况。
- 然后要枚举它的旋转、翻转变换。然而旋转操作并不方便编程直接实现，要考虑怎么实现。画一画图发现一共有 8 种不同摆法（如图）
- 发现可以通过三种简单的变换得到：交换横纵坐标、横坐标取相反数、纵坐标取相反数。
把枚举出来的所有状态全扔进一个矩阵，DLX 跑个精确覆盖，完事了。统计答案有很多种方法，可以各路神仙八仙过海一下（
对于全空拼盘，理论上矩阵最大行数 $12 \times 55 \times 8 = 5280$，实际上由于零件形状的制约只有 2730 种，矩阵中的 $1$ 最多 15083 个。

code

对于上面的三种变换，这里通过状态 $S \in [0, 8)$ 来枚举，$S$ 的最低三位从高到低依次表示是否交换横纵坐标、横坐标是否相反数、纵坐标是否取相反数。

部分变量解释：

siz[i] 表示第 $i$ 种零件的珠子个数
sh[i][][] 表示第 $i$ 种零件的形状（见上文）


constexpr int siz[] = {
  0, 3, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5
}, sh[13][5][2] = {{},
  {{0, 0}, {1, 0}, {0, 1}},
  {{0, 0}, {0, 1}, {0, 2}, {0, 3}},
  {{0, 0}, {1, 0}, {0, 1}, {0, 2}},
  {{0, 0}, {1, 0}, {0, 1}, {1, 1}},
  {{0, 0}, {1, 0}, {2, 0}, {0, 1}, {0, 2}},
  {{0, 0}, {-1, 0}, {0, 1}, {1, 0}, {2, 0}},
  {{0, 0}, {0, 1}, {1, 0}, {2, 0}, {2, 1}},
  {{0, 0}, {0, 1}, {1, 0}, {2, 0}, {1, 1}},
  {{0, 0}, {1, 0}, {2, 0}, {2, 1}, {3, 1}},
  {{0, 0}, {-1, 0}, {0, -1}, {1, 0}, {0, 1}},
  {{0, 0}, {0, 1}, {1, 1}, {1, 2}, {2, 2}},
  {{0, 0}, {0, 1}, {1, 0}, {2, 0}, {3, 0}}
};

using Data = std::vector<std::pair<int, int> >;

std::bitset<13> used;
char z[13][13];
Data v[13];

namespace DLX {
  using iarr = int[15085];

  std::bitset<69> mp[2733];
  int tot, ans[15], cnt[69], h[2733];
  iarr l, r, u, d, x, y;

  void init(int M);
  inline void link(int R, int C);
  void remove(int C);
  void resume(int C);

  void dance(int dep) {
    if (!r[0]) {
      /* output */
      exit(0);
    }
    /* ... */
  }
} // DLX

inline int pos(int i, int j) {
  static int to[] = {0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55};
  return to[i - 1] + j;
}

void insert(int c, const Data &v) {
  // v 中的元素为颜色 c 某种状态占据的坐标
  static int M = 0;
  DLX::link(++M, c);
  for (const auto &x : v)
    DLX::link(M, 12 + pos(x.first, x.second));
}

inline bool valid(int x, int y) {
  return 1 <= x && x <= 10 && 1 <= y && y <= x && z[x][y] == '.';
}

int main() {
  for (int i = 1; i <= 10; ++i)
    for (int j = 1; j <= i; ++j) {
      std::cin >> z[i][j];
      if (z[i][j] != '.') {
        used.set(z[i][j] & 31);
        v[z[i][j] & 31].push_back({i, j});
      }
    }

  DLX::init(67);
  for (int i = 1; i <= 12; ++i) {
    if (used[i]) {
      insert(i, v[i]);
      continue;
    }
    v[i].resize(siz[i]);
    for (int pi = 1; pi <= 10; ++pi)
      for (int pj = 1; pj <= pi; ++pj)
        for (int S = 0; S < 8; ++S) {
          bool g = 0;
          for (int j = 0, nx, ny; j < siz[i]; ++j) {
            nx = sh[i][j][0], ny = sh[i][j][1];
            if (S & 4) std::swap(nx, ny);
            if (S & 2) nx = -nx;
            if (S & 1) ny = -ny;
            nx += pi, ny += pj;
            if (!valid(nx, ny)) { g = 1; break; }
            v[i][j] = {nx, ny};
          }
          if (!g) insert(i, v[i]);
        }
  }

  DLX::dance(1);
  puts("No solution"); // 有解的情况在 DLX 中直接退出
  return 0;
}

数学

在 CSG 那里学习了一个。

本来想搞多项式，结果被数论水平限制了，去从零重学数论，没想到一学就是好长好长好长时间 /kk

学习杜教筛

双亲数发现了自己以前题解的错误并重写

类欧几里得算法花两天边颓边抄了一遍式子又花一天写出来代码，~~因为懒~~研究了各种不定参的写法。

BSGS

学了 BSGS 和 exBSGS。~~好难啊直接背吧~~

内附一个 Hash 表板子，机制类似于 std::unordered_map 但是跑得快（

using ll = long long;

constexpr int Hmod = 1000003;

struct Hash {
  std::bitset<Hmod> vis;
  int tv[Hmod], tp[Hmod];

  void clear() {
    vis.reset();
  }

  void insert(int p, int v) {
    int k = p % Hmod;
    for (; vis[k]; ++k %= Hmod)
      if (tp[k] == p) {
        tv[k] = v;
        return;
      }
    vis.set(k), tv[k] = v, tp[k] = p;
  }

  int query(int p) {
    for (int k = p % Hmod; vis[k]; ++k %= Hmod)
      if (tp[k] == p) return tv[k];
    return -1;
  }
} z;

int BSGS(int a, int b, int p, int sd) {
  z.clear();
  int v = 1, t = ceil(sqrt(p));

  for (int i = 0; i < t; ++i, v = 1ll * v * a % p)
    z.insert(1ll * v * b % p, i);
  for (int i = 0, j, v_ = v, v = sd; i <= t; ++i, v = 1ll * v * v_ % p)
    if ((j = z.query(v)) != -1 && i * t >= j) return i * t - j;
  return -1;
}

int exBSGS(int a, int b, int p) {
  a %= p, b %= p;
  if (b == 1 || p == 1) return 0;
  int cnt = 0, sd = 1;
  for (int d; (d = std::__gcd(a, p)) != 1; ) {
    if (b % d) return -1;
    ++cnt, b /= d, p /= d;
    sd = 1ll * sd * a / d % p;
    if (sd == b) return cnt;
  }
  int ans = BSGS(a, b, p, sd);
  if (ans != -1) ans += cnt;
  return ans;
}

原根

关于原根的笔记在这里

[SDOI2010]古代猪文

problem

给定正整数 $n, g \ (n, g \le 10^9)$，求 $g^{ \sum_{d \mid n}{ \binom nd}} \bmod 999911659$。

solution

本能考虑欧拉定理，则答案等于 $g^{ \sum_{d \mid n}{ \binom nd} \bmod 999911658} \bmod 999911659$。

设 $x = \sum_{d \mid n}{ \binom nd} \bmod 999911658$，重点在于求出 $x$。怎么求呢？

Extended Lucas

这个指数可以直接无脑上 Extended Lucas 嘛对不对（

这个复杂度上天的做法在 LOJ 是最劣解第一，在谷直接 TLE 成了 70 分（

Lucas + CRT

考虑对模数进行质因数分解，$999911658 = 2 \times 3 \times 4679 \times 35617$。发现四个因子都是质数而不是质数的幂，因此可以使用 Lucas 计算出

$\begin{cases} x \equiv r_1 \pmod 2 \\ x \equiv r_2 \pmod 3 \\ x \equiv r_3 \pmod{4679} \\ x \equiv r_4 \pmod{35617} \end{cases}$

中的 $r_1, r_2, r_3, r_4$，再使用 CRT 合并得到 $x$。

code

这个做法我的实现在 LOJ 得到了最优解第二。~~如果我早交 11 个小时就能当第一~~

using ll = long long;

constexpr int m[5] = {999911659, 2, 3, 4679, 35617};

int r[5], fac[5][35619], invF[5][35619];

int Pow(ll b, ll p, ll m);

void prework(); // 在 fac[k][i] 和 invF[k][i] 中求出模 m[k] 意义下 i 的阶乘及其逆元
int C(int N, int M, int k);
int Lucas(int N, int M, int k); // 求 C(N, M) % m[k]

inline void calc(int N, int i) {
  for (int k = 1; k <= 4; ++k)
    r[k] = (r[k] + Lucas(N, i, k)) % m[k];
}

int CRT(int N, int *r, const int *m) {
  constexpr int p = *m - 1;
  int ans = 0;
  for (int k = 1; k <= N; ++k)
    ans = (ans + 1ll * r[k] * (p / m[k]) % p * Pow(p / m[k], m[k] - 2, m[k]) % p) % p;
  return ans;
}

int solve(int N, int g) {
  if (g % *m == 0) {
    std::cout << 0;
    return 0;
  }

  prework();

  for (int i = 1; i * i <= N; ++i)
    if (N % i == 0) {
      calc(N, i);
      if (i * i != N) calc(N, N / i);
    }
  return Pow(g, CRT(4, r, m), *m);
}

[HAOI2011]Problem b

problem

$T \ (1 \le T \le 5 \times 10^4)$ 组数据，给定 $a, b, c, d, k \ (1 \le a, b, c, d, k \le 5 \times 10^4)$，求 $\sum{i = a}^b{ \sum{j = c}^d{[\gcd(i, j) = k]}}$。

solution

我以前写过这道题的题解对吧（

一模一样，套个傻逼容斥就完了，总复杂度 $O(n + T \sqrt n)$，带一个不小的常数。

LCMSUM

problem

$T \ (1 \le T \le 3 \times 10^5)$ 组数据，给定 $n \ (1 \le n \le 10^6)$，求 $\sum_{i = 0}^n{ \mathrm{lcm}(i, n)}$。

solution

设 $d’ = n / d$，则

$\begin{align} \sum_{i = 0}^n{ \mathrm{lcm}(i, n)} &= \sum_{i = 1}^n \frac{in}{ \gcd(i, n)} \\ &= \frac 12 \left( \sum_{i = 1}^{n - 1}{ \frac{in}{ \gcd(i, n)} + \sum_{i = n - 1}^1{ \frac{in}{ \gcd(i, n)}}} \right) + n \\ &= \frac 12 \left( \sum_{i = 1}^{n - 1}{ \frac{in}{ \gcd(i, n)} + \sum_{i = n - 1}^1{ \frac{in}{ \gcd(n - i, n)}}} \right) + n \\ &= \frac 12 \sum_{i = 1}^{n - 1}{ \frac{n^2}{ \gcd(i, n)}} + n \\ &= \frac 12 \sum_{i = 1}^n{ \frac{n^2}{ \gcd(i, n)}} + \frac n2 \\ &= \frac 12 \sum_{d \mid n}{\frac{n^2 \varphi(d')}d} + \frac n2 \\ &= \frac n2 \left( \sum_{d \mid n}{d' \varphi(d') + 1} \right) \end{align}$

设 $f(n) = \sum_{d \mid n} d \ \varphi(d)$。

$f$ 是积性函数，可以 $O(n) - O(1)$，总复杂度 $O(n + T)$。$O(n)$ 筛这个函数的操作好像很骚的样子，不过也有 $O(n \log n)$ 的用 $\varphi$ 间接求的算法。

code

std::bitset<1000003> v;
int p[78501];
long long f[1000003];

static int sieve = ([]() {
  constexpr int N = 1000000;
  int M = 0;
  v.set(1), f[1] = 1;
  for (int i = 2; i <= N; ++i) {
    if (!v[i]) p[++M] = i, f[i] = 1ll * i * (i - 1) + 1;
    for (int j = 1; j <= M && i * p[j] <= N; ++j) {
      v.set(i * p[j]);
      if (i % p[j] == 0) {
        if (i / p[j] % p[j] == 0)
          f[i * p[j]] = f[i] + (f[i] - f[i / p[j]]) * p[j] * p[j];
        else
          f[i * p[j]] = f[i] + f[i / p[j]] * (p[j] - 1) * p[j] * p[j] * p[j];
        break;
      }
      f[i * p[j]] = f[i] * f[p[j]];
    }
  }
  return 0;
})();

int solve(int N) {
  return (f[N] * N + N) >> 1;
}

Crash的数字表格

problem

给定 $n, m \ (1 \le n, m, \le 10^7)$，求 $\sum{i = 1}^n{ \sum{j = 1}^m{ \mathrm{lcm}(i, j)}}$ 在模意义下的值。

solution

不妨设 $n \le m$。

$\begin{align} \sum_{i = 1}^n{ \sum_{j = 1}^m{ \mathrm{lcm}(i, j)}} &= \sum_{i = 1}^n{ \sum_{j = 1}^m{ \frac{ij}{ \gcd(i, j)}}} \\ &= \sum_{i = 1}^n{ \sum_{j = 1}^m{ \sum_{d \mid \gcd(i, j), \left[ \frac id \perp \frac jd \right]}{ \frac{ij}d}}} \\ &= \sum_{d = 1}^n{d \sum_{i = 1}^{ \lfloor n / d \rfloor}{ \sum_{j = 1}^{ \lfloor m / d \rfloor}{[i \perp j] \ ij}}} \end{align}$

令 $S(n, m) = \sum{i = 1}^n{ \sum{j = 1}^m{[i \perp j] \ ij}}$，则

$\begin{align} S(n, m) &= \sum_{i = 1}^n{ \sum_{j = 1}^m{ \epsilon( \gcd(i, j)) ij}} \\ &= \sum_{d = 1}^n{ \sum_{d \mid i}^n{ \sum_{d \mid j}^m{ \mu(d)ij}}} \\ &= \sum_{d = 1}^n{ \mu(d)d^2 \sum_{i = 1}^{n / d}{ \sum_{j = 1}^{m / d}{ij}}} \end{align}$

令 $f(n, m) = \sum{i = 1}^n{ \sum{j = 1}^m ij} = \frac{n(n + 1)}2 \times \frac{m(m + 1)}2$，则

$S(n, m) = \sum_{d = 1}^n{ \mu(d) · d^2 · f \left( \left\lfloor \frac nd \right\rfloor, \left\lfloor \frac md \right\rfloor \right)}$

可以数论分块 $O \left( \left\lfloor \frac{n}{ \lfloor n / d \rfloor} \right\rfloor \right)$ 求 $S(n, m)$。则原答案为

$\sum_{d = 1}^n{d · S \left( \left\lfloor \frac nd \right\rfloor, \left\lfloor \frac md \right\rfloor \right)}$

code

注意代码中有两处 assert，改变其数量，在评测鸭实测速度 2 > 0 > 1，不知为何。~~可能是玄学~~

constexpr int mod = 20101009;

std::bitset<10000003> v;
int p[664581], mu[10000003], g[10000003];

void sieve(int N) {
  /* 筛出 mu[1..N] */
  for (int i = 1; i <= N; ++i)
    g[i] = (g[i - 1] + 1ll * i * i % mod * (mu[i] + mod) % mod) % mod;
}

int f(int N, int M) {
  return (1ll * N * (N + 1) / 2 % mod) * (1ll * M * (M + 1) / 2 % mod) % mod;
}

int S(int N, int M) {
  assert(N <= M);
  int ans = 0;
  for (int l = 1, r, n, m; l <= N; l = r + 1) {
    n = N / l, m = M / l;
    r = std::min(N / n, M / m);
    ans = (ans + 1ll * (g[r] - g[l - 1] + mod) * f(n, m) % mod) % mod;
  }
  return ans;
}

int solve(int N, int M) {
  assert(N <= M);
  int ans = 0;
  for (int l = 1, r, n, m; l <= N; l = r + 1) {
    n = N / l, m = M / l;
    r = std::min(N / n, M / m);
    ans = (ans + 1ll * (r - l + 1) * (l + r) / 2 % mod * S(n, m) % mod) % mod;
  }
  return ans;
}

int main() {
  int N, M;

  std::cin >> N >> M;
  if (N > M) std::swap(N, M);
  sieve(N);
  std::cout << solve(N, M);
  return 0;
}

[SDOI2015约数个数和]

problem

$T \ (1 \le T \le 50000)$ 组数据，给定 $n, m \ (1 \le n, m \le 50000)$，求 $\sum{i = 1}^n{ \sum{j = 1}^m{ \mathrm d(ij)}}$。其中 $\mathrm d$ 为约数个数函数。

solution

首先有个很少人会证的结论

$\mathrm d(ij) = \sum_{x \mid i}{ \sum_{y \mid j}{[x \perp y]}}$

然后可以通过一系列简单的操作得到答案为（设 $n \le m$）

$\sum_{d = 1}^n{ \mu(d) · S_\mathrm d \left( \left\lfloor \frac nd \right\rfloor \right) · S_\mathrm d \left( \left\lfloor \frac md \right\rfloor \right)}$

其中 $S_\mathrm d$ 为 $\mathrm d$ 的前缀和函数。外层显然数论分块，总复杂度是 $O(n + T \sqrt n)$。

[ZJOI2014]力

problem

给出 $n \ (1 \le n \le 10^5)$ 个数 $q_i \ (1 \le q_i < 10^9)$，定义

$F_i = \sum_{j = 1}^{i - 1}{\frac{q_i q_j}{(i - j)^2}} - \sum_{j = i + 1}^n \frac{q_i q_j}{(i - j)^2} \\ E_i = \frac{F_i}{q_i}$

求 $E_{1 \dots n}$。

solution

$E_i = \frac{F_i}{q_i} = \sum_{j = 1}^{i - 1}{\frac{q_j}{(i - j)^2}} - \sum_{j = i + 1}^n \frac{q_j}{(i - j)^2}$

令

$f_i = \begin{cases} 0 & , i = 0 \\ q_i & , i > 0 \end{cases}, \quad g_i = \begin{cases} 0 & , i = 0 \\ \frac 1{i^2} & , i > 0 \end{cases}$

则

$\begin{align} E_i &= \sum_{j = 0}^{i}{f_j g_{i - j}} - \sum_{j = i}^n{f_j g_{j - i}} \\ &= \sum_{j = 0}^{i}{f_j g_{i - j}} - \sum_{j = 0}^{n - i}{f_{i + j} g_{j}} \end{align}$

令 $f’i = f{n - i}, t = n - i$，则

$E_i = \sum_{j = 0}^{i}{f_j g_{i - j}} - \sum_{j = 0}^t{f'_{t - j} g_{j}}$

式子推完了，开始整多项式。令（注意此处 $F$ 与题目给出的 $F$ 不同）

$\begin{align} A(x) &= \sum_{i = 0}^n{f_i x^i} \\ B(x) &= \sum_{i = 0}^n{f'_i x^i} \\ C(x) &= \sum_{i = 0}^n{g_i x^i} \\ F(x) &= A(x) \ast C(x) \\ G(x) &= B(x) \ast C(x) \\ E_i &= F_i - G_{n - i} \end{align}$

code

double ans[50003];

double* solve(int N, double *p) {
  int lim, l;

  for (lim = 1, l = -1; lim <= N << 1; lim <<= 1, ++l);
  for (int i = 0; i < lim; ++i)
    to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

  for (int i = 1; i <= N; ++i) {
    A[i].x = B[N - i].x = p[i];
    C[i].x = 1. / (1. * i * i);
  }

  FFT(lim, A, 1), FFT(lim, B, 1), FFT(lim, C, 1);
  for (int i = 0; i < lim; ++i)
    A[i] = A[i] * C[i], B[i] = B[i] * C[i];
  FFT(lim, A, -1), FFT(lim, B, -1);

  for (int i = 1; i <= N; ++i)
    ans[i] = A[i].x - B[N - i].x;
}

[AHOI2017/HNOI2017]礼物

solution

枚举对第一个手环整体加 $c \ (-m < c < m)$，那么两个手环之间的差异值为

$\begin{align} & \sum_{i = 1}^n{(a_i + c - b_i)^2} \\ =& \sum_{i = 1}^n{(a_i^2 + b_i^2 + c^2 + 2 a_i c - 2 b_i c - 2 a_i b_i)} \\ =& \sum_{i = 1}^n{(a_i^2 + b_i^2)} + 2 c \sum_{i = 1}^n{(a_i - b_i)} - 2 \sum_{i = 1}^n{a_i b_i} + nc^2 \end{align}$

发现除了最后一项以外都是确定的。那么只要令 $\sum_{i = 1}^n{a_i b_i}$ 最大就可以得到最小值了。

首先把 $a$ 序列反向一下，得到 $\sum{i = 1}^n{a{n - i + 1}^n b_i}$，然后在反向后的 $a$ 倍长破环成链，与 $b$ 卷起来得到的新多项式的第 $n + 1\sim 2n$ 项中找最大值即可。

code

int solve(int N, int *a, int M, int *b) {
  int lim, l, ans = 0x3f3f3f3f, ans1 = 0, ans2 = 0, ans3 = 0;

  for (lim = 1, l = -1; lim <= N * 3; lim <<= 1, ++l);
  for (int i = 0; i < lim; ++i)
    to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

  for (int i = 0, t; i < N; ++i) {
    t = a[i + 1];
    A[N - i].x = A[(N << 1) - i].x = t;
    ans1 += t * t, ans2 += t << 1;
  }
  for (int i = 1, t; i <= N; ++i) {
    t = b[i];
    B[i].x = t;
    ans1 += t * t, ans2 -= t << 1;
  }

  FFT(lim, A, 1), FFT(lim, B, 1);
  for (int i = 0; i < lim; ++i)
    A[i] = A[i] * B[i];
  FFT(lim, A, -1);
  for (int i = N + 1; i <= N << 1; ++i)
    ans3 = std::max(ans3, int(A[i].x + 0.5));
  for (int i = -M; i <= M; ++i)
    ans = std::min(ans, i * ans2 + N * i * i);
  ans += ans1 - (ans3 << 1);
  printf("%d\n", ans);
  return 0;
}

SP4882 DAGCNT2 - Counting in a DAG

problem

$T \ (1 \le T \le 2)$ 组数据，给定 $n \ (1 \le n \le 2 \times 10^4)$ 个点 $m \ (1 \le m \le 5 \times 10^5)$ 条边的 DAG，点带权，求每个点出发所能到达的点的点权和。

solution

啊，这个问题的复杂度已经被证明下界是 $O(n^2)$ 的了（

我非常好奇怎么过这个数据范围，于是在网上找到一篇神仙状压的代码并抄了一遍（

code

using iarrN = int[20003];

std::bitset<20003> vis;
int top, f[20003][627];
iarrN mp, rmp, s, w, outd, ans;
std::vector<int> G[20003], R[20003], D[20003];

void init(int N) {
  vis.reset();
  memset(outd + 1, 0, sizeof(int) * N);
  memset(ans + 1, 0, sizeof(int) * N);
  for (int i = 1; i <= N; ++i) {
    G[i].clear(), R[i].clear(), D[i].clear();
    memset(f[i], 0, sizeof(int) * ((i >> 5) + 1));
  }
}

void dfs(int u) {
  vis.set(u);
  f[u][u >> 5] |= 1 << (u & 31);
  for (int v : D[u])
    if (!(f[u][v >> 5] & (1 << (v & 31)))) {
      if (!vis[v]) dfs(v);
      for (int j = 0; j <= v >> 5; ++j)
        f[u][j] |= f[v][j];
    }
  for (int j = 0; j <= u >> 5; ++j)
    for (int i = 0; i < 32; ++i)
      if (f[u][j] & (1 << i))
        ans[rmp[u]] += w[rmp[(j << 5) + i]];
}

int main() {
  int T, N, M;

  cin >> T;
  for (int _T = 0; _T < T; ++_T) {
    cin >> N >> M;
    if (_T) init(N);
    for (int i = 1; i <= N; ++i)
      cin >> w[i];
    for (int i = 1, u, v; i <= M; ++i) {
      cin >> u >> v;
      G[u].push_back(v);
      R[v].push_back(u);
      ++outd[u];
    }

    for (int i = 1; i <= N; ++i)
      if (!outd[i]) s[++top] = i;
    for (int i = 1, u; top; ) {
      u = s[top--], mp[u] = i, rmp[i] = u, ++i;
      for (int v : R[u])
        if (!--outd[v]) s[++top] = v;
    }
    for (int u = 1; u <= N; ++u)
      for (int v : G[u])
        D[mp[u]].push_back(mp[v]);
    for (int i = 1; i <= N; ++i)
      std::sort(D[i].rbegin(), D[i].rend());
    for (int i = N; i; --i)
      if (!vis[i]) dfs(i);

    for (int i = 1; i <= N; ++i)
      cout << ans[i] << ' ';
    cout << '\n';
  }
  return 0;
}

残缺的字符串

problem

应该叫【模板】带通配符的单模式串匹配（

给定长为 $m$ 的模式串 $A$ 和长为 $n$ 的模式串 $B \ (1 \le m, n \le 3 \times 10^5)$。字符集小写字母和*，*可以匹配任意字符。

求 $B$ 串每个可以完全匹配 $A$ 串的位置。

solution - 单模式串匹配

如果不带通配符，kmp 算法可以在 $O(n + m)$ 的时间内解决这个问题。对于带通配符的问题，我们使用多项式卷积来解决。

首先从不带通配符开始考虑。方便起见，字符串的下标均从 $0$ 开始。

定义完全匹配函数 $P(x) = \sum{i = 0}^{m - 1}{(A_i - B{x - m + i + 1})^2}$，其中平方是为了保证每一位的取值非负；那么若 $P(x) = 0$，则 $B_{x - m + 1 \dots x}$ 与 $A$ 完全匹配。为了制造卷积形式，把 $A$ 串翻转得到串 $S$，那么匹配函数

$\begin{align} P(x) &= \sum_{i = 0}^{m - 1}{(S_{m - i - 1} - B_{x - m + i + 1})^2} \\ &= \sum_{i = 0}^{m - 1}{S_{m - i - 1}^2} + \sum_{i = 0}^{m - 1}{B_{x - m + i + 1}^2} - 2 \sum_{i = 0}^{m - 1}{S_{m - i - 1} B_{x - m + i + 1}} \end{align}$

第一项可以 $O(m)$ 求得，第二项可以 $O(n)$ 预处理前缀和，需要求的只有第三项。发现第三项两个下标的和恰好等于 $x$，因此可以写成 $- 2 \sum_{i + j = x}{S_i B_j}$（想一想，为什么？）

设 $T = \sum{i = 0}^{m - 1} S_i^2, f(x) = \sum{i = 0}^x{Bi^2}, g(x) = \sum{i + j = x}{S_i B_j}$，则有 $P(x) = T + f(x) - f(x - m) - 2 g(x)$。这个直接 FFT 一下就好了。

solution - 带通配符的单模式串匹配

由于出现了通配符，需要定义一个新的匹配函数。设通配符的数值为 0，定义完全匹配函数 $P(x) = \sum{i = 0}^{m - 1}(A_i - B{x - m + i + 1})^2 Ai B{x - m + i + 1}$。把 $A$ 串翻转得到串 $S$，那么匹配函数

$\begin{align} P(x) &= \sum_{i = 0}^{m - 1}{(S_{m - i - 1} - B_{x - m + i + 1})^2 S_{m - i - 1} B_{x - m + i + 1}} \\ &= \sum_{i = 0}^{m - 1}{S_{m - i - 1}^3 B_{x - m + i + 1}} + \sum_{i = 0}^{m - 1}{S_{m - i - 1} B_{x - m + i + 1}^3} - 2 \sum_{i = 0}^{m - 1}{S_{m - i - 1}^2 B_{x - m + i + 1}^2} \end{align}$

发现所有的项的两个下标的和都为 $x$，则有

$P(x) = \sum_{i + j = x}{S_i^3 B_j} + \sum_{i + j = x}{S_i B_j^3} - 2 \sum_{i + j = x}{S_i^2 B_j^2}$

6 次 DFT 和 1 次 IDFT 即刻求出。

code

Complex F[1<<20|1], G[1<<20|1], P[1<<20|1];
int A[1<<20|1], B[1<<20|1], to[1<<20|1];
std::vector<int> ans;

std::vector<int> solve(int M, char *s1, int N, char *s2) {
  int lim, l, tot = 0;

  for (lim = 1, l = -1; lim <= M + N; lim <<= 1, ++l);
  for (int i = 0; i < lim; ++i)
    to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

  for (int i = 0; i < M; ++i)
    A[M - i - 1] = s1[i] == '*' ? 0 : s1[i] & 31;
  for (int i = 0; i < N; ++i)
    B[i] = s2[i] == '*' ? 0 : s2[i] & 31;

  for (int i = 0; i < lim; ++i)
    F[i] = Complex(A[i] * A[i] * A[i], 0), G[i] = Complex(B[i], 0);
  FFT(lim, F, 1), FFT(lim, G, 1);
  for (int i = 0; i < lim; ++i)
    P[i] = P[i] + F[i] * G[i];

  for (int i = 0; i < lim; ++i)
    F[i] = Complex(A[i], 0), G[i] = Complex(B[i] * B[i] * B[i], 0);
  FFT(lim, F, 1), FFT(lim, G, 1);
  for (int i = 0; i < lim; ++i)
    P[i] = P[i] + F[i] * G[i];

  for (int i = 0; i < lim; ++i)
    F[i] = Complex(A[i] * A[i], 0), G[i] = Complex(B[i] * B[i], 0);
  FFT(lim, F, 1), FFT(lim, G, 1);
  for (int i = 0; i < lim; ++i)
    P[i] = P[i] - Complex(2, 0) * F[i] * G[i];

  FFT(lim, P, -1);
  for (int i = M - 1; i < N; ++i)
    if (fabs(P[i].x) < eps) // 本题 eps 取 1e-1 可过
      ans.push_back(i - M + 2);
  return ans;
}

[集训队作业2013]城市规划

problem

求 $n \ (1 \le n \le 130000)$ 个点的有标号简单无向连通图的个数，答案对 $1004535809$ 取模。

solution

设 $f(n)$ 表示有 $n$ 个点的有标号简单无向连通图的个数，$g(n)$ 表示 $n$ 个点的有标号简单无向图的个数，显然有 $g(n) = 2^{\binom n2}$。

一个有标号简单无向图是由若干个连通分量构成的，为了避免重复计数，我们枚举 $1$ 号点所在的连通块的大小（其他点任意连边，不会重复）

$\begin{align} g(n) &= \sum_{i = 1}^n{\binom{n - 1}{i - 1} f(i) g(n - i)} \\ 2^{\binom n2} &= \sum_{i = 1}^n{\binom{n - 1}{i - 1} f(i) 2^{\binom{n - i}2}} \\ \frac{2^{\binom n2}}{(n - 1)!} &= \sum_{i = 1}^n{\frac{f(i)}{(i - 1)!} \times \frac{2^{ \binom{n - i}2}}{(n - i)!}} \end{align}$

定义形式幂级数

$\begin{align} F(x) &= \sum_{n = 0}^\infty \frac{f(n)}{(n - 1)!} x^n \\ G(x) &= \sum_{n = 0}^\infty \frac{2^{ \binom n2}}{n!} x^n \\ H(x) &= \sum_{n = 0}^\infty \frac{2^{ \binom n2}}{(n - 1)!} x^n \end{align}$

那么有

$\begin{align} H(x) & \equiv F(x) G(x) \pmod{x^{n + 1}} \\ F(x) & \equiv H(x) G^{-1}(x) \pmod{x^{n + 1}} \end{align}$

把 $H(x)$ 和 $G(x)$ 的逆卷一块就得到了 $F(x)$，再乘上一个 $(n - 1)!$ 就完事了。总复杂度 $O(n \log n)$。

code

using ll = long long;
using Poly_t = int[1<<18|1];

constexpr int mod = 1004535809, g = 3, invG = 334845270;

inline int Pow(ll b, ll p = mod - 2, ll m = mod);

Poly_t F, G, G_, to;
int *H = F, fac[130003], invF[130003];

void prework(int N) {
  fac[0] = 1;
  for (int i = 1; i <= N; ++i)
    fac[i] = 1ll * fac[i - 1] * i % mod;
  invF[N] = Pow(fac[N]);
  for (int i = N; i; --i)
    invF[i - 1] = 1ll * invF[i] * i % mod;
}

void NTT(int N, int *a, int type);

void PolyInv(int N, int *a, int *b);

int solve(int N) {
  int lim, l;

  std::cin >> N;
  prework(N);
  G[0] = 1;
  for (int i = 1; i <= N; ++i) {
    G[i] = 1ll * Pow(2, 1ll * i * (i - 1) >> 1) * invF[i] % mod;
    H[i] = 1ll * Pow(2, 1ll * i * (i - 1) >> 1) * invF[i - 1] % mod;
  }
  PolyInv(N + 1, G, G_);

  for (lim = 1, l = -1; lim <= N << 1; lim <<= 1, ++l);
  for (int i = 0; i < lim; ++i)
    to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

  NTT(lim, H, 1), NTT(lim, G_, 1);
  for (int i = 0; i < lim; ++i)
    F[i] = 1ll * G_[i] * H[i] % mod;
  NTT(lim, F, -1);
  
  return 1ll * F[N] * fac[N - 1] % mod;
}

CF662C Binary Table

problem

给定一个 $n \ (1 \le n \le 20)$ 行 $m \ (1 \le m \le 10^6)$ 列的 01 矩阵，每次操作可以选定一行或一列并将其中所有元素取反，求若干次操作后可能的最少的 1 的个数。

solution

暴力是这样子的：$2^n$ 枚举每一行是否翻转，每种状态的答案是每一列 0/1 个数的较小值之和。总复杂度 $O(m 2^n)$。

枚举状态 $S$ 表示每一行是否翻转，设 $a_i$ 表示初始矩阵中有多少列状态为 $i$，$b_i$ 表示状态 $i$ 的 0/1 的个数的较小值，$f_i$ 表示枚举到状态 $i$ 时的答案，那么

$f_S = \sum_{i = 0}^{2^n - 1}{a_i b_{i \oplus S}}$

也即

$f_S = \sum_{i \oplus j = S}{a_i b_j}$

发现这是异或卷积的形式，可以用 FWT 在 $O(n \log n)$ 的时间内求出 $f_{1 \cdots 2^n}$，取最小值就是答案。

code

FWT 实现看多项式笔记。本题不需要也不应该取模。

using ll = long long;
using Poly_t = ll[1<<20|1];

int z[100003];
Poly_t a, b;

int main() {
  int N, N_, M;
  ll* f = a;

  cin >> N >> M;
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < M; ++j)
      if (cin.get())
        z[j] |= 1 << i;
  for (int i = 0; i < M; ++i)
    ++a[z[i]];
  N_ = N, N = 1 << N;
  for (int i = 0; i < N; ++i)
    b[i] = b[i >> 1] + (i & 1);
  for (int i = 0; i < N; ++i)
    b[i] = std::min(b[i], N_ - b[i]);
  
  XOR(N, a, 1), XOR(N, b, 1);
  for (int i = 0; i < N; ++i)
    f[i] = a[i] * b[i];
  XOR(N, f, 0);

  std::cout << *std::min_element(f, f + N);
  return 0;
}

CF850E Random Elections

problem

A, B, C 三人比赛，$2^n \ (1 \le n \le 20)$ 个人投票。每个人有一个确定的顺序，他们把票投给顺序靠前的人。

有一个函数 $f(S)$，接受一个长度 $n = 2^k$ 的二进制串 $S$ 并返回 true/false。保证 $f(\neg S_1, \neg S_2, \cdots, \neg S_n) = \neg f(S_1, S_2, \cdots, S_n)$。

比赛结果由 $n$ 个人投票决定，给定 $S$，$S_i$ 表示 $f(i)$ 的取值。询问单循环的三场比赛满足有同一个人连胜两场的方案数。

solution

假设 A 连胜两场，答案乘以 $3$ 即可。

设某人两场比赛的决策是 $(x_1, x_2)$，那么可以得到他的决策序列：

$(0, 0) \to \text{BCA, CBA}$
$(0, 1) \to \text{BAC}$
$(1, 0) \to \text{CAB}$
$(1, 1) \to \text{ABC, ACB}$

发现如果 $x_1 \oplus x_2 = 0$ 就可以对方案数造成 $2$ 的贡献。所有位压在一起那么方案就是 $2^{n - \operatorname{popcount}(S_1 \oplus S_2)}$。

直接枚举是 $O(4^n)$ 的，可以 FWT 预处理出异或值等于 $1 \cdots 2^n$ 的数对个数再统计贡献。

code

constexpr int mod = 1'000'000'007, inv2 = 500'000'004;

int solve(int N_, int* F) { // F[0...2^N_)
  int N = 1 << N_, ans = 0, *popcnt = new int[N];

  XOR(N, F, 1);
  for (int i = 0; i < N; ++i)
    F[i] = 1ll * F[i] * F[i] % mod;
  XOR(N, F, inv2);

  for (int i = 0; i < N; ++i) {
    popcnt[i] = popcnt[i >> 1] + (i & 1);
    ans = (ans + 1ll * F[i] * (1 << (N_ - popcnt[i]))) % mod;
  }
  return 3ll * ans % mod;
}

P5816 [CQOI2010]内部白点

problem

初始时网格图里有 $n \ (1 \le n \le 10^5)$ 个黑点，所有其他整数点都是白色的。每秒钟所有内部白点同时变黑，直到不存在内部白点为止。计算最终网格图中的黑点个数或判断变化不会结束。

内部白点的定义：一个白色的整点是内部白点当且仅当其所在的水平线的左边和右边各至少有一个黑点，且所在竖直线的上边和下边各至少有一个黑点。形式化地，点 $P(x, y)$ 是黑点当且仅当存在黑点 $(x_1, y), (x_2, y), (x, y_1), (x, y_2)$ 满足 $x_1 < x < x_2, y_1 < x < y_2$。

solution

👴又回来整 DS 🌶！

首先这题变化是一定会结束的，并且会在一次变化内结束~~（好像很显然的样子）~~，因为新产生的黑点能够产生的黑点一定是原有黑点能产生的黑点的子集。

我们对于每个 $x$ 坐标，记录 $y$ 坐标最小和最大的点，对于 $y$ 最小/最大的点打标记 $1 / -1$，再对每对 $x$ 坐标相同的点连一条水平线段。

这样从下向上扫描，对于每个标记修改相应位置的点值；对于每条线段统计区间和，即可得到答案。

code

这个做法除去 ds 外时间常数并不好；空间常数会好一些；在碰到重复点时会暴毙，好在这题保证了没有重复点。

struct Point {
  int x, y;
} z[100003];

bool cmpx(Point lhs, Point rhs) {
  return std::tie(lhs.x, lhs.y) < std::tie(rhs.x, rhs.y);
}

bool cmpy(Point lhs, Point rhs) {
  return std::tie(lhs.y, lhs.x) < std::tie(rhs.y, rhs.x);
}

struct Line {
  int l, r, y, k;

  bool operator<(Line rhs) {
    return std::tie(y, k) < std::tie(rhs.y, rhs.k);
  }
} e[300003];

int N_, z_[100003];

namespace BIT {
  int t[100003];
  void Modify(int p, int v);
  int Query(int p);
} // BIT

int main() {
  int N, M = 0, ans = 0;

  cin >> N;
  for (int i = 1; i <= N; ++i) {
    cin >> z[i].x >> z[i].y;
    z_[i] = z[i].x;
  }

  std::sort(z_ + 1, z_ + N + 1);
  N_ = std::unique(z_ + 1, z_ + N + 1) - z_ - 1;
  for (int i = 1; i <= N; ++i)
    z[i].x = std::lower_bound(z_ + 1, z_ + N_ + 1, z[i].x) - z_;

  std::sort(z + 1, z + N + 1, cmpx);
  for (int i = 1; i < N; ++i)
    if (z[i].x == z[i + 1].x) {
      e[++M] = {z[i].x, 0, z[i].y, 1};
      e[++M] = {z[i + 1].x, 0, z[i + 1].y, -1};
    }
  std::sort(z + 1, z + N + 1, cmpy);
  for (int i = 1; i < N; ++i)
    if (z[i].y == z[i + 1].y)
      e[++M] = {z[i].x, z[i + 1].x, z[i].y, 0};
  std::sort(e + 1, e + M + 1);

  for (int i = 1; i <= M; ++i)
    if (e[i].k) BIT::Modify(e[i].l, e[i].k);
    else ans += BIT::Query(e[i].r - 1) - BIT::Query(e[i].l);

  std::cout << ans + N << '\n';
  return 0;
}

总结

写于 2020/4/30 22:46

这个月结束了，可以说是放假来最不颓的一个月（其实颓得也蛮多的）

先整了几道题，又重学了数论、从零开始学了多项式，统共码了 2w+ 字的博文。

基础算法和思维比较薄弱，希望下个月能加强下。

下个月这时候说不定就要退役了。

希望能够不辜负自己吧。