2025 7 月闲话合集

这一部分都是我 7 月份（中后旬）写的题的总结，如果对你有用，真是倍感荣幸。为了轻松一点，每四道题目我将放一段歌词，以供欣赏。

[Medium] P10664 BZOJ3328 PYXFIB

给定 $n,k,p$，计算：

$$\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor} \binom{n}{ik} \mathrm{fib}(i) \pmod{p}$$

其中 $\mathrm{fib}(i)$ 表示 Fibonacci 数列，定义为 $\mathrm{fib}(0)=\color{red}\mathbf{1}\color{default}, \mathrm{fib}(1)=1, \mathrm{fib}(n)=\mathrm{fib}(n-1)+\mathrm{fib}(n-2)$。

$T$ 组数据，$1\leq T\leq 20,1\leq n\leq 10^{18},1\leq k\leq 2\times 10^4,1\leq p\leq 10^9,p\in\mathbb{P},k\mid p-1$。

那么我们可以改为求：

$$\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ik}\left(\begin{bmatrix}1&1\end{bmatrix}\cdot \begin{bmatrix}1&1\\1&0\end{bmatrix}^i\right)=\begin{bmatrix}1&1\end{bmatrix}\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor} \binom{n}{ik}\begin{bmatrix}1&1\\1&0\end{bmatrix}^i$$

接下来就是平凡的推式子了：

$$\begin{aligned} &\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor} \binom{n}{ik}\begin{bmatrix}1&1\\1&0\end{bmatrix}^i=\sum_{i=0}^{n}\binom{n}{i}\begin{bmatrix}1&1\\1&0\end{bmatrix}^i[k\mid i]\\ =&\frac{1}{k}\sum_{i=0}^{n}\binom{n}{i}\begin{bmatrix}1&1\\1&0\end{bmatrix}^i\sum_{j=0}^{k-1} \omega_k^{ij}=\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^{n}\binom{n}{i}\left(\omega_k^{j}\begin{bmatrix}1&1\\1&0\end{bmatrix}\right)^i\\ =&\frac{1}{k}\sum_{j=0}^{k-1}\left(\begin{bmatrix}1&0\\0&1\end{bmatrix}+\omega_k^{j}\begin{bmatrix}1&1\\1&0\end{bmatrix}\right)^n \end{aligned}$$

使用矩阵快速幂实现，时间复杂度 $O(k\log n)$。最后我们还需要找出有限域 $\mathbb{Z}/p\mathbb{Z}$ 中的本原单位根 $\omega_k$。只需要找到一个模 $p$ 意义下原根 $g$，就可以计算 $g^{\frac{p-1}{k}}$ 并用它替代 $k$ 次本原单位根。寻找原根的时间复杂度是 $O(\sqrt{p}+g_{\min}(p)\omega(p-1)\log p)$。

碎碎念

• 根据 Hardy–Ramanujan Theorem，$\omega(p-1)$ 平均规模为 $O(\log \log p)$，值得注意点是，最大值为 $O(\frac{\log p}{\log\log p})$。
• 如果你承认 GRH，那么根据 Shoup, V. (1992). Searching for primitive roots in finite fields. Mathematics of Computation, 58(197), 369-380 的结果，$g_{\min}(p)\sim o(\log^6 p)$。
• Kevin J. McGown, Jonathan P. Sorenson Computation of the least primitive root 给出了当 $1\leq p\leq 10^9$ 时，$g_{\min}(p)\leq 151$。

实现

#include <bits/stdc++.h>
#define _all(x) std::begin(x), std::end(x)
using std::cin, std::cout;

int mod;
int A(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int S(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int M(int x, int y){ return 1ll * x * y % mod; }
int P(int x, int y){ int ans = 1; for(;y;x=M(x,x),y>>=1) if(y & 1) ans = M(ans, x); return ans; }

template<int N>
struct matrix {
    int a[N][N], n, m;
    matrix(int n, int m) : n(n), m(m) { memset(a, 0, sizeof(a)); }
    int* operator[](int i){ return a[i]; }
    const int* operator[](int i) const { return a[i]; }
    matrix operator+(const matrix& b){
        matrix c(n, m);
        for(int i=0;i<n;i++) for(int j=0;j<m;j++) c[i][j] = A(a[i][j], b[i][j]);
        return c;
    }
    matrix operator*(int b){
        matrix c(n, m);
        for(int i=0;i<n;i++) for(int j=0;j<m;j++) c[i][j] = M(a[i][j], b);
        return c;
    }
    matrix operator*(const matrix& b){
        matrix c(n, b.m);
        for(int i=0;i<n;i++){
            for(int k=0;k<m;k++) for(int j=0;j<b.m;j++) c[i][j] = A(c[i][j], M(a[i][k], b[k][j]));
        }
        return c;
    }
    static matrix idt(int n){ matrix I(n, n); for(int i=0;i<n;i++) I[i][i] = 1; return I; }
    template<class T> matrix operator^(T x){
        matrix ans = idt(n), a = *this;
        for(;x;x>>=1,a=a*a) if(x & 1) ans = ans * a;
        return ans;
    }
};

namespace solution {
    using i64 = int64_t;
    i64 n; int k;
    std::basic_string<int> mp;
    using matrix = ::matrix<2>;
    void solve(){
        cin >> n >> k >> mod; int tmp = mod - 1, g = 2;
        for(int i=2;i*i<=tmp;i++) if(!(tmp % i)){ mp += i; while(!(tmp % i)) tmp /= i; }
        if(tmp > 1) mp += tmp;
        for(;;g++) if(std::all_of(_all(mp), [=](int x){ return P(g, (mod - 1) / x) != 1; })) break;
        matrix ret(2, 2), idt = matrix::idt(2), mat(2, 2);
        mat[0][0] = mat[0][1] = mat[1][0] = 1;
        int omega = P(g, (mod - 1) / k);
        for(int i=0;i<k;i++){
            matrix tmp = idt + mat * P(omega, i);
            ret = ret + (tmp ^ n);
        }
        cout << M(A(ret[0][1], ret[1][1]), P(k, mod - 2)) << '\n';
    }
    void clear(){ mp.clear(); }
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_10664_BZOJ_3328_PYXFIB.cpp] Think twice, code once.
// Written by xiezheyuan

[Hard] P10084 [GDKOI2024 提高组] 计算

令 $F(x,a,b)=\gcd(x^a-1,x^b-1)+1$，特别地，当 $ab=0$ 时规定 $F(x,a,b)=0$。

给定 $m,a,b,c,d$，令 $L=F(m,a,b)+1,R=F(m,c,d)$。

计算有多少个 $S\subseteq [L,R]\cap\mathbb{Z}$ 满足 $m\mid \sum\limits_{i\in S} i$，答案对 $\mathrm{mod}=998,244,353$ 取模。

$T$ 组数据，$1\leq T\leq 10^4,1\leq m\leq 10^7,0\leq a,b,c,d\leq 10^3,0\leq L<R\leq 10^{18}$。

引理 1：$F(x,a,b)=x^{\gcd(a,b)}$（$ab\neq 0$）。

证明：钦定 $a\geq b$，则有：

$$\begin{aligned} F(x,a,b)&=\gcd(x^a-1,x^b-1)+1=\gcd(x^a-1-(x^b-1),x^b)+1\\ &=\gcd(x^b(x^{a-b}-1),x^b-1)=\gcd(x^{a-b}-1,x^b-1) \end{aligned}$$

对指数数学归纳，则立即有：

$$F(x,a,b)=x^{\gcd(a,b)}-1+1=x^{\gcd(a,b)}$$

Q.E.D.

那么就有：

$$L = m^{\gcd(a,b)}+1,R=m^{\gcd(c,d)}$$

则 $m\mid (L-1),m\mid R$。不妨将所有数全部减去 $L-1$，并令 $n=R-L+1$，则将原本的集合 $[L,R]\cap\mathbb{N}$ 归约到 $[1,n]\cap\mathbb{N}$。

考虑如何刻画子集和，常见的方法就是生成函数，根据 OGF 理论，我们很容易写出答案的一个表达式：

$$\sum_{i\geq 0} [m\mid i]\cdot [z^i]\prod_{k=1}^{n}(1+z^k)$$

后面的 $n$ 个二项式的乘积目前是较难处理的，不妨先令 $F(z)=\prod\limits_{k=1}^{n}(1+z^k)$，将答案改写为 $\sum\limits_{i\geq 0}[m\mid i]\cdot [z^i]F(z)$​。

不妨进行单位根反演：

$$\begin{aligned} &\sum_{i\geq 0}[m\mid i]\cdot [z^i]F(z)=\frac{1}{m}\sum_{i\geq 0}[z^i]F(z)\sum_{k=0}^{m-1}\omega_m^{ik}\\ =&\frac{1}{m}\sum_{k=0}^{m-1}\sum_{i\geq 0}[z^i]F(z)\cdot \left(\omega_m^k\right)^i=\frac{1}{m}\sum_{k=0}^{m-1} F(\omega_m^k) \end{aligned}$$

因此我们得到一个经典引理，这是单位根反演公式的直接推论：

引理 2：

$$\sum_{i\geq 0}[m\mid i]\cdot [z^i]F(z)=\frac{1}{m}\sum_{k=0}^{m-1} F(\omega_m^k)$$

那么我们现在就需要研究 $F(\omega_m^k)$ 的取值，将 $F(z)=\prod\limits_{k=1}^{n}(1+z^k)$ 代入：

$$F(\omega_m^k)=\prod_{i=1}^{n}(1+\omega_m^{ik})$$

由于 $m\mid (L-1),m\mid R$，因此 $m\mid R-(L-1)$，即 $m\mid n$，又因为根据本原单位根的性质，$\forall i,j,k\in\mathbb{Z}^+,\omega_m^{im+k}=\omega_m^{jm+k}$，所以可以改写：

$$F(\omega_m^k)=\prod_{i=1}^{n}(1+\omega_m^{ik})=\left(\prod_{i=0}^{m-1}(1+\omega_m^{ik})\right)^{\frac{n}{m}}$$

取 $p=\gcd(m,k)$，根据本原单位根的另一个性质 $\omega_{mk}^{ak}=\omega_m^a$，得到：

$$\begin{aligned} &\sqrt[\frac{n}{m}\cdot p]{F(\omega_m^k)}=\prod_{i=0}^{\frac{m}{p}-1}(1+\omega_{\frac{m}{p}}^{i\frac{k}{p}})=\prod_{i=0}^{\frac{m}{p}-1}(1+\omega_{\frac{m}{p}}^i)\\ =&(-1)^{\frac{m}{p}}\prod_{i=0}^{\frac{m}{p}-1}((-1)-\omega_{\frac{m}{p}}^i)=(-1)^{\frac{m}{p}}((-1)^{\frac{m}{p}}-1)=2\left[2\nmid \frac{m}{p}\right] \end{aligned}$$

倒数第二步用到的是根据单位根的定义的因式分解 $z^n-1=\prod_{i=0}^{n-1} (z-\omega_{n}^i)$。将结果代入：

$$\begin{aligned} &\frac{1}{m}\sum_{k=0}^{m-1} F(\omega_m^k)=\frac{1}{m}\sum_{k=0}^{m-1}\left[2\nmid \frac{m}{\gcd(m,k)}\right]2^{\frac{n\gcd(m,k)}{m}}\\ =&\frac{1}{m}\sum_{d\mid m}\left[2\nmid \frac{m}{d}\right]2^{\frac{nd}{m}}\sum_{k=0}^{m-1}[\gcd(m,k)=d]=\frac{1}{m}\sum_{d\mid m}2^{\frac{nd}{m}}\left[2\nmid \frac{m}{d}\right]\sum_{k=1}^{\frac{m}{d}}\left[\gcd\left(\frac{m}{d},k\right)=1\right]\\ =&\frac{1}{m}\sum_{d\mid m}2^{\frac{nd}{m}}\left[2\nmid \frac{m}{d}\right]\varphi\left(\frac{m}{d}\right)=\frac{1}{m}\sum_{d\mid m}2^{\frac{n}{d}}[2\nmid d]\varphi(d) \end{aligned}$$

$O(m)$ 线性筛出 $\varphi$，然后每组数据暴力枚举 $m$ 的因子即可，时间复杂度 $O(m+T\sqrt{m}\log \mathrm{mod})$，勉强能过。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

constexpr int mod = 998244353;
int A(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int S(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int M(int x, int y){ return 1ll * x * y % mod; }
int P(int x, int y){ int ans = 1; for(;y;x=M(x,x),y>>=1) if(y & 1) ans = M(ans, x); return ans; }
int I(int x){ return P(x, mod - 2); }
template<class T> int Fit(T x){ return (x % mod + mod) % mod; }

const int N = 1e7 + 5;
int pri[N], phi[N], tot;
bool vis[N];

void sieve(int n){
    phi[1] = 1;
    for(int i=2;i<=n;i++){
        if(!vis[i]) pri[++tot] = i, phi[i] = i - 1;
        for(int j=1;j<=tot&&1ll*i*pri[j]<=n;j++){
            vis[i * pri[j]] = 1;
            if(!(i % pri[j])){
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
    }
}

namespace solution {
    using i64 = long long;
    int m, a, b, c, d;
    std::basic_string<int> ds;

    i64 naive_pow(int x, int y){
        i64 ret = 1;
        for(int i=0;i<y;i++) ret = ret * x;
        return ret;
    }

    i64 F(int m, int a, int b){ return (!a || !b) ? 0 : naive_pow(m, std::gcd(a, b)); }

    void solve(){
        cin >> m >> a >> b >> c >> d;
        i64 L = F(m, a, b) + 1, R = F(m, c, d), n = R - L + 1;
        for(int i=1;i*i<=m;i++){
            if(m % i) continue;
            ds += i;
            if(i * i != m) ds += m / i;
        }
        int ans = 0;
        for(int x : ds) if(x & 1) ans = A(ans, M(phi[x], P(2, (n / x) % (mod - 1))));
        cout << M(ans, I(m)) << '\n';
    }

    void clear(){ ds.clear(); }
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    sieve(1e7);
    int T; cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_10084_GDKOI_2024_提高组_计算.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] P4068 [SDOI2016] 数字配对

给定 $n$ 种数字，第 $i$ 种数字的值为 $a_i$，有 $b_i$ 个，权值为 $c_i$。

如果存在两种数字 $a_i$ 和 $a_j$，满足 $a_i$ 是 $a_j$ 的倍数，且 $\frac{a_i}{a_j}$ 是一个质数，那么这两种数字可以进行配对，并获得 $c_i \times c_j$ 的价值。

每个数字只能参与一次配对，也可以不参与配对。在总价值不小于 $0$ 的前提下，求最多可以进行多少次配对。

$1\leq n\leq 200,1\leq a_i\leq 10^9$。

注意到这个形式与匹配非常相似，考虑将其套上二分图匹配模型，采用网络流技术解决。

注意到若记 $\Omega(x)$ 表示 $x$ 的质因子数量，则 $a_i,a_j$ 可以匹配的一个必要条件是 $\Omega(a_i),\Omega(a_j)$ 奇偶性不同，因此按照能否匹配连边，就可以得到一个二分图。

具体地，对于可以匹配的 $a_i,a_j$，连边 $(i,j,\min(b_i,b_j),c_i c_j)$。然后我们需要限制费用 $\geq 0$ 的情况下求出最大流。这个我们可以通过二分最大流量，然后建立一个超级源点 $S'$，向原来的源点 $S$ 连边 $(S',S,\mathrm{flow},0)$ 限制总流量来实现。

时间复杂度 $O(n\sqrt{a_i}+\mathrm{mcmf}(n,n^2)\log b_i)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int, class Cost = int>
struct MCMF {
    struct edge {
        int nxt, to;
        Cap cap;
        Cost cost;
    } g[M << 1];

    int head[N], ec, S, T, cur[N];
    Cost dis[N], mincost;
    bool vis[N];

    MCMF() : ec(1) {}

    void add(int u, int v, Cap cap, Cost cost){
        g[++ec] = {head[u], v, cap, cost}, head[u] = ec;
        g[++ec] = {head[v], u, 0, -cost}, head[v] = ec;
    }

    bool spfa(){
        for(int i=0;i<N;i++) dis[i] = std::numeric_limits<Cost>::max();
        Cost inf = std::numeric_limits<Cost>::max();
        std::queue<int> q; q.push(S);
        dis[S] = 0;
        while(!q.empty()){
            int u = q.front();
            q.pop(), vis[u] = 0;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(g[i].cap && dis[u] + g[i].cost < dis[v]){
                    dis[v] = dis[u] + g[i].cost;
                    if(!vis[v]) q.push(v), vis[v] = 1;
                }
            }
        }
        return dis[T] < inf;
    }

    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        vis[u] = 1;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(!vis[v] && g[i].cap && dis[u] + g[i].cost == dis[v]){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t){
                    mincost += t * g[i].cost;
                    g[i].cap -= t, g[i ^ 1].cap += t;
                    return vis[u] = 0, t;
                }
            }
        }
        return vis[u] = 0, 0;
    }

    std::pair<Cost, Cap> operator()(int s, int t){
        S = s, T = t, mincost = 0;
        Cap maxflow = 0, tmp;
        while(spfa()){
            std::memcpy(cur, head, sizeof(head));
            while((tmp = dfs(s, std::numeric_limits<Cap>::max()))) maxflow += tmp;
        }
        return {mincost, maxflow};
    }
};

namespace solution {
    const int N = 205;
    MCMF<N, N * N, int, long long> mcmf;
    int n, a[N], b[N], c[N], f[N];

    std::pair<long long, int> calc(long long penalty){
        for(int i=2;i<=mcmf.ec;i+=2){
            mcmf.g[i].cap += mcmf.g[i ^ 1].cap;
            mcmf.g[i ^ 1].cap = 0;
        }
        mcmf.g[mcmf.ec - 1].cap = penalty;
        auto [cost, flow] = mcmf(0, n + 1);
        return {cost, flow};
    }

    void solve(){
        cin >> n;
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=n;i++) cin >> b[i];
        for(int i=1;i<=n;i++) cin >> c[i];
        for(int i=1;i<=n;i++){
            int tmp = a[i];
            for(int j=2;j*j<=tmp;j++){
                if(tmp % j) continue;
                while(!(tmp % j)) f[i]++, tmp /= j;
            }
            if(tmp > 1) f[i]++;
        }
        int S = n + 2, T = n + 1;
        for(int i=1;i<=n;i++){
            if(f[i] & 1) mcmf.add(S, i, b[i], 0);
            else mcmf.add(i, T, b[i], 0);
        }
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                if(!(a[i] % a[j]) && f[i] == f[j] + 1){
                    if(f[i] & 1) mcmf.add(i, j, INT_MAX, -1ll * c[i] * c[j]);
                    else mcmf.add(j, i, INT_MAX, -1ll * c[i] * c[j]);
                }
            }
        }
        mcmf.add(0, S, INT_MAX, 0);
        int L = 0, R = 1e9;
        while(L < R){
            int mid = (L + R + 1) >> 1;
            if(-calc(mid).first >= 0) L = mid;
            else R = mid - 1;
        }
        cout << calc(L).second << '\n';
    }

    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_4068_SDOI_2016_数字配对.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium-Hard] P5331 [SNOI2019] 通信3000 ms

$n$ 个排成一列的哨站要进行通信。第 $i$ 个哨站的频段为 $a_i$。每个哨站 $i$ 需要选择以下二者之一：

1. 直接连接到控制中心，代价为 $W$；
2. 连接到前面的某个哨站 $j$ ($j<i$)，代价为 $|a_i-a_j|$。
   每个哨站只能被后面的至多一个哨站连接。

请你求出最小可能的代价和，$1\leq n\leq 10^3$。

有一个非常简单的网络流建模：

• 对于第 $i$ 个哨站，建立点 $u_i,v_i$，连边 $(S,u_i,1,0),(v_i,T,1,0)$ 表示每个哨站只能被使用一次。
• 连边 $(S,v_i,1,W)$ 表示哨站 $i$ 可以直接连接到控制中心。
• 对于 $i>j$，连边 $(u_i,v_j,1,|a_i-a_j|)$ 表示哨站 $i$ 可以连接到哨站 $j$。

然后跑最小费用最大流即可。不过这样做需要连接 $O(n^2)$ 条边，无法接受。并且让人难受的是，边权是绝对值，众所周知，绝对值和 $\max$ 很兼容，但是和 $\min$ 不太兼容，所以没有办法拆绝对值后直接做，只能通过大小关系把绝对值强行拆掉。

我们发现这个东西比较像二维偏序，可以考虑数据结构优化建图。这里我们使用 cdq 分治优化建图。具体地，我们将 $u_i$ 拆成三个点 $u_{i1},u_{i2},u_{i3}$，连边 $(S,u_{i1},1,0),(u_{i1},u_{i2},1,a_i),(u_{i1},u_{i3},1,-a_i)$ 用于强行拆绝对值，$v_i$ 也是类似的。这样我们只需要在对应到 $u_{ix},v_{iy}$ 连边就可以体现绝对值。如何连边？我们 cdq 分治，考虑右侧区间对左侧区间的贡献，对于右侧区间的每一个位置，用双指针找出左侧区间第一个大于它的位置，然后需要向一个左侧区间前缀 / 后缀连边。这可以使用前缀和 / 后缀和优化建图来实现（具体来说就是拉一条毛毛虫那样的链出来）。这样我们就将连边数缩减到了 $O(n\log n)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int, class Cost = int>
struct MCMF {
    struct edge {
        int nxt, to;
        Cap cap;
        Cost cost;
    } g[M << 1];

    int head[N], ec, S, T, cur[N];
    Cost dis[N], mincost, phi[N];
    bool vis[N];

    MCMF() : ec(1) {}

    void add(int u, int v, Cap cap, Cost cost){
        g[++ec] = {head[u], v, cap, cost}, head[u] = ec;
        g[++ec] = {head[v], u, 0, -cost}, head[v] = ec;
    }

    bool spfa(){
        for(int i=0;i<N;i++) phi[i] = std::numeric_limits<Cost>::max();
        Cost inf = std::numeric_limits<Cost>::max();
        std::queue<int> q; q.push(S);
        phi[S] = 0;
        while(!q.empty()){
            int u = q.front();
            q.pop(), vis[u] = 0;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(g[i].cap && phi[u] + g[i].cost < phi[v]){
                    phi[v] = phi[u] + g[i].cost;
                    if(!vis[v]) q.push(v), vis[v] = 1;
                }
            }
        }
        return phi[T] < inf;
    }

    Cost cost(int id){ return phi[g[id ^ 1].to] - phi[g[id].to] + g[id].cost; }

    bool dijkstra(){
        for(int i=0;i<N;i++) dis[i] = std::numeric_limits<Cost>::max();
        using pii = std::pair<Cost, int>;
        std::priority_queue<pii, std::vector<pii>, std::greater<pii>> q;
        q.emplace(dis[S] = 0, S);
        while(!q.empty()){
            int u = q.top().second;
            q.pop(), vis[u] = 1;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(!vis[v] && g[i].cap && dis[u] + cost(i) < dis[v]){
                    dis[v] = dis[u] + cost(i), q.emplace(dis[v], v);
                }
            }
        }
        bool ans = vis[T];
        std::fill(vis, vis + N, 0);
        return ans;
    }

    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        vis[u] = 1;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(!vis[v] && g[i].cap && dis[u] + cost(i) == dis[v]){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t){
                    mincost += t * g[i].cost;
                    g[i].cap -= t, g[i ^ 1].cap += t;
                    return vis[u] = 0, t;
                }
            }
        }
        return vis[u] = 0, 0;
    }

    std::pair<Cost, Cap> operator()(int s, int t){
        S = s, T = t, mincost = 0;
        Cap maxflow = 0, tmp;
        spfa();
        while(dijkstra()){
            memcpy(cur, head, sizeof(head));
            while((tmp = dfs(s, std::numeric_limits<Cap>::max()))) maxflow += tmp;
            for(int i=0;i<N;i++){
                if(dis[i] != std::numeric_limits<Cost>::max()) phi[i] += dis[i];
            }
        }
        return {mincost, maxflow};
    }
};

namespace solution {
    const int N = 1005;
    int n, W, a[N]; MCMF<int(5e5), int(1e6), int, long long> mcmf;
    int node(int x, int y){ return (x - 1) * 6 + y + 1; }
    int tot;

    template<class T, class Cmp = std::greater<>>
    using heap = std::priority_queue<T, std::vector<T>, Cmp>;

    void cdq(int l, int r){
        if(l >= r) return;
        int mid = (l + r) >> 1;
        cdq(l, mid), cdq(mid + 1, r);
        heap<std::pair<int,int>> h;
        for(int i=l;i<=mid;i++) h.push({a[i], i});
        std::vector<int> buf(mid - l + 1), bufid(mid - l + 1); int bt = tot;
        for(int i=0;i<=mid-l;i++){
            int x = h.top().second; h.pop();
            if(i) mcmf.add(tot, tot + 1, INT_MAX, 0);
            mcmf.add(node(x, 3), ++tot, 1, 0);
            buf[i] = a[x], bufid[i] = x;
        }
        int bt2 = tot;
        for(int i=mid;i>=l;i--){
            int x = bufid[i - l];
            if(i != mid) mcmf.add(tot, tot + 1, INT_MAX, 0);
            mcmf.add(node(x, 2), ++tot, 1, 0);
        }
        for(int i=mid+1;i<=r;i++){
            int rnk = std::lower_bound(buf.begin(), buf.end(), a[i]) - buf.begin();
            if(rnk) mcmf.add(bt + rnk, node(i, 4), 1, 0);
            if(rnk != buf.size()) mcmf.add(bt2 + buf.size() - rnk, node(i, 5), 1, 0);
        }
    }

    void solve(){
        cin >> n >> W; int S = 0, T = 1;
        for(int i=1;i<=n;i++){
            cin >> a[i];
            mcmf.add(S, node(i, 1), 1, 0);
            mcmf.add(node(i, 1), node(i, 2), 1, a[i]);
            mcmf.add(node(i, 1), node(i, 3), 1, -a[i]);
            mcmf.add(node(i, 4), node(i, 6), 1, a[i]);
            mcmf.add(node(i, 5), node(i, 6), 1, -a[i]);
            mcmf.add(node(i, 6), T, 1, 0), mcmf.add(S, node(i, 6), 1, W);
        }
        tot = node(n, 6);
        cdq(1, n);
        cout << mcmf(S, T).first << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_5331_SNOI_2019_通信.cpp] Think twice, code once.
// Written by xiezheyuan

Have a Break Round I

太逼真的恶耗缠着 太坚稳的隔膜围着 太崩紧的脑部何事又烦恼着

你可知这个地球上 有千亿种怪异模样 快出走观照着迷幻 告别惆怅

雪国里 那雪兽 渗奇香 斜角里 发觉有 麻瓜演奏一场

暗笑我 这岁数 爱幻想 你却变作顽石 心不知痒

陪我迷上荒谬 越过幻变山丘 拾块木碎伸手学施小魔咒

踏遍无数星斗 云雨外喝清酒 忘记俗世几多繁琐的争斗

—— 《特伦斯梦游仙境》 林家谦 SEVEN

[Medium-Hard] P6967 [NEERC 2016] Delight for a Cat

每小时，猫可以选择睡觉或吃东西。猫不能在同一小时内同时进行这两种活动，并且猫在整小时内只能进行其中一种活动。

对于接下来的 $n$ 小时，已知猫在第 $i$ 小时睡觉的快乐值为 $s_i$，吃东西的快乐值为 $e_i$。

你还知道一个整数时间段 $k$。在每 $k$ 个连续的小时中，至少有 $m_{s}$ 小时猫在睡觉，至少有 $m_{e}$ 小时猫在吃东西。因此，有正好 $n - k + 1$ 个 $k$ 小时的时间段需要满足这个条件。

求猫在接下来的 $n$ 小时内所能获得的最大总快乐值，$1\leq n\leq 10^3$。

Reference：ix35 的题解。

考虑一个更加普适的问题：

区间覆盖问题：给定 $[1,m]$ 中的 $n$ 个区间，每个区间选择一次的代价是 $w_i$，最多可以选 $c_i$ 次。要求对于任意点 $j$ 被覆盖的次数位于区间 $[l_j,r_j]$ 中。求最小 / 最大代价。

不妨连边 $(1,2),(2,3),\cdots,(m,m+1)$，然后用 $(i,i+1)$ 的边的流量来刻画 $i$ 被覆盖的次数。

考虑刻画区间覆盖，不妨对于区间 $[L_i,R_i]$，连边 $(L_i,R_i+1,c_i,w_i)$，当有 $1$ 单位流量流过该边时，表示选择该区间一次。那么 $(i,i+1)$ 的容量设计就很明显了，取一个极大值 $N$，连边 $(i,i+1,[N-r_i,N-l_i],0)$（这里流量有了上下界）。然后只需要连边 $(S,1,\infty,0),(m+1,T,\infty,0)$ 跑费用流就是答案。

回到本题。我们改为钦定猫总是在睡觉，然后选择一些小时改成吃饭。那么连续 $k$ 个小时猫吃饭的次数要介于 $[m_e,k-m_s]$。接下来考虑如何选择吃饭，我们可以将第 $i$ 小时吃饭改成对区间 $[i,i+k-1]$ 做一次覆盖，每个点的覆盖次数介于 $[m_e,k-m_s]$ 那么就转成了区间覆盖问题。

取 $N=k-m_s$，这样流量下界就变成了 $0$，就可以用普通的最小费用最大流来做了。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int, class Cost = int>
struct MCMF {
    struct edge {
        int nxt, to;
        Cap cap;
        Cost cost;
    } g[M << 1];

    int head[N], ec, S, T, cur[N];
    Cost dis[N], mincost;
    bool vis[N];

    MCMF() : ec(1) {}

    void add(int u, int v, Cap cap, Cost cost){
        g[++ec] = {head[u], v, cap, cost}, head[u] = ec;
        g[++ec] = {head[v], u, 0, -cost}, head[v] = ec;
    }

    bool spfa(){
        for(int i=0;i<N;i++) dis[i] = std::numeric_limits<Cost>::max();
        Cost inf = std::numeric_limits<Cost>::max();
        std::queue<int> q; q.push(S);
        dis[S] = 0;
        while(!q.empty()){
            int u = q.front();
            q.pop(), vis[u] = 0;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(g[i].cap && dis[u] + g[i].cost < dis[v]){
                    dis[v] = dis[u] + g[i].cost;
                    if(!vis[v]) q.push(v), vis[v] = 1;
                }
            }
        }
        return dis[T] < inf;
    }

    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        vis[u] = 1;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(!vis[v] && g[i].cap && dis[u] + g[i].cost == dis[v]){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t){
                    mincost += t * g[i].cost;
                    g[i].cap -= t, g[i ^ 1].cap += t;
                    return vis[u] = 0, t;
                }
            }
        }
        return vis[u] = 0, 0;
    }

    std::pair<Cost, Cap> operator()(int s, int t){
        S = s, T = t, mincost = 0;
        Cap maxflow = 0, tmp;
        while(spfa()){
            std::memcpy(cur, head, sizeof(head));
            while((tmp = dfs(s, std::numeric_limits<Cap>::max()))) maxflow += tmp;
        }
        return {mincost, maxflow};
    }
};

namespace solution {
    using i64 = long long;
    const int N = 1005;
    int n, k, ms, me, s[N], e[N], chk[N];
    MCMF<N, N << 1, int, i64> mcmf;

    void solve(){
        cin >> n >> k >> ms >> me;
        for(int i=1;i<=n;i++) cin >> s[i];
        for(int i=1;i<=n;i++) cin >> e[i];
        int S = 0, T = n + 1, total = k - ms;
        mcmf.add(S, 1, total, 0);
        for(int i=1;i<=n;i++) mcmf.add(i, i + 1, total - (i >= k ? me : 0), 0);
        for(int i=1;i<=n;i++) mcmf.add(i, std::min(i + k, T), 1, s[i] - e[i]), chk[i] = mcmf.ec;
        i64 cost = mcmf(S, T).first;
        cout << (std::accumulate(s + 1, s + n + 1, 0ll) - cost) << '\n';
        for(int i=1;i<=n;i++) cout << (mcmf.g[chk[i]].cap ? 'E' : 'S');
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_6967_NEERC_2016_Delight_for_a_Cat.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] ABC214H Collecting4000 ms

有一个 $n$ 个点 $m$ 条边的有向图。每个点有 $x_i$ 个遗物。现在有 $k$ 个人分别从 $1$ 出发，沿着任意路径走，收集经过的所有遗物。询问可以收集到多少个遗物。

$1\leq n,m\leq 10^5,1\leq k\leq 10$。

对于每一个 SCC，如果我们到达了其中的某个点，那么我们一定可以到达该 SCC 中的所有点，因此可以缩点，将原图缩成一个 DAG。然后考虑网络流建模，将每个点拆成 $u_i,v_i$，连边 $(u_i,v_i,1,x_i)$ 表示有 $x_i$ 个遗物在点 $i$ 上，只能收集一次。连边 $(u_i,v_i,\infty,0)$，表示可以经过该点任意多次，但是没有遗物。然后对于原图的边，连边 $(v_i,u_j,\infty,0)$，并且连边 $(S,u_1,\infty,0),(v_i,T,\infty,0)$，跑最大费用最大流即可。

一般的费用流，时间复杂度 $O(nmk)$（本图建模，总流量就是 $k$）无法通过本题。改为 Primal-Dual 原始对偶算法 + Dijkstra 求解费用流。但是初始势能如果用 SPFA 来求时间复杂度就会达到 $O(nm)$ 而无法通过。不过考虑到原图是一个 DAG，因此势能也可以通过拓扑排序在 DAG 上 dp 求出来，这样求初始势能的时间复杂度就达到了 $O(n+m)$。因此总时间复杂度为 $O(km\log m)$ 可以通过本题。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int, class Cost = int>
struct MCMF {
    struct edge {
        int nxt, to;
        Cap cap;
        Cost cost;
    } g[M << 1];

    int head[N], ec, S, T, cur[N];
    Cost dis[N], maxcost, phi[N];
    bool vis[N];

    MCMF() : ec(1) {}

    void add(int u, int v, Cap cap, Cost cost){
        g[++ec] = {head[u], v, cap, cost}, head[u] = ec;
        g[++ec] = {head[v], u, 0, -cost}, head[v] = ec;
    }

    int deg[N];

    void toposort(){
        for(int i=0;i<N;i++) phi[i] = -1e18;
        for(int i=2;i<=ec;i+=2) deg[g[i].to]++;
        std::queue<int> q;
        for(int i=0;i<N;i++) if(!deg[i]) q.push(i);
        phi[S] = 0;
        while(!q.empty()){
            int u = q.front(); q.pop();
            for(int i=head[u];i;i=g[i].nxt){
                if(i & 1) continue;
                int v = g[i].to;
                if(!deg[v]) continue;
                if(!(--deg[v])){
                    q.push(v);
                    for(int j=head[v];j;j=g[j].nxt){
                        if(j & 1) phi[v] = std::max(phi[v], phi[g[j].to] - g[j].cost);
                    }
                }
            }
        }
    }

    Cost cost(int id){ return phi[g[id ^ 1].to] - phi[g[id].to] + g[id].cost; }

    bool dijkstra(){
        for(int i=0;i<N;i++) dis[i] = std::numeric_limits<Cost>::min();
        std::priority_queue<std::pair<Cost, int>> q;
        q.emplace(dis[S] = 0, S);
        while(!q.empty()){
            int u = q.top().second; q.pop();
            if(vis[u]) continue;
            vis[u] = 1;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(g[i].cap && dis[u] + cost(i) > dis[v]){
                    dis[v] = dis[u] + cost(i), q.emplace(dis[v], v);
                }
            }
        }
        bool ans = vis[T];
        std::fill(vis, vis + N, 0);
        return ans;
    }

    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        vis[u] = 1;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(!vis[v] && g[i].cap && dis[u] + cost(i) == dis[v]){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t){
                    maxcost += t * g[i].cost;
                    g[i].cap -= t, g[i ^ 1].cap += t;
                    return vis[u] = 0, t;
                }
            }
        }
        return vis[u] = 0, 0;
    }

    std::pair<Cost, Cap> operator()(int s, int t){
        S = s, T = t, maxcost = 0;
        Cap maxflow = 0, tmp;
        toposort();
        while(dijkstra()){
            std::memcpy(cur, head, sizeof(head));
            while((tmp = dfs(s, std::numeric_limits<Cap>::max()))) maxflow += tmp;
            for(int i=0;i<N;i++){
                if(dis[i] != std::numeric_limits<Cost>::min()) phi[i] += dis[i];
            }
        }
        return {maxcost, maxflow};
    }
};

namespace solution{
    const int N = 2e5 + 5;
    int n, m, k, a[N];
    long long f[N];
    std::vector<int> g[N];
    MCMF<N * 2, N * 4, int, long long> mcmf;

    int dfn[N], low[N], stk[N], top, scc[N];
    bool vis[N];

    void tarjan(int u){
        dfn[u] = low[u] = ++dfn[0], stk[++top] = u, vis[u] = 1;
        for(int v : g[u]){
            if(!dfn[v]) tarjan(v), low[u] = std::min(low[u], low[v]);
            else if(vis[v]) low[u] = std::min(low[u], dfn[v]);
        }
        if(dfn[u] == low[u]){
            scc[u] = ++scc[0], vis[u] = 0;
            while(stk[top] != u) vis[stk[top]] = 0, scc[stk[top--]] = scc[u];
            top--;
        }
    }

    void solve(){
        cin >> n >> m >> k;
        for(int i=1,u,v;i<=m;i++) cin >> u >> v, g[u].push_back(v);
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=n;i++) if(!dfn[i]) tarjan(i);
        for(int i=1;i<=n;i++) f[scc[i]] += a[i];
        for(int i=1;i<=scc[0];i++){
            mcmf.add(i, i + scc[0], 1, f[i]), mcmf.add(i, i + scc[0], INT_MAX, 0);
            mcmf.add(i + scc[0], scc[0] * 2 + 1, INT_MAX, 0);
        }
        for(int u=1;u<=n;u++){
            for(int v : g[u]){
                if(scc[u] == scc[v]) continue;
                mcmf.add(scc[u] + scc[0], scc[v], INT_MAX, 0);
            }
        }
        mcmf.add(0, scc[1], k, 0);
        cout << mcmf(0, scc[0] * 2 + 1).first << '\n';
    }

    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_abc_214_h_ABC_214_H_Collecting.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium-Hard] P4480 [BJWC2018] 餐巾计划问题

一个餐厅在相继的 $n$ 天里，每天需用的餐巾数不尽相同。假设第 $i$ 天 $(i=1, 2, ..., n)$ 需要 $r_i$ 块餐巾。餐厅可以在任意时刻购买新的餐巾，每块餐巾的费用为 $p$ 。使用过的旧餐巾，则需要经过清洗才能重新使用。把一块旧餐巾送到清洗店A，需要等待 $m_1$ 天后才能拿到新餐巾，其费用为 $c_1$ ；把一块旧餐巾送到清洗店B，需要等待 $m_2$ 天后才能拿到新餐巾，其费用为 $c_2$ 。例如，将一块第 $k$ 天使用过的餐巾送到清洗店A清洗，则可以在第 $k+m_1$ 天使用。

请为餐厅合理地安排好 $n$ 天中餐巾使用计划，使总的花费最小。$1\leq n\leq 2\times 10^5$。

这道题有一个非常简单的网络流做法，但是如果拘泥于网络流做法的话是没有优化的空间的，考虑一般的策略。

如果我们钦定了购买了多少块餐巾，那么我们有一个显然的贪心做法：优先使用购买来的餐巾，如果购买来的餐巾不够用，那么再使用旧餐巾。旧餐巾中，我们总是优先使用最早的那块，并且能慢洗就慢洗，如果慢洗不够用，那么再快洗。这个过程的时间复杂度是 $O(n)$ 的。

考虑如何钦定购买多少块餐巾，直接枚举肯定是不行的。不过由于这道题有一个费用流做法，购买餐巾的数量是流量，那么由于费用流的流函数是凸的，因此可以用三分法找到最优的餐巾购买数量，时间复杂度 $O(n\log V)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution {
    const int N = 2e5 + 5;
    int n, m1, m2, c1, c2, p, r[N];

    int query(int init){
        int ans = init * p;
        std::deque<std::pair<int, int>> qWashing, qFastWashed, qSlowWashed;
        for(int i=1;i<=n;i++){
            // 快洗洗完了的，丢到对应的队列里
            while(!qWashing.empty() && qWashing.front().first <= i - m1){
                qFastWashed.push_back(qWashing.front());
                qWashing.pop_front();
            }
            // 慢洗洗完的，可以从快洗中回退
            while(!qFastWashed.empty() && qFastWashed.front().first <= i - m2){
                qSlowWashed.push_back(qFastWashed.front());
                qFastWashed.pop_front();
            }
            // 如果还有新买的毛巾，就直接用
            int tmp = r[i];
            int use = std::min(tmp, init);
            tmp -= use, init -= use;
            // 如果还需要毛巾，就先从慢洗中取，因为它更加便宜
            while(tmp && !qSlowWashed.empty()){
                auto& [time, cnt] = qSlowWashed.front();
                int use = std::min(tmp, cnt);
                tmp -= use, cnt -= use, ans += use * c2;
                if(!cnt) qSlowWashed.pop_front();
            }
            // 如果依然需要毛巾，就从快洗中取，注意先取新的，因为旧的有机会回退
            while(tmp && !qFastWashed.empty()){
                auto& [time, cnt] = qFastWashed.back();
                int use = std::min(tmp, cnt);
                tmp -= use, cnt -= use, ans += use * c1;
                if(!cnt) qFastWashed.pop_back();
            }
            // 检查一下还需要毛巾，如果需要那么初始购买的毛巾 init 不够
            if(tmp) return INT_MAX;
            // 今天的毛巾得洗了
            qWashing.push_back({i, r[i]});
        }
        return ans;
    }

    void solve(){
        cin >> n >> m1 >> m2 >> c1 >> c2 >> p;
        if(m1 > m2) std::swap(m1, m2), std::swap(c1, c2); // 保证慢洗耗时大于快洗
        if(c1 < c2) m2 = m1, c2 = c1; // 如果快洗更便宜，那么就全部用快洗
        for(int i=1;i<=n;i++) cin >> r[i];
        int L = 0, R = std::accumulate(r + 1, r + n + 1, 0);
        while(L < R){
            int mid = (L + R) >> 1;
            if(query(mid) < query(mid + 1)) R = mid;
            else L = mid + 1;
        }
        cout << query(L) << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_4480_BJWC_2018_餐巾计划问题.cpp] Think twice, code once.
// Written by xiezheyuan

[Easy-Medium] P3317 [SDOI2014] 重建

有 $n$ 个点，点 $(i,j)$ 间有无向边的概率为 $G_{i,j}$。求这 $n$ 个点恰好形成一棵树的概率。

$1<n\leq 50$。保证答案非零时，答案不小于 $10^{-4}$。

考虑答案即：

$$\begin{aligned} \sum_{T} \prod_{e\in T} G_e \prod_{e\not\in T} (1-G_e)=\prod_{e} (1-G_e) \sum_{T} \prod_{e\in T} \frac{G_e}{1-G_e} \end{aligned}$$

这是可以用矩阵树定理处理的形式，计算其 Krichhoff 矩阵的行列式即可，不过需要注意 $G_e$ 可能为 $0$ 的情况需要精细处理。时间复杂度 $O(n^3)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution {
    const int N = 55;
    int n;
    double G[N][N], deg[N];

    double det(){
        double res = 1;
        for(int k=1;k<n;k++){
            int pivot = k;
            for(int i=k+1;i<n;i++){
                if(std::fabs(G[i][k]) > std::fabs(G[pivot][k])) pivot = i;
            }
            if(pivot != k){
                for(int j=1;j<n;j++) std::swap(G[k][j], G[pivot][j]);
                res = -res;
            }
            if(std::fabs(G[k][k]) < 1e-10) return 0;
            for(int i=k+1;i<n;i++){
                double f = G[i][k] / G[k][k];
                for(int j=k;j<n;j++) G[i][j] -= f * G[k][j];
            }
            res *= G[k][k];
        }
        return res;
    }

    void solve(){
        cin >> n;
        double _ = 1;
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                cin >> G[i][j];
                if(std::fabs(G[i][j] - 1) < 1e-8) G[i][j] = 1 - 1e-8;
                if(i < j) _ *= (1 - G[i][j]);
                G[i][j] /= (1 - G[i][j]);
                if(i < j) deg[i] += G[i][j], deg[j] += G[i][j];
                G[i][j] = -G[i][j];
            }
        }
        for(int i=1;i<=n;i++) G[i][i] = deg[i];
        cout << std::fixed << std::setprecision(5) << _ * det() << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_3317_SDOI_2014_重建.cpp] Think twice, code once.
// Written by xiezheyuan

Have a Break Round II

太多太多构思踫壁 太多太多理想压抑 沉默对应着岁月无声

天空尚有片夕阳带领

踏上这无尽旅途 过去飘散消散失散花火 重燃起 重燃点起鼓舞

或许到最后没有完美句号 仍然倔强冒险一一去征讨

踏上这无尽旅途 谁又能鉴定你的丑恶与美好 来提步 壮阔跑道

合十双手去祷告 人生梦一场革命至苍老 难得梦一场革命不老

—— 《无尽》Supper Moment 世界变了样

[Medium] CF863F Almost Permutation3000 ms

给定一个长度为 $n$ 的未知序列 $a$，值域为 $[1,n]$。有如下 $q$ 条限制，形如：

• 区间 $[l,r]$ 的每个数均 $\leq v$。
• 区间 $[l,r]$ 的每个数均 $\geq v$。

记 $\mathrm{cnt}(i)$ 表示 $i$ 在 $a$ 中的出现次数。求 $\sum_{i=1}^n \mathrm{cnt}(i)^2$ 的最小值。若无解，输出 $-1$。

$1\leq n\leq 50,0\leq q\leq 100$。

通过这 $q$ 条限制，我们可以求出来每个数所在的一个区间 $[L_i,R_i]$。若 $[L_i,R_i]=\varnothing$ 则无解。

不妨钦定由位置去选择值，那么对于每一个 $a_i$ 建立一个点 $u_i$，对于每一个数 $x$ 建立一个点 $v_i$。对于 $x\in [L_i,R_i]$，连边 $(u_i,v_x,1,0)$ 表示 $a_i=x$。连边 $(S,u_i,1,0)$ 表示每个点只能选一个数。

那么我们现在需要连边 $(v_x,T)$，使得流过 $t$ 单位的流量会获得 $x^2$ 的代价，这看起来比较困难。

这里有一个被称为“凸费用”的模型，具体来说 $x^2$ 就是一个凸函数，我们取其差分 $\Delta(x)=x^2-(x-1)^2=2x-1$，那么 $2x-1$ 是递增的。我们可以连一组边 $(v_x,T,1,2i-1)$。这样第一个流量就会获得 $1$ 的代价，第 $i$ 个流量就会获得 $2i-1$ 的代价，这样从整体来看就是 $t^2$ 的代价。时间复杂度 $O(\mathrm{mcmf}(n,n^2))$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int, class Cost = int>
struct MCMF {
    struct edge {
        int nxt, to;
        Cap cap;
        Cost cost;
    } g[M << 1];

    int head[N], ec, S, T, cur[N];
    Cost dis[N], mincost;
    bool vis[N];

    MCMF() : ec(1) {}

    void add(int u, int v, Cap cap, Cost cost){
        g[++ec] = {head[u], v, cap, cost}, head[u] = ec;
        g[++ec] = {head[v], u, 0, -cost}, head[v] = ec;
    }

    bool spfa(){
        for(int i=0;i<N;i++) dis[i] = std::numeric_limits<Cost>::max();
        Cost inf = std::numeric_limits<Cost>::max();
        std::queue<int> q; q.push(S);
        dis[S] = 0;
        while(!q.empty()){
            int u = q.front();
            q.pop(), vis[u] = 0;
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(g[i].cap && dis[u] + g[i].cost < dis[v]){
                    dis[v] = dis[u] + g[i].cost;
                    if(!vis[v]) q.push(v), vis[v] = 1;
                }
            }
        }
        return dis[T] < inf;
    }

    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        vis[u] = 1;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(!vis[v] && g[i].cap && dis[u] + g[i].cost == dis[v]){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t){
                    mincost += t * g[i].cost;
                    g[i].cap -= t, g[i ^ 1].cap += t;
                    return vis[u] = 0, t;
                }
            }
        }
        return vis[u] = 0, 0;
    }

    std::pair<Cost, Cap> operator()(int s, int t){
        S = s, T = t, mincost = 0;
        Cap maxflow = 0, tmp;
        while(spfa()){
            std::memcpy(cur, head, sizeof(head));
            while((tmp = dfs(s, std::numeric_limits<Cap>::max()))) maxflow += tmp;
        }
        return {mincost, maxflow};
    }
};

namespace solution{
    const int N = 55;
    MCMF<N << 1, (N * N) << 1> mcmf;
    int n, q, L[N], R[N];

    void solve(){
        cin >> n >> q;
        std::fill(R + 1, R + n + 1, n), std::fill(L + 1, L + n + 1, 1);
        while(q--){
            int op, l, r, x; cin >> op >> l >> r >> x;
            for(int i=l;i<=r;i++){
                if(op == 1) L[i] = std::max(L[i], x);
                else R[i] = std::min(R[i], x);
            }
        }
        int S = 0, T = n + n + 1;
        for(int i=1;i<=n;i++){
            if(R[i] < L[i]) return cout << "-1\n", void();
            mcmf.add(S, i, 1, 0);
            for(int j=L[i];j<=R[i];j++) mcmf.add(i, j + n, 1, 0);
        }
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++) mcmf.add(i + n, T, 1, (j << 1) - 1);
        }
        cout << mcmf(S, T).first << '\n';
    }

    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [Almost_Permutation.cpp] Think twice, code once.
// Written by xiezheyuan

[Hard] P5470 [NOI2019] 序列

你需要确定两个长度为 $K$ 的序列 $\{c_i\}, \{d_i\}$，其中
$1 \leq c_1 < c_2 < \cdots < c_K \leq n , 1 \leq d_1 < d_2 < \cdots < d_K \leq n$

并要求 $\left \lvert \{c_1, c_2, \cdots , c_K\} \cap \{d_1, d_2, · · · , d_K\} \right \rvert \geq L$

目标是最大化 $\sum^{K}_{i=1} a_{c_i} +\sum^{K}_{i=1} b_{d_i}$。

$T$ 组数据。$1\leq T\leq 10,1\leq n\leq 2\times 10^5$。

这道题相当有难度。

考虑费用流建模。至少 $L$ 个相同的元素，我们可以看成至多 $K-L$ 个不同的元素，给这方面打上流量限制。而递增可以看成每个数只能选择一次。

具体地，我们钦定由 $c_i$ 选择 $d_i$，并且尽量选择相同的配对。连边 $(S,u_i,1,a_i),(v_i,T,1,b_i)$ 表示每个数只能选一次，且会带来 $a_i$ 或 $b_i$ 的代价。连边 $(S,u_i,v_i,1,0)$ 表示选择相同的数没有任何限制。对于选择不同的数，我们需要限制总流量。因此自然想到建立两个虚点 $U,V$，连边 $(u_i,U,1,0),(V,v_i,1,0)$。然后连边 $(U,V,K-L,0)$ 表示至多 $K-L$ 个不同的数。那么这样求一个最大费用最大流就是答案。建出来的图大致长这样：

时间复杂度爆炸。我们尝试模拟费用流。可以发现我们其实只有下面几种增广路：

Case 1：$S\to u_x\to v_x\to T$，即选择一对相同的数 $u_x,v_x$，费用为 $a_x+b_x$。

Case 2：$S\to u_x\to U\to V\to v_y\to T$，即选择一对不同的数 $u_x,v_y$，费用为 $a_x+b_y$。注意该操作会消耗边 $(U,V)$ 的 $1$ 个单位的容量（下面称为自由流量）：

Case 3：$S\to u_x\to v_x\to V\to v_y\to T$，即将原先与 $v_x$ 配对的改为与 $v_y$ 配对，然后将 $u_x$ 与 $v_x$ 配对，费用为 $a_x+b_y$：

Case 4：$S\to u_x\to U\to u_y\to v_y\to T$，即将原本与 $v_y$ 配对的改为与 $u_x$ 配对，然后将 $u_y$ 与 $v_y$ 配对，费用为 $a_x+b_y$：

Case 5：最后也是最难想到的一种，$S\to u_x\to v_x\to V\to U\to u_y\to v_y\to T$。这种选择不好解释，如果硬要解释的话可以理解成将原本与 $v_x$ 配对的与原本与 $u_y$ 配对的配对，然后让 $u_x,v_x$ 配对，让 $u_y,v_y$ 配对，费用为 $a_x+b_y$。这种情况会增加 $1$ 个单位的自由流量。

可以证明只有这五种情况（其他情况不优），然后自然考虑用若干个堆来维护增广。可以用五个堆 $q,q_u,q_v,q_u',q_v'$ 分别表示 $a_x+b_x$，$u_x,v_x$ 均未流过的 $a_x$ 和 $b_y$，以及只有另一侧流过的 $a_x$ 和 $b_y$。然后每次选出这些堆的堆顶，分类讨论哪一种增广最优即可。注意需要额外维护自由流量。

容易证明最多增广 $K$ 次。时间复杂度 $O((n+K)\log n)$，带有大常数。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution{
    using i64 = long long;
    const int N = 2e5 + 5;
    int n, L, K, a[N], b[N];
    bool visA[N], visB[N];
    std::priority_queue<std::pair<int, int>, std::vector<std::pair<int,int>>> qAwithB, qA, qB, qAsingle, qBsingle;

    template<class T> int topid(const T& pq){ return pq.top().second; }
    template<class T> int topv(const T& pq){ return pq.size() ? pq.top().first : -1e9; }
    void useA(int x){ visA[x] = 1, (!visB[x] ? qBsingle.emplace(b[x], x) : void()); }
    void useB(int x){ visB[x] = 1, (!visA[x] ? qAsingle.emplace(a[x], x) : void()); }

    void solve(){
        cin >> n >> K >> L, L = K - L;
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=n;i++) cin >> b[i];
        for(int i=1;i<=n;i++) qAwithB.emplace(a[i] + b[i], i), qA.emplace(a[i], i), qB.emplace(b[i], i);
        i64 ans = 0;
        for(int i=1;i<=K;i++){
            while(qAwithB.size() && (visA[topid(qAwithB)] || visB[topid(qAwithB)])) qAwithB.pop();
            while(qA.size() && visA[topid(qA)]) qA.pop();
            while(qB.size() && visB[topid(qB)]) qB.pop();
            while(qAsingle.size() && visA[topid(qAsingle)]) qAsingle.pop();
            while(qBsingle.size() && visB[topid(qBsingle)]) qBsingle.pop();
            int ret = INT_MIN, op = 0, AwithB = topv(qAwithB), A = topv(qA), B = topv(qB), Asingle = topv(qAsingle), Bsingle = topv(qBsingle);
            auto update = [&](int x, int id){ ret < x ? (ret = x, op = id) : 0; };
            update(AwithB, 1), update(Asingle + Bsingle, 3);
            update(Asingle + B, 4), update(A + Bsingle, 5);
            (L ? update(A + B, 2) : void()), ans += ret;
            if(op == 1) useA(topid(qAwithB)), useB(topid(qAwithB));
            else if(op == 2) useA(topid(qA)), useB(topid(qB)), L--;
            else if(op == 3) useA(topid(qAsingle)), useB(topid(qBsingle)), L++;
            else if(op == 4) useA(topid(qAsingle)), useB(topid(qB));
            else useA(topid(qA)), useB(topid(qBsingle));
        }
        cout << ans << '\n';
    }

    void clear(){
        auto clrq = [](auto& pq){ while(pq.size()) pq.pop(); };
        clrq(qAwithB), clrq(qA), clrq(qB), clrq(qAsingle), clrq(qBsingle);
        std::fill(visA, visA + n + 1, 0), std::fill(visB, visB + n + 1, 0);
    }
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T; cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_5470_NOI_2019_序列.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] P3153 [CQOI2009] 跳舞

有 $n$ 个男孩和 $n$ 个女孩，男孩和女孩之间有双向的喜欢关系。现在开始跳舞，每次播放舞曲时所有男女结成 $n$ 对跳舞。每一对男女最多跳一支舞，每个男孩最多和 $k$ 个不喜欢的女孩跳舞，每个女孩最多和 $k$ 个不喜欢的男孩跳舞。求总共最多能播放多少次舞曲。

$1\leq n\leq 50,0\leq k\leq 30$。

这是一个完美匹配的形式，直接做比较困难，考虑二分转成判定性问题。那么可以看成每个人都最多跳 $x$ 支舞，且总共结成了 $nx$ 对跳舞。那么可以网络流建模了，将每个男孩拆成三个点 $u_{i1},u_{i2},u_{i3}$，连边 $(u_{i1},u_{i2},\infty),(u_{i1},u_{i3},k)$，女孩也是类似的。然后对于相互喜欢的一对，连边 $(u_{i2},v_{j2},1)$，否则连边 $(u_{i3},v_{j3},1)$。最后从源点向每个 $u_{i1}$ 连边，从每个 $v_{j1}$ 向汇点连边。那么如果最大流为 $nx$，则是合法的。时间复杂度 $O(\mathrm{maxflow}(n,n^2)\log n)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int>
struct MaxFlow {
    struct node{
        int nxt, to;
        Cap cap;
    } g[M << 1];
    int head[N], ec, dep[N], cur[N], S, T;

    MaxFlow() : ec(1) {}

    void add(int u, int v, Cap c){
        g[++ec] = {head[u], v, c}, head[u] = ec;
        g[++ec] = {head[v], u, 0}, head[v] = ec;
    }

    bool bfs(){
        std::memset(dep, 0, sizeof(dep));
        std::queue<int> q;
        q.push(S), dep[S] = 1;
        while(!q.empty()){
            int u = q.front();
            q.pop();
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(!dep[v] && g[i].cap) q.push(v), dep[v] = dep[u] + 1;
            }
        }
        return dep[T];
    }
    
    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(dep[v] == dep[u] + 1 && g[i].cap){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t) return g[i].cap -= t, g[i ^ 1].cap += t, t;
            }
        }
        return 0;
    }

    Cap operator()(int s, int t){
        S = s, T = t;
        Cap mxflow = 0;
        while(bfs()){
            std::memcpy(cur, head, sizeof(cur));
            Cap x = 0;
            while((x = dfs(S, std::numeric_limits<Cap>::max()))) mxflow += x;
        }
        return mxflow;
    }
};

namespace solution {
    const int N = 55;
    int n, k, S, T;
    std::string s[N]; std::vector<int> vct;
    MaxFlow<(N << 1) * 3, N * N + N * 6> mf;
    int node(int x, int y){ return (x - 1) * 3 + y; }
    bool check(int x){
        for(int i=2;i<=mf.ec;i+=2) mf.g[i].cap += mf.g[i ^ 1].cap, mf.g[i ^ 1].cap = 0;
        for(int i : vct) mf.g[i].cap = x;
        return mf(S, T) == n * x;
    }
    void solve(){
        cin >> n >> k;
        for(int i=1;i<=n;i++) cin >> s[i], s[i] = ' ' + s[i];
        S = 0, T = node(n << 1, 3) + 1;
        for(int i=1;i<=n;i++){
            mf.add(S, node(i, 1), 1), vct.push_back(mf.ec ^ 1);
            mf.add(node(i, 1), node(i, 2), INT_MAX);
            mf.add(node(i, 2), node(i, 3), k);
            mf.add(node(i + n, 1), T, 1), vct.push_back(mf.ec ^ 1);
            mf.add(node(i + n, 2), node(i + n, 1), INT_MAX);
            mf.add(node(i + n, 3), node(i + n, 1), k);
            for(int j=1;j<=n;j++){
                if(s[i][j] == 'Y') mf.add(node(i, 2), node(j + n, 2), 1);
                else mf.add(node(i, 3), node(j + n, 3), 1);
            }
        }
        int L = 0, R = n;
        while(L < R){
            int mid = (L + R + 1) >> 1;
            if(check(mid)) L = mid;
            else R = mid - 1;
        }
        cout << L << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_3153_CQOI_2009_跳舞.cpp] Think twice, code once.
// Written by xiezheyuan

[Easy-Medium] ABC274G Security Camera 31024 MB

有一个高为 $H$、宽为 $W$ 的网格。第 $i$ 行第 $j$ 列的格子记作格子 $(i,j)$。
当 $S_{i,j}=$ # 时，格子上有障碍物；当 $S_{i,j}=$ . 时，格子上没有任何东西。

高桥君打算在网格中安装若干台监视摄像头。监视摄像头只能放在没有障碍物的格子上，并且可以朝上下左右四个方向中的任意一个方向放置。同一个格子上可以放置多台监视摄像头。

每台摄像头可以监视其所在的格子，以及其朝向方向上、在没有被障碍物阻挡的直线上所有格子。

要想监视所有没有障碍物的格子，最少需要安装多少台监视摄像头？
$1\leq H,W\leq 300$。

可以发现我们的摄像头覆盖的其实是以这个点出发，某个方向的最长连续段。那么容易注意到每次我们只会选择连续段的边界放置摄像头。由于在左侧和右侧，上方和下方放置摄像头是等价的，所以只需要纵向和横向分别只需要考虑一个方向。

每个空地在两个方向的连续段中，选择一个连续段放摄像头可以覆盖其中的所有空地，这个模型很像二分图最小点覆盖。我们只需要对横向连续段建立左部点，对纵向连续段建立右部点，对于每个空地，向其所在的横向和纵向连续段连边，那么最小点覆盖就是答案。根据 König 定理，最小点覆盖等于最大匹配，跑二分图最大匹配即可。时间复杂度 $O(HW\sqrt{HW})$。

实现

#include <bits/stdc++.h>
#include <atcoder/maxflow>
using std::cin, std::cout;

namespace solution {
    const int N = 305;
    int n, m, horiz[N][N], vert[N][N]; bool type[N * N];
    std::string s[N];

    void solve(){
        cin >> n >> m;
        for(int i=1;i<=n;i++) cin >> s[i], s[i] = ' ' + s[i];
        int t = 0;
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                if(s[i][j] == '#') continue;
                if(!vert[i][j]){
                    type[++t] = 1; int c = i;
                    while(c <= n && s[c][j] == '.') vert[c++][j] = t;
                }
                if(!horiz[i][j]){
                    t++; int c = j;
                    while(c <= m && s[i][c] == '.') horiz[i][c++] = t;
                }
            }
        }
        atcoder::mf_graph<int> g(t + 2);
        int S = 0, T = t + 1;
        for(int i=1;i<=t;i++) type[i] ? g.add_edge(S, i, 1) : g.add_edge(i, T, 1);
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                if(s[i][j] == '.') g.add_edge(vert[i][j], horiz[i][j], 1);
            }
        }
        cout << g.flow(S, T) << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_abc_274_g_ABC_274_G_Security_Camera_3.cpp] Think twice, code once.
// Written by xiezheyuan

Have a Break Round III

刻版的传记会给淡忘 想小说悲壮 人物那心理曝了光 弦外那声音会解放

从感情出发 现实梦境穿插 不必百分之百 如实去记低那份平白

人生如画册 我有我的浪漫手法 剪贴着记忆 将最完美的定格

难得有颜色笔填色册 来让我努力乱填也别留白

从感情出发 现实梦境穿插 想怎拍便怎拍 桥段我觉得对便成立

人生如画册 我有我的浪漫手法 剪贴着记忆 将最完美的定格

—— 《心照》郑伊健 Friends For Life

[Easy] P2172 [国家集训队] 部落战争

给定一个 $M \times N$ 的网格，其中包含城镇（'.'）和高山深涧（'x'）。需要从城镇出发，向下行进，每次移动可以跳跃 $R$ 行和 $C$ 列，即从 $(r, c)$ 可以移动到 $(r+R, c+C),(r+C, c-R),(r+R,c-C),(r+C,c+R)$。军队只能经过城镇，且一个城镇只能被一支军队经过。目标是求最少需要多少支军队才能覆盖所有城镇。

$1\leq M,N\leq 50$。

显然这是一个 DAG 最小不交路径覆盖问题，这里讲讲如何解决这类问题以防忘记。只需要将每个点拆成一个入点一个出点，然后正常建图，那么你会得到一个二分图，最后的点数 - 匹配数就是答案。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

template<int N, int M, class Cap = int>
struct MaxFlow {
    struct node{
        int nxt, to;
        Cap cap;
    } g[M << 1];
    int head[N], ec, dep[N], cur[N], S, T;

    MaxFlow() : ec(1) {}

    void add(int u, int v, Cap c){
        g[++ec] = {head[u], v, c}, head[u] = ec;
        g[++ec] = {head[v], u, 0}, head[v] = ec;
    }

    bool bfs(){
        std::memset(dep, 0, sizeof(dep));
        std::queue<int> q;
        q.push(S), dep[S] = 1;
        while(!q.empty()){
            int u = q.front();
            q.pop();
            for(int i=head[u];i;i=g[i].nxt){
                int v = g[i].to;
                if(!dep[v] && g[i].cap) q.push(v), dep[v] = dep[u] + 1;
            }
        }
        return dep[T];
    }
    
    Cap dfs(int u, Cap flow){
        if(u == T) return flow;
        for(int &i=cur[u];i;i=g[i].nxt){
            int v = g[i].to;
            if(dep[v] == dep[u] + 1 && g[i].cap){
                Cap t = dfs(v, std::min(flow, g[i].cap));
                if(t) return g[i].cap -= t, g[i ^ 1].cap += t, t;
            }
        }
        return 0;
    }

    Cap operator()(int s, int t){
        S = s, T = t;
        Cap mxflow = 0;
        while(bfs()){
            std::memcpy(cur, head, sizeof(cur));
            Cap x = 0;
            while((x = dfs(S, std::numeric_limits<Cap>::max()))) mxflow += x;
        }
        return mxflow;
    }
};

namespace solution{
    const int N = 55;
    int n, m, R, C;
    MaxFlow<N * N * 2, N * N * 6> G;
    std::string s[N];
    int node(int x, int y, int op){ return ((x - 1) * m + (y - 1)) * 2 + op; }
    void chkadd(int x1, int y1, int x2, int y2){
        if(x2 < 1 || x2 > n || y2 < 1 || y2 > m || s[x2][y2] == 'x') return;
        G.add(node(x1, y1, 1), node(x2, y2, 2), 1);
    }
    void solve(){
        cin >> n >> m >> R >> C;
        for(int i=1;i<=n;i++) cin >> s[i], s[i] = ' ' + s[i];
        int S = 0, T = node(n, m, 2) + 1, cnt = 0;
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                if(s[i][j] == 'x') continue; cnt++;
                chkadd(i, j, i + R, j + C), chkadd(i, j, i + C, j - R);
                chkadd(i, j, i + R, j - C), chkadd(i, j, i + C, j + R);
                G.add(S, node(i, j, 1), 1), G.add(node(i, j, 2), T, 1);
            }
        }
        cout << cnt - G(S, T) << '\n';
    }
    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_2172_国家集训队_部落战争.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] ABC320G Slot Strategy 2 (Hard)1024 MB

给定 $N$ 个字符串 $S_1, S_2, \dots, S_N$，每个字符串长度为 $M$，只包含数字。现在有 $N$ 个转轮，第 $i$ 个转轮对应字符串 $S_i$。在时刻 $t$ 按下第 $i$ 个转轮的按钮，它会停在 $S_i$ 的第 $(t \pmod M) + 1$ 个字符上。要求选择 $N$ 个不同的时刻 $t_1, t_2, \dots, t_N$，分别按下 $N$ 个转轮的按钮，使得所有转轮停止后显示的字符都相同。求满足条件的最早的时刻，即 $\max(t_1, t_2, \dots, t_N)$ 的最小值。如果无法实现，输出 $-1$。

$1\leq N\leq 100,1\leq M\leq 10^5$。

感觉这题真的不简单。我们先枚举最后显示的字符是什么，然后考虑用转盘去选择时刻，那么这是一个完美匹配型问题，二分这个时间最小值，如果在一个时刻按下按钮可以让这个转轮转到这个现实的字符，就连一条边，最后如果匹配数为 $N$ 则说明可行。

这个做法乍一看好像没有问题，不过有一个很严重的问题是时间复杂度为 $O((N+M)\sqrt{N+M}\log NM)$ 而 $M$ 很大。这里给出一个引理（这个引理看起来很玄学，但是细细想来确实是正确的）：

引理：对于一个二分图，如果存在左部点到右部点的完美匹配，那么我们选择一个度数大于 $n$ 的左部点并删除它的一条边，重复上述操作任意次，那么最终得到的二分图仍然存在左部点到右部点的完美匹配。

这个引理很容易用 Hall 定理证明，这里不再赘述。应用这个定理，我们便可以将右部点缩减到 $O(N^2)$。时间复杂度 $O(N^3\log NM)$，带有一个 $10$ 倍常数。

实现

#include <bits/stdc++.h>
#include <atcoder/maxflow>
using std::cin, std::cout;

namespace solution {
    const int N = 105;
    int n, m;
    std::string s[N];
    std::basic_string<int> g[10][N];
    bool check(int num, int tim){
        std::map<int, int> mp; int S = 0, T = 1, tot = n + 1;
        for(int i=1;i<=n;i++){
            int cnt = 0;
            for(int delta=0;delta<=tim;delta+=m){
                for(int j : g[num][i]){
                    int t = delta + j;
                    if(t > tim) break;
                    if(!mp.count(t)) mp[t] = ++tot;
                    cnt++;
                    if(cnt == n) break;
                }
                if(cnt == n) break;
            }
        }
        atcoder::mf_graph<int> G(tot + 1);
        for(int i=1;i<=n;i++) G.add_edge(S, i + 1, 1);
        for(int i=n+2;i<=tot;i++) G.add_edge(i, T, 1);
        for(int i=1;i<=n;i++){
            int cnt = 0;
            for(int delta=0;delta<=tim;delta+=m){
                for(int j : g[num][i]){
                    int t = delta + j;
                    if(t > tim) break;
                    G.add_edge(i + 1, mp[t], 1);
                    cnt++;
                    if(cnt == n) break;
                }
                if(cnt == n) break;
            }
        }
        return G.flow(S, T) == n;
    }
    void solve(){
        cin >> n >> m;
        for(int i=1;i<=n;i++){
            cin >> s[i];
            for(int j=0;j<m;j++) g[s[i][j] - '0'][i] += j;
        }
        int ans = INT_MAX;
        for(int i=0;i<9;i++){
            int L = 0, R = n * m + 1;
            while(L < R){
                int mid = (L + R) >> 1;
                if(check(i, mid)) R = mid;
                else L = mid + 1;
            }
            ans = std::min(ans, L);
        }
        cout << (ans > n * m ? -1 : ans) << "\n";
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [Slot_Strategy_2_Hard.cpp] Think twice, code once.
// Written by xiezheyuan

[Hard] AT_arc080_d [ARC080F] Prime Flip

有一个无限长的序列，一开始 $n$ 个位置是 $1$，其他位置是 $0$。每次操作可以选择一个不小于 $3$ 的素数 $p$，然后选择序列中的 $p$ 个连续位置进行翻转（$0$ 变 $1$，$1$ 变 $0$）。求把所有 $1$ 都变为 $0$ 的最少操作次数。

$1\leq n\leq 100$。最大为 $1$ 的位置 $\leq 10^7$。

先对整个序列异或差分，这样选择一个长度为 $p$ 的连续区间 flip，等价于选择一个 $i$，将 $a_i$ 和 $a_{i+p}$ flip。然后我们可以很方便地求出来在差分后的序列中，哪些位置是 $1$。然后我们不妨枚举两个为 $1$ 的位置 $i,j(i<j)$，考察最少进行多少次操作可以将这两个位置变成 $0$。

Case 1：$j-i\in\mathbb{P}\land j-i\neq 2$，此时只需要一次操作即可。

Case 2：$2\mid j-i$，此时：

• 根据 Goldbach 猜想，若 $j-i>2$，则可以写成两个素数的和。而这里我们要求奇素数，对 Goldbach 猜想的命题稍加转换可以得到任意大于 $4$ 的偶数都可以写成两个奇素数之和。假设 $j-i=p_1+p_2$，那么我们只需要对 $i,i+p_1$ 进行一次操作，$i+p_1,j$ 进行一次操作即可。
• 若 $j-i=2$，我们注意到 $2=5-3$，所以可以先对 $i,i+5$ 进行一次操作，再对 $j,i+5$ 进行一次操作即可。
• 若 $j-i=4$，我们注意到 $4=7-3$，所以可以先对 $i,i+7$ 进行一次操作，再对 $j,i+7$ 进行一次操作即可。

因此当 $2\mid j-i$ 时，我们只需要两次操作。

Case 3：$2\nmid j-i\land j-i\not\in\mathbb{P}$，由于偶数个奇数相加 / 相减，结果一定是偶数，所以对于该情况，我们至少需要三次操作，下面证明三次操作总是可行的。对于任意一个奇数 $x$，我们一定可以选择一个奇素数 $p$，使得 $x+p$ 为一个 $>4$ 偶数（这个结论看起来很平凡，实际上也很平凡）。那么根据 Goldbach 猜想，我们仍然可以 $x+p=p_1+p_2\implies x=p_1+p_2-p$。因此我们只需要对 $i,i+p_1$ 进行一次操作，对 $i+p_1,i+p_1+p_2$ 进行一次操作，最后对 $j,i+p_1+p_2$ 进行一次操作即可。因此，对于该情况，我们只需要三次操作。

综上所述，我们一定要尽量执行 Case 1 的操作，注意到若 $j-i\in\mathbb{P}\land j-i\neq 2$，则 $j\not\equiv i\pmod{2}$。那么我们可以对于 $j-i$ 满足该条件，且 $i,j$ 差分后均为 $1$ 的情况连一条边 $(i,j)$，那么这是一个二分图，求出最大匹配就是 Case 1 的最大进行次数。

现在我们让一些位置达到了 $0$，对于剩下的操作，我们要优先执行 Case 2。那么非常显然地，我们可以将剩下的为 $1$ 的位置按照奇偶性分成两类，每类之间进行 Cae 2 的操作。如何统计每一类有多少个数？只需要预先处理为 $1$ 的奇数点、偶数点数量，然后跑完最大匹配后，每个数量减去匹配数即可（因为匹配的每条边，端点奇偶性不同）。

剩下的我们就要执行 Case 3。不过在此之前，我们必须证明一个引理，保证我们以前的操作没有胡来：

初始时异或差分数组为 $1$ 的点数是偶数。

证明起来也很简单，考虑原序列的一个 $1$，它理想情况下会给异或差分后的数组贡献两个 $1$，如果只贡献了一个 $1$，那么说明这个 $1$ 的位置在异或差分数组上以前就是 $1$，现在他们俩抵消了，相当于减去了 $2$。所以总是偶数。

于是尽可能执行 Case 2 后，剩下的要么没有数了，要么偶数类和奇数类各剩下一个数，那么我们只需要对剩下的两个数执行 Case 3 的操作即可。最后将所有情况对应的操作数加起来就是答案。时间复杂度上界为 $O(\mathrm{maxflow}(n,n^2))$。

实现

#include <bits/stdc++.h>
#include <atcoder/math>
#include <atcoder/maxflow>
using std::cin, std::cout, atcoder::internal::is_prime_constexpr;

namespace solution {
    int n, t[2], tot; std::map<int, int> mp, idt;

    void solve(){
        cin >> n;
        for(int i=1,x;i<=n;i++) cin >> x, mp[x - 1] ^= 1, mp[x] ^= 1;
        for(auto [p, v] : mp) if(v) idt[p] = ++tot, t[p & 1]++;
        int S = 0, T = tot + 1; atcoder::mf_graph<int> g(tot + 2);
        for(auto [p, v] : mp) if(v) p & 1 ? g.add_edge(S, idt[p], 1) : g.add_edge(idt[p], T, 1);
        for(auto [p1, v1] : mp){
            if(!(p1 & 1) || !v1) continue;
            for(auto [p2, v2] : mp){
                if(!(p2 & 1) && v2 && is_prime_constexpr(abs(p2 - p1))) g.add_edge(idt[p1], idt[p2], 1);
            }
        }
        int flow = g.flow(S, T);
        flow += (((t[0] - flow) >> 1) + ((t[1] - flow) >> 1)) << 1;
        cout << flow + (t[1] - flow & 1) * 3 << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_arc_080_d_ARC_080_F_Prime_Flip.cpp] Think twice, code once.
// Written by xiezheyuan

[Hard] ARC156F Make Same Set4000 ms1024 MB

给定三个长度为 $N$ 的整数序列 $A, B, C$。

需要找到一个最大的集合，使得该集合可以由两种方式生成：

1. 从空集开始，对于 $i = 1, \dots, N$，在每一步中添加 $A_i$ 或 $B_i$。
2. 从空集开始，对于 $i = 1, \dots, N$，在每一步中添加 $A_i$ 或 $C_i$。

求出该最大集合的元素个数以及集合中的元素。$1\leq N\leq 5\times 10^3,1\leq A_i,B_i,C_i\leq 10^4$。

这道题真让人脑洞大开。我们先来考虑一个弱化版：

你需要找到一个最大的集合，使得该集合可以由两种方式生成：

1. 从空集开始，对于 $i = 1, \dots, N$，在每一步中添加 $A_i$ 或 $B_i$，或者什么也不做。
2. 从空集开始，对于 $i = 1, \dots, N$，在每一步中添加 $A_i$ 或 $C_i$，或者什么也不做。

对于该问题，我们有一个平凡的网络流做法。对于每一个 $i=1,\cdots,N$，我们建立两个点 $AB_i,AC_i$。然后连边 $(S,AB_i,1),(AC_i,T,1)$。然后对于每一个数，拆成两个点 $u_i,v_i$，然后连边 $(AB_i,u_{a_i},1),(AB_i,u_{b_i},1),(v_{a_i},AC_i,1),(v_{c_i},AC_i,1)$。最后连边 $(u_i,v_i,1)$ 表示集合中每个数只能出现一次。跑最大流即可。容易证明这个做法对于弱化版是正确的，且对于原问题，求出的是答案的一个上界。

考虑为什么说这样求出的上界，而不能直接说求出的就是答案呢？先摆出我们建出来的图：

如果对于一个 $AB_x,AC_y$，如果它连接的中间点没有流过任何流，那么我们就不知道它应该选什么。我们称这样的点为待定点。然后不妨分析一下在最大流增广过程中，待定点是如何产生的。我们可以证明这样的一个事实：

如果我们在一次增广中出现了一个待定点，那么一定存在一条更短的增广路径。

放一张官方题解的图：

在上图中，$AB_5$ 变成了待定点，这是因为我们撤回了 $u_3,v_3$ 的增广而改为增广 $u_4,v_4$。那么我们完全可以改为只撤销 $(AB_3,3)$ 的增广，这样就不会出现待定点。

因此如果我们可以建出一张没有待定点的图，那么采用一种最大流算法，该算法的过程是每次选择最长的一条增广路增广，那么求出的上界就是答案。这样的最大流算法其实就是 Dinic。建图的方法也很简单，强制钦定每个 $a_i$ 都必选即可，从这样的残差网络开始跑 Dinic 即可。

实现

#include <bits/stdc++.h>
#include <atcoder/maxflow>
using std::cin, std::cout;

namespace solution {
    const int N = 5e3 + 5, M = 1e4 + 5;
    int n, a[N], b[N], c[N];
    atcoder::mf_graph<int> g(N * 2 + M * 2);
    void solve(){
        cin >> n;
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=n;i++) cin >> b[i];
        for(int i=1;i<=n;i++) cin >> c[i];
        int S = 0, T = N * 2 + M * 2 - 1;
        for(int i=1;i<=n;i++){
            g.add_edge(S, i, 1);
            g.add_edge(i, (a[i] + (n << 1)), 1);
            g.add_edge(i + n, T, 1);
            g.add_edge(a[i] + (n << 1) + M, i + n, 1);
        }
        for(int i=1;i<=M;i++) g.add_edge(i + (n << 1), i + (n << 1) + M, 1);
        g.flow(S, T);
        for(int i=1;i<=n;i++){
            g.add_edge(i, b[i] + (n << 1), 1);
            g.add_edge(c[i] + (n << 1) + M, i + n, 1);
        }
        std::map<std::pair<int, int>, int> mp;
        g.flow(S, T);
        for(auto [u, v, cap, flow] : g.edges()) mp[{u, v}] = flow;
        std::basic_string<int> vct;
        for(int i=1;i<=M;i++){
            if(mp[{i + (n << 1), i + (n << 1) + M}]) vct += i;
        }
        cout << vct.size() << '\n';
        for(int x : vct) cout << x << ' ';
        cout << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_arc_156_f_ARC_156_F_Make_Same_Set.cpp] Think twice, code once.
// Written by xiezheyuan

Have a Break Round IV

狂风 猛风 怒似狂龙

巨浪乱撞乱动 雷电用疾劲暴力划破这世界 震裂天空

狂奔雨中 视野蒙眬 狂叫着 Lorelei Sweet Lorelei

神秘的 Lorelei 居于风眼中 活着只喜欢破坏 事后无踪

曾在那揭不开的漆黑中 压不低的恋火中 猛风当中她拥抱我 Oh 这刻身边只有风

停在那抹不清的追忆中 割不穿的空虚中 猛风当中只有我

跌跌碰碰 再次发觉她于风中

—— 《暴风女神 LORELEI》 谭咏麟 电影《龙兄虎弟》主题曲

[Medium] ARC125E Snack1024 MB

有 $n$ 种零食，第 $i$ 种零食有 $a_i$ 个；有 $m$ 个小朋友，第 $i$ 个小朋友每种零食最多吃 $b_i$ 个，总共最多吃 $c_i$ 个。求小朋友们总共最多能吃多少个零食。

$1\leq n,m\leq 2\times 10^5$。

首先本题有一个比较简单的网络流做法，对于每一个零食建一个点 $u_i$，连边 $(S,u_i,a_i)$，对于每一个小朋友建一个点 $v_i$，连边 $(v_i,T,c_i)$，然后连边 $(u_i,v_j,b_j)$。跑最大流即可，时间复杂度爆炸。

我们似乎要模拟最大流，但是这比较麻烦，一种更加好的方法是模拟最小割。具体地，我们利用最大流最小割定理将最大流转成最小割，然后尝试刻画这个图的最小割的形式。

由于小朋友的贡献由两个量控制，而零食的贡献只有一个量控制，这些量都是独立的。不妨枚举我们割断 $x$ 多少条 $(S,u_i,a_i)$ 的边，那么贪心地，我们总是选择 $a_i$ 前 $x$ 小的边割断。那么剩下的 $n-x$ 种零食都在割集 $c(S)$ 种。对于每一个小朋友，它要么选择割断 $(v_j,T,c_j)$，要么选择割断所有 $(u_i,v_j,b_j)$。前者的代价是 $c_j$，后者的代价是 $(n-x)b_j$。

因此我们可以将所有小朋友按照 $\frac{c_j}{b_j}$ 排序，然后双指针扫一遍就可以求出最优的 $x$ 和此时的最小割，取个最小值即可。时间复杂度 $O(n\log n+m\log m)$，瓶颈在排序。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution {
    using i64 = long long;
    const int N = 2e5 + 5;
    int n, m;
    i64 a[N], c[N];
    struct node { i64 b, c; } b[N];

    void solve(){
        cin >> n >> m;
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=m;i++) cin >> b[i].b;
        for(int i=1;i<=m;i++) cin >> b[i].c;
        std::sort(a + 1, a + n + 1);
        std::sort(b + 1, b + m + 1, [](node x, node y) { return x.c * y.b < y.c * x.b;; });
        for(int i=1;i<=m;i++) c[i] = c[i - 1] + b[i].c;
        i64 ans = LLONG_MAX, pre = 0, sum = 0; int R = m;
        for(int i=0;i<=n;i++){
            while(R && (n - i) * b[R].b < b[R].c) sum += b[R--].b;
            ans = std::min(ans, (pre += a[i]) + sum * (n - i) + c[R]);
        }
        cout << ans << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_arc_125_e_ARC_125_E_Snack.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium-Hard] ABC332G Not Too Many Balls3000 ms1024 MB

有若干个球。每个球的颜色为颜色 $1$、颜色 $2$、$\ldots$、颜色 $N$ 之一，对于 $i=1,2,\ldots,N$，颜色 $i$ 的球共有 $A_i$ 个。另外，有 $M$ 个箱子。对于 $j=1,2,\ldots,M$，第 $j$ 个箱子最多可以放 $B_j$ 个球。

但是，对于所有满足 $1 \leq i \leq N$ 且 $1 \leq j \leq M$ 的整数对 $(i, j)$，将颜色 $i$ 的球放入第 $j$ 个箱子的数量不能超过 $i \times j$ 个。

请输出可以放入 $M$ 个箱子中的球的最大总数。

$1\leq N\leq 500,1\leq M\leq 5\times 10^5$。

与上一题类似地，我们也考虑一个最大流做法。很容易想到对于每一个颜色的球 $i$，建立一个点 $u_i$，连边 $(S,u_i,A_i)$，然后对于每一个箱子 $j$，建立一个点 $v_j$，连边 $(v_j,T,B_j)$。然后对于每一个颜色 $i$ 的球，连边 $(u_i,v_j,i\times j)$。跑最大流即可，不过这样无法通过本题。

考虑用最大流最小割定理转成最小割问题，考虑最小割怎么做。不妨将 $c(S)$ 中的球集合设为 $P$，$c(T)$ 中的箱集合设为 $Q$，那么我们求的其实是：

$$\min_{P}\min_{Q} \left(\sum_{i\not\in P} A_i+\sum_{i\not\in Q} B_i+\sum_{i\in P}\sum_{j\in Q} ij\right)=\min_P\min_Q\left(\sum_{i\not\in P}A_i+\sum_{i\not\in Q}B_i+\sum_{i\in P}i\sum_{j\in Q}j\right)$$

由于 $\sum_{i\in P}i$ 的规模只有 $O(N^2)$ 级别，所以可以尝试枚举它：

$$\begin{aligned} &\min_{k=0}^{\frac{1}{2} n(n+1)}\min_{\sum_{i\in P}i=k}\min_{Q}\left(\sum_{i\not\in P} A_i+\sum_{i\not\in Q}B_i+k\sum_{j\in Q}j\right)=\min_{k=0}^{\frac{1}{2} n(n+1)}\min_{\sum_{i\in P}i=k}\left(\sum_{i\not\in P} A_i + \sum_{j=1}^{m} \min(B_j,kj)\right)\\ &=\min_{k=0}^{\frac{1}{2} n(n+1)}\left(\sum_{j=1}^{m} \min(B_j,kj) + \min_{\sum_{i\in P}i=k}\sum_{i\not\in P} A_i\right) \end{aligned}$$

枚举 $k$ 后，我们将解决的问题分成了两个独立的部分：

• $\sum_{j=1}^{m} \min(B_j,kj)$ 只需要将所有 $B_j$ 按照 $\frac{B_j}{j}$ 排序后双指针即可。
• $\min_{\sum_{i\in P}i=k}\sum_{i\not\in P} A_i$ 很容易用 01 背包状的 dp 解决。

时间复杂度 $O(n^3+m\log m)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution{
    using i64 = long long;
    const int N = 505, M = 5e5 + 5;
    int n, m; i64 a[N], f[N * N], pre[M];
    struct node { int id; i64 v; } b[M];

    void solve(){
        cin >> n >> m; int L = n * (n + 1) / 2;
        for(int i=1;i<=n;i++) cin >> a[i];
        for(int i=1;i<=m;i++) cin >> b[i].v, b[i].id = i;
        std::fill(f + 1, f + L + 1, 1e18);
        for(int i=1;i<=n;i++){
            for(int j=L;~j;j--) f[j] = std::min(f[j] + a[i], j >= i ? f[j - i] : LLONG_MAX);
        }
        std::sort(b + 1, b + m + 1, [](node x, node y){ return x.v * y.id < y.v * x.id; });
        for(int i=1;i<=m;i++) pre[i] = pre[i - 1] + b[i].v;
        i64 ans = LLONG_MAX, rhs = 1ll * m * (m + 1) / 2, l = 1;
        for(int i=0;i<=L;i++){
            while(l <= m + 1 && b[l].v < 1ll * b[l].id * i) rhs -= b[l++].id;
            ans = std::min(ans, f[i] + rhs * i + pre[l - 1]);
        }
        cout << ans << '\n';
    }

    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [G_Not_Too_Many_Balls.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] AT_arc085_c [ARC085E] MUL

有 $N$ 个物品，编号为 $1, 2, \dots, N$。对于每个物品 $i$，如果你保留它，你会得到 $a_i$ 的价值；如果你不保留它，则价值为 $0$。你可以执行任意次操作：选择一个正整数 $x$，然后将所有编号为 $x$ 的倍数的物品全部不保留。请你决定最终哪些物品保留，哪些不保留，使得总价值最大。

$1\leq N\leq 100,|a_i|\leq 10^9$。

最大权闭合子图题。考虑原题的描述转置一下，可以变成执行任意次操作，每次选择一个正整数 $x$，然后将所有编号为 $x$ 的因数的物品全部保留。这样可以考虑对于 $i\mid j$，连边 $(j,i)$，这样我们就变成了最大权闭合子图模型（选择一个点，必须选择它的出点）。

这个东西也很经典，不过值得记下来。考虑将点分成正权点和负权点。对于正权点，连边 $(S,i,a_i)$，对于负权点，连边 $(i,T,-a_i)$，对于原本的边 $(i,j)$，连边 $(i,j,+\infty)$。这样，最大权闭合子图就是正权点点权总和，减去最大流。

实现

#include <bits/stdc++.h>
#include <atcoder/maxflow>
using std::cin, std::cout;

namespace solution{
    using i64 = long long;
    const int N = 105;
    i64 a[N], ans; int n, S, T;
    void solve(){
        cin >> n; atcoder::mf_graph<i64> g(n + 2); S = 0, T = n + 1;
        for(int i=1;i<=n;i++){
            cin >> a[i];
            if(a[i] > 0) g.add_edge(S, i, a[i]), ans += a[i];
            else g.add_edge(i, T, -a[i]);
        }
        for(int i=1;i<=n;i++){
            for(int j=i;j<=n;j+=i) g.add_edge(j, i, LLONG_MAX);
        }
        cout << ans - g.flow(S, T) << '\n';
    }
    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [E_MUL.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] ABC259G Grid Card Game1024 MB

有一个 $H \times W$ 的网格，每个单元格 $(i, j)$ 有一个值 $A_{i, j}$。你可以选择一些行和一些列。如果一个单元格所在的行或列被选中，则该单元格被标记。

如果一个值为负数的单元格同时所在的行和列都被选中，则得分为一个非常小的负数。否则，得分为所有被标记的单元格的值的总和。求最大可能的得分。

$1\leq H,W\leq 100$。

首先答案有一个下界 $0$，所以我们不应该让选中的行和列的交集中有负数。考虑我们先将所有正数都选上，然后考虑扣分：

• 如果行 $i$ 和列 $j$ 都没有选中，且 $A_{i,j}>0$，那么我们扣去 $A_{i,j}$。
• 那么如果选取了行 $i$，那么扣除 $\sum_{j=1}^W -\min(A_{i,j},0)$ 分，选取列 $j$ 同理。

这个问题怎么做呢？你可以考虑对于每个行建一个点 $u_i$，每一列建一个点 $v_i$。连边 $(S,u_i,\sum_{j=1}^W -\min(A_{i,j},0))$ 和 $(v_i,T,\sum_{j=1}^{H} -\min(A_{j,i},0))$ 表示选取了行 $i$ 与列 $j$ 的代价。对于 $A_{i,j}>0$，连边 $(u_i,v_j,A_{i,j})$ 表示如果 $(S,u_i),(v_i,T)$ 都没有割断，那么 $A_{i,j}$ 无法贡献，被扣除。

考虑对于 $A_{i,j}<0$，我们要求 $(S,u_i),(v_j,T)$ 都不能割断，可以尝试连边 $(v_j,u_i,+\infty)$。

实现

#include <bits/stdc++.h>
#include <atcoder/maxflow>
using std::cin, std::cout;

namespace solution{
    using i64 = long long;
    const int N = 105;
    int n, m; i64 a[N][N], row[N], col[N];
    void solve(){
        cin >> n >> m; atcoder::mf_graph<i64> g(n + m + 2);
        int S = 0, T = n + m + 1; i64 ans = 0;
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                cin >> a[i][j];
                if(a[i][j] > 0) ans += a[i][j], g.add_edge(i, j + n, a[i][j]);
                else row[i] -= a[i][j], col[j] -= a[i][j], g.add_edge(j + n, i, LLONG_MAX);
            }
        }
        for(int i=1;i<=n;i++) g.add_edge(S, i, row[i]);
        for(int i=1;i<=m;i++) g.add_edge(i + n, T, col[i]);
        cout << ans - g.flow(S, T) << '\n';
    }
    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [AT_abc_259_g_ABC_259_G_Grid_Card_Game.cpp] Think twice, code once.
// Written by xiezheyuan

[Hard] P13275 [NOI2025] 集合

给定 $n$ 和一个长度为 $2^n$ 的序列 $a$，定义 $U=[0,2^n)\cap\mathbb{Z}$，计算：

$$\sum_{\substack{P\subseteq U,Q\subseteq U\\P\cap Q=\varnothing}} \left[\bigwedge_{k\in P} k=\bigwedge_{k\in Q}k\right]\prod_{k\in P\cup Q} a_k$$

答案对 $p=998,244,353$ 取模。

$T$ 组数据。$1\leq T\leq 3,2\leq n\leq 20,0\leq a_i<p$。

记号约定

• $\land,\lor$ 分别表示按位与、或运算。
• $S\sqsubseteq T$，当且仅当 $S\land T=S$。

下面的过程为了避免引起混乱，将区分集合运算与关系 $\cup,\cap,\subseteq$ 和 bitmask 的运算与关系 $\land,\lor,\sqsubseteq$，不过 $|x|$ 将同时表示集合 $x$ 的大小和 bitmask $x$ 的 popcount。

刻画两个集合的按位 and 相等极其困难，不妨容斥，不妨设：

$$\begin{aligned} &F(S,T)=\sum_{\substack{P\subseteq U,Q\subseteq U\\P\cap Q=\varnothing}} \left[\bigwedge_{k\in P} k=S,\bigwedge_{k\in Q}k=T\right]\prod_{k\in P\cup Q} a_k\\ &G(S,T)=\sum_{\substack{P\subseteq U,Q\subseteq U\\P\cap Q=\varnothing}} \left[S\sqsubseteq \left(\bigwedge_{k\in P} k\right),T\sqsubseteq\left(\bigwedge_{k\in Q}k\right)\right]\prod_{k\in P\cup Q} a_k\\ &H(S,T)=\sum_{\substack{P\subseteq U,Q\subseteq U\\P\cap Q=\varnothing}} \left[\left(\bigwedge_{k\in P} k\right)=S,T\sqsubseteq\left(\bigwedge_{k\in Q}k\right)\right]\prod_{k\in P\cup Q} a_k \end{aligned}$$

那么利用超集 / 子集反演，可以得到：

$$\begin{aligned} &G(S,T)=\sum_{S\sqsubseteq P}H(P,T)\implies H(S,T)=\sum_{S\sqsubseteq P} (-1)^{|P|-|S|} G(P,T)\\ &H(S,T)=\sum_{T\sqsubseteq Q} F(S,Q)\implies F(S,T)=\sum_{T\sqsubseteq Q} (-1)^{|Q|-|T|} H(S,Q)\\ &\implies F(S,T)=\sum_{S\sqsubseteq P}\sum_{T\sqsubseteq Q} (-1)^{|P|+|Q|-|S|-|T|} G(P,Q)\\ &\implies F(S,S)=\sum_{S\sqsubseteq P}\sum_{S\sqsubseteq Q} (-1)^{|P|+|Q|} G(P,Q) \end{aligned}$$

答案就是：

$$\sum_{S=0}^{2^n-1}F(S,S)=\sum_{S=0}^{2^n-1}\sum_{S\sqsubseteq P}\sum_{S\sqsubseteq Q}(-1)^{|P|+|Q|} G(P,Q)=\sum_{P=0}^{2^n-1}\sum_{Q=0}^{2^n-1} (-1)^{|P|+|Q|} 2^{|P\land Q|} G(P,Q)$$

考察 $G(S,T)$，设 $\tau(S)=\{x\mid S\sqsubseteq x,x\in U\}$ 也就是位运算意义下的“超集”集合。

那么对于 $\tau(A)\cap \tau(B)=\tau(A\lor B)$ 集合的数，既可以放在 $P$，也可以放在 $Q$（但是为了维持不交性质，不能两个都不放），其余的只能放到对应的集合中去。因此可以立即写出式子：

$$G(S,T)=\prod_{k\in \tau(S)}(1+a_k)\cdot \prod_{k\in \tau(T)}(1+a_k)\cdot \prod_{k\in\tau(S\lor T)} \frac{1+2a_k}{(1+a_k)^2}$$

换元 $f(S)=\prod_{k\in S}(1+a_k),g(S)=\prod_{k\in S} \frac{1+2a_k}{(1+a_k)^2}$，那么可以将答案化简：

$$\begin{aligned} &\sum_{P=0}^{2^n-1}\sum_{Q=0}^{2^n-1} (-1)^{|P|+|Q|} 2^{|S\land T|} G(P,Q)\\ &=\sum_{S=0}^{2^n-1}\sum_{T=0}^{2^n-1} (-1)^{|S|+|T|}2^{|S\land T|}f(\tau(S))f(\tau(T))g(\tau(S\lor T)) \end{aligned}$$

注意到 $|S\land T|+|S\lor T|=|S|+|T|$，因此得到：

$$\begin{aligned} &\sum_{S=0}^{2^n-1}\sum_{T=0}^{2^n-1}(-1)^{|S|+|T|}2^{|S|+|T|-|S\lor T|}f(\tau(S))f(\tau(T))g(\tau(S\lor T))\\ =&\sum_{A=0}^{2^n-1}2^{-|A|}g(\tau(A))\sum_{S\lor T=A} (-2)^{|S|}f(\tau(S))\cdot (-2)^{|T|}f(\tau(T)) \end{aligned}$$

对于 $f\circ \tau,g\circ \tau$，可以看成高维后缀积，利用 Sum over Subsets DP $O(n2^n)$ 计算出来，而原式可以看成一个 OR 卷积的形式，使用 FMT 做即可，时间复杂度也是 $O(n 2^n)$，故总时间复杂度 $O(Tn 2^n)$。

你以为你就做完了吗？不！由于不保证 $a_i\neq p-1$，所以 $g(S)$ 可能无意义。这时不妨将所有 $f(S)$ 中的 $0$ 和 $g(S)$ 除数中的 $0$ 换成一个变量 $x$，这样就变成了一个多项式，照常做高维后缀积以及 FMT，那么我们最终会归约到下面的一个问题：

$$\lim_{x\to 0} \frac{F(x)}{G(x)}$$

其中 $F(x),G(x)$ 是多项式函数。在本题中 $F(x)$ 的最低次项一定比 $G(x)$ 的高或相等（如果更低就极限不存在了，而在本题中极限是存在的），那么设 $G(x)$ 的最低次项为 $p$，则可以做 $p$ 次 L'Hospital 法则，得到：

$$\lim_{x\to 0} \frac{F(x)}{G(x)}=\lim_{x\to 0} \frac{F^{(p)} (x)}{G^{(p)} (x)}=\lim_{x\to 0} \frac{F^{(p)}(x)}{p!\cdot [x^p]G(x)}$$

若 $F(x)$ 的最低次项 $>p$，那么 $\lim_{x\to 0} F^{(p)}(x)=0$，原式 $=0$，否则原式 $=\frac{[x^p] F(x)}{[x^p] G(x)}$。因此实际上我们不需要计算出整个多项式，只需要记录多项式的最低次项即可，这样时间复杂度不变。

实现

#include <bits/stdc++.h>
#define popcnt __builtin_popcount
using std::cin, std::cout;

constexpr int mod = 998244353, inv2 = (mod + 1) >> 1;
int A(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int S(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int M(int x, int y){ return 1ll * x * y % mod; }
int P(int x, int y){ int ans = 1; for(;y;x=M(x,x),y>>=1) if(y & 1) ans = M(ans, x); return ans; }
int I(int x){ return P(x, mod - 2); }
template<class T> int Fit(T x){ return (x % mod + mod) % mod; }

namespace solution{
    const int N = 25, N_ = (1 << 20) + 5;
    int n, a[N_];
    struct ZP {
        int val, deg;
        ZP(int v = 0, int d = 0) : val(v), deg(d) {}
        ZP operator-() const { return {S(0, val), deg}; }
        ZP operator+(const ZP& rhs) const {
            ZP res = deg < rhs.deg ? *this : rhs;
            if(deg == rhs.deg) res.val = A(val, rhs.val);
            return res;
        }
        ZP operator*(const ZP& rhs) const { return {M(val, rhs.val), deg + rhs.deg}; }
    } f[N_], g[N_];
    void solve(){
        cin >> n;
        for(int i=0;i<(1<<n);i++){
            cin >> a[i];
            if(a[i] != mod - 1){
                int p = a[i] + 1, pi = I(p);
                f[i] = {p, 0}, g[i] = {M(A(p, a[i]), M(pi, pi)), 0};
            }
            else f[i] = {1, 1}, g[i] = {mod - 1, -2};
        }
        for(int i=1;i<(1<<n);i<<=1){
            for(int j=0;j<(1<<n);j++) if(j & i) f[j ^ i] = f[j ^ i] * f[j], g[j ^ i] = g[j ^ i] * g[j];
        }
        for(int i=0;i<(1<<n);i++) f[i] = f[i] * ZP(P(mod - 2, popcnt(i))), g[i] = g[i] * ZP(P(inv2, popcnt(i)));
        for(int i=1;i<(1<<n);i<<=1){
            for(int j=0;j<(1<<n);j++) if(j & i) f[j] = f[j] + f[j ^ i];
        }
        for(int i=0;i<(1<<n);i++) f[i] = f[i] * f[i];
        for(int i=1;i<(1<<n);i<<=1){
            for(int j=0;j<(1<<n);j++) if(j & i) f[j] = f[j] + (-f[j ^ i]);
        }
        int ans = 0;
        for(int i=0;i<(1<<n);i++) if(f[i].deg == -g[i].deg) ans = A(ans, M(f[i].val, g[i].val));
        cout << ans << '\n';
    }
    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int c, T; cin >> c >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_13275_NOI_2025_集合.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] P13352 [GDCPC 2024] 腊肠披萨1024 MB

给定 $L$ 个字符串 $s_1, \dots, s_L$ 和两个整数 $C, P$。求 $\sum_{i=1}^L \sum_{j=1}^L C^{|LCPS(s_i, s_j)|} \bmod P$ 的值。其中 $LCPS(s_i, s_j)$ 表示字符串 $s_i$ 的最长前缀和字符串 $s_j$ 的最长后缀相同且最长的字符串。

$1\leq L,\sum |s_i| \leq 10^6$。

一个非常直观的做法是每次我们枚举字符串的一个前缀，然后尝试找到这 $L_i$ 个字符串中有多少个后缀，这个前缀恰好等于这个后缀，那么这个东西显然可以通过 Hash 实现。

但是这个做法会有一个问题，我们会算重。完全可能存在 $S[1,i]=T[|T|-i+1,|T|]$ 同时又有 $S[1,j]=T[|T|-j+1,|T|]$ 的情况（假设 $i<j$）。但是这样就能导出一个事实，即 $S[1,i]=T[|T|-j+1,|T|-j+i]$，或者说，$S[1,j]$ 的长度为 $i$ 的前缀与长度为 $i$ 的后缀是相等的。也就是说，$S[1,j]$ 有一个长度为 $i$ 的 Border。刻画 Border 的首选技术是 KMP / ACAM。当然这道题如果用 ACAM 就是有点大炮打蚊子的意思。

用 KMP 求出每一个字符串的前缀函数，那么我们就知道每个字符串的前缀的最长 Border 长度，直接容斥，算这一个前缀的贡献时，减去最长 Border 的贡献（因为它肯定算过了），就不会算重了。时间复杂度 $O(\sum |s_i|)$。

实现

#include <bits/stdc++.h>
#include <bits/extc++.h>
using std::cin, std::cout;

namespace solution{
    using u64 = unsigned long long;
    constexpr int N = 3e6 + 5;
    constexpr u64 hash_base = 1145141;
    int n, C, mod, Cpw[N];
    std::string s[N];
    std::vector<int> fail[N];
    u64 base[N]; __gnu_pbds::gp_hash_table<u64, int> cnt;
    int A(int x, int y){ return (x + y) >= mod ? (x + y) - mod : x + y; }
    int S(int x, int y){ return (x - y) < 0 ? (x - y) + mod : x - y; }
    void solve(){
        cin >> n >> C >> mod;
        base[0] = Cpw[0] = 1;
        for(int i=1;i<N;i++) base[i] = base[i - 1] * hash_base, Cpw[i] = 1ll * Cpw[i - 1] * C % mod;
        for(int i=1;i<=n;i++){
            cin >> s[i], fail[i].resize(s[i].size() + 1), s[i] = ' ' + s[i];
            for(size_t r=2,l=0;r<s[i].size();r++){
                while(l && s[i][l + 1] != s[i][r]) l = fail[i][l];
                if(s[i][l + 1] == s[i][r]) l++;
                fail[i][r] = l;
            }
            u64 hsh = 0; cnt[0]++;
            for(size_t j=s[i].size()-1;j>=1;j--) cnt[hsh = hsh + base[s[i].size() - j - 1] * (s[i][j] - 'a' + 1)]++;
        }
        int ans = 0;
        for(int i=1;i<=n;i++){
            u64 hsh = 0;
            for(size_t j=0;j<s[i].size();j++){
                if(j) hsh = hsh * hash_base + (s[i][j] - 'a' + 1), ans = S(ans, 1ll * cnt[hsh] * Cpw[fail[i][j]] % mod);
                ans = A(ans, 1ll * cnt[hsh] * Cpw[j] % mod);
            }
        }
        cout << ans << '\n';
    }

    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [P_13352_GDCPC_2024_腊肠披萨.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] P13360 [GDCPC 2024] 另一个计数问题3000 ms

给定一个 $n - 1$ 个点的无向图，点的编号为 $2 \sim n$。对于所有的 $2 \le u < v \le n$，边 $(u, v)$ 存在当且仅当 $v$ 是 $u$ 的正整数倍。定义 $f(u, v)$ 表示 $u$ 与 $v$ 是否连通：当 $u, v$ 连通时 $f(u, v) = 1$，否则 $f(u, v) = 0$。求：

$$\left(\sum_{u = 2} ^ {n - 1} \sum_{v = u + 1} ^ n f(u, v) \cdot u \cdot v\right) \bmod {998244353}$$

$4\leq n\leq 10^{11}$。

下面记：

$$\begin{aligned} &S_1(n)=\sum_{i=1}^{n} i=\frac{1}{2}n(n+1)\\ &S_2(n)=\sum_{i=1}^{n} i^2=\frac{1}{6}n(n+1)(2n+1)\\ &P_k(n)=\sum_{p\in\mathbb{P},p\leq n} p^k\\ \end{aligned}$$

注意对于任意 $x\not\in\mathbb{P}$，取 $p$ 为 $x$ 的最小质因子，那么 $p,x$ 连通。又取任意 $p\in\mathbb{P},p\leq \frac{n}{2}$，则有 $2p\leq n$，$p$ 与 $2p$ 连通，而 $2p$ 与 $2$ 连通，故所有合数以及 $\leq \frac{n}{2}$ 的质数，都与 $2$ 连通。而 $>\frac{n}{2}$ 的质数，在 $[1,n]$ 中不存在大于本身的倍数，所以它们是孤立点。于是原问题就相当于我们要求（记 $F_k=P_k(n)-P_k(\left\lfloor \frac{n}{2} \right\rfloor$）：

$$\begin{aligned} &\sum_{u=2}^{n-1}\sum_{v=u+1}^n uv-\sum_{u\in (\frac{n}{2},n]\cap\mathbb{P}} u(S_1(n)-1-F_1)-\sum_{u\in (\frac{n}{2},n]\cap\mathbb{P}}\sum_{v\in (u,n]\cap\mathbb{P}} uv\\ =&\frac{1}{2}\left(S_1^2(n)-S_2(n)-2S_1(n)+1\right)-F_1(S_1(n)-1-F_1)-\frac{1}{2}(F_1^2-F_2) \end{aligned}$$

因此只要快速求出来 $F_1,F_2$ 就可以快速求出本题的答案。而 $P_1,P_2$ 是经典的质数前缀统计问题，可以用 Min_25 筛的前半部分做，时间复杂度 $O(\frac{n^{3/4}}{\log n})$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

using i64 = long long;
constexpr int mod = 998244353, inv6 = 166374059, inv2 = 499122177;
int A(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int S(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int M(int x, int y){ return 1ll * x * y % mod; }
int P(int x, int y){ int ans = 1; for(;y;x=M(x,x),y>>=1) if(y & 1) ans = M(ans, x); return ans; }
int I(int x){ return P(x, mod - 2); }
template<class T> int Fit(T x){ return (x % mod + mod) % mod; }

namespace solution {
    using i64 = long long;
    const int N = 1e6 + 5;
    int pri[N], tot, F[N], G[N];
    bool vis[N];

    void sieve(int n){
        for(int i=2;i<=n;i++){
            if(!vis[i]){
                pri[++tot] = i;
                F[tot] = A(F[tot - 1], i);
                G[tot] = A(G[tot - 1], M(i, i));
            }
            for(int j=1;j<=tot&&1ll*i*pri[j]<=n;j++){
                vis[i * pri[j]] = 1;
                if(!(i % pri[j])) break;
            }
        }
    }

    i64 n, w[N];
    int f[N], g[N], sqn, wtt, rev[N], revt[N];
    int get(i64 x){ return x <= sqn ? rev[x] : revt[n / x]; }
    int S1(i64 n){ return n = Fit(n), M(M(n, A(n, 1)), inv2); }
    int S2(i64 n){ return n = Fit(n), M(M(M(n, A(n, 1)), A(A(n, n), 1)), inv6); }

    void min25(){
        sqn = std::sqrt(n), sieve(sqn);
        for(i64 l=1,r;l<=n;l=r+1){
            r = n / (n / l), w[++wtt] = n / l;
            if(w[wtt] <= sqn) rev[w[wtt]] = wtt;
            else revt[n / w[wtt]] = wtt;
            f[wtt] = S(S1(w[wtt]), 1),  g[wtt] = S(S2(w[wtt]), 1);
        }
        for(int j=1;j<=tot;j++){
            for(int i=1;i<=wtt&&1ll*pri[j]*pri[j]<=w[i];i++){
                f[i] = S(f[i], M(pri[j], S(f[get(w[i] / pri[j])], F[j - 1])));
                g[i] = S(g[i], M(M(pri[j], pri[j]), S(g[get(w[i] / pri[j])], G[j - 1])));
            }
        }
    }

    void solve(){
        cin >> n, min25();
        int ans = A(S(M(S(M(S1(n), S1(n)), S2(n)), inv2), S1(n)), 1);
        int F = S(f[1], f[2]), G = S(g[1], g[2]);
        ans = S(ans, M(S(S(S1(n), 1), F), F));
        ans = S(ans, M(S(M(F, F), G), inv2));
        cout << ans << '\n';
    }

    void clear() {}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [J_另一个计数问题.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium] P13351 [GDCPC 2024] 田字格3000 ms

给定 $n$ 条线段，仅包含垂直线段（$x_1=x_2$）和水平线段（$y_1=y_2$），保证任意两条垂直线段 / 水平线段无交。

你需要计算有多少个田字格。定义一个田字格 $(x_0,y_0,d)$ 表示存在点 $(x_0,y_0)$ 和边长 $d$，使得区域 $U=[x_0-d,x_0+d]\times [y_0-d,y_0+d]$ 存在且仅存在六条线段，且这六条线段分别为 $x=x_0-d,x=x_0,x=x_0+d,y=y_0-d,y=y_0,y=y_0+d$ 分别与 $U$ 的交。

$1\leq n\leq 3\times 10^5$。

田字格可以拆成两个部分：

• 三条垂直直线 $x=x_0-d,x=x_0,x=x_0+d$ 与 $U$ 的交。
• 三条水平直线 $y=y_0-d,y=y_0,y=y_0+d$ 与 $U$ 的交。

先来考虑怎样的三条垂直的线段可能可以贡献答案，可以发现它应当满足：

1. 两两相邻线段的距离应当均为定值 $d$。
2. 任意两条相邻线段间没有其他垂直线段。
3. 设三条线段的上端点的最小纵坐标为 $R$，下断点的最大纵坐标是 $L$，则 $R-L\geq 2d$。

考虑按照纵坐标扫描线，当一条垂直线段被加入 / 删除时统计三条垂直线段组。那么这样的线段组可以用 std::set 维护，这样我们可以方便地寻找前驱后继来满足条件 $2$，找到后验证是否满足条件 $1$ 即可。可以发现扫描线的每一轮我们只会找出 $O(1)$ 个（准确来说是 $5$ 个）三条垂直线段组来判断是否满足条件，所以总共只有 $O(n)$ 个满足条件的线段组。

注意这里有一个技巧，可以在生长线和消亡线时都将线段组放进同一个向量中，这样奇数位就是 $L$，偶数为就是 $R$ 非常方便。

水平线段是类似的，现在只需要考虑两个线段组是否可以贡献答案。我们不妨对于每一个线段组，记录 $(d,x,L,R)$ 表示距离、该线段组中间的那条的横坐标 / 纵坐标，以及存在区间 $[L,R]$（上面求出的）。

那么你会发现水平线段组 $(d_0,x_0,L_0,R_0)$ 和垂直线段组 $(d_1,x_1,L_1,R_1)$ 可以贡献，当且仅当：

• $x_1\in [L_0+d_0,R_0-d_0]$。
• $x_0\in [L_1+d_1,R_1-d_1]$。
• $d_0=d_1$。

于是可以按照 $d$ 分组，组内离散化后扫描线 + 树状数组即可。时间复杂度 $O(n\log n)$。

实现

#include <bits/stdc++.h>
using std::cin, std::cout;

namespace solution{
    using i64 = long long;
    const int N = 3e5 + 5;
    int n;
    struct triple {
        int x, y, z;
        bool operator<(const triple& rhs) const {
            return std::make_tuple(x, y, z) < std::make_tuple(rhs.x, rhs.y, rhs.z);
        }
    };
    std::basic_string<triple> a[2], L;
    std::map<int, std::basic_string<triple>> Rg[2];
    struct BIT {
        i64 t[N];
        void update(int p, int v){
            while(p <= n) t[p] += v, p += p & -p;
        }
        i64 query(int p){
            i64 res = 0;
            while(p) res += t[p], p -= p & -p;
            return res;
        }
    } T;
    void solve(){
        cin >> n;
        for(int i=1,x1,y1,x2,y2;i<=n;i++){
            cin >> x1 >> y1 >> x2 >> y2;
            if(x1 == x2) a[0] += {x1, y1, y2};
            else a[1] += {y1, x1, x2};
        }
        for(int _ : {0, 1}){
            L.clear();
            for(auto [i, l, r] : a[_]) L += {l, 1, i}, L += {r, -1, i};
            std::sort(L.begin(), L.end());
            std::set<int> S;
            std::map<std::pair<int, int>, std::basic_string<int>> semi;
            for(auto tr : L){
                int x = tr.x, o = tr.y, i = tr.z;
                auto chk = [&](int l, int m, int r){ if(r - m == m - l) semi[{m, m - l}] += x; };
                auto stg0 = [&](auto it){
                    if(it != S.begin() && it != S.end() && std::next(it) != S.end()) chk(*std::prev(it), *it, *std::next(it));
                    if(it != S.end() && it != S.begin() && std::prev(it) != S.begin()) chk(*std::prev(it, 2), *std::prev(it), *it);
                };
                auto stg1  = [&](auto it){
                    if(it != S.begin() && std::next(it) != S.end()) chk(*std::prev(it), *it, *std::next(it));
                    if(it != S.begin() && std::prev(it) != S.begin()) chk(*std::prev(it, 2), *std::prev(it), *it);
                    if(std::next(it) != S.end() && std::next(it, 2) != S.end()) chk(*it, *std::next(it), *std::next(it, 2));
                };
                if(~o) stg0(S.upper_bound(i)), stg1(S.insert(i).first);
                else stg1(S.find(i)), S.erase(i), stg0(S.upper_bound(i));
            }
            for(auto [pd, v] : semi){
                auto [p, d] = pd;
                for(std::size_t i=0;i<v.size();i+=2){
                    if(v[i + 1] - v[i] >= (d << 1)) Rg[_][d] += {p, v[i], v[i + 1]};
                }
            }
        }
        i64 ans = 0;
        for(auto [d, m0] : Rg[0]){
            if(!Rg[1].count(d)) continue;
            auto m1 = Rg[1][d];
            L.clear();
            std::basic_string<int> buf;
            for(auto [p, l, r] : m0) buf += p;
            for(auto [p, l, r] : m1) buf += r - d, buf += l + d - 1;
            std::sort(buf.begin(), buf.end());
            buf.erase(std::unique(buf.begin(), buf.end()), buf.end());
            for(auto [p, l, r] : m0) p = std::lower_bound(buf.begin(), buf.end(), p) - buf.begin() + 1, L += {l + d, 1, p}, L += {r - d + 1, -1, p};
            for(auto [p, l, r] : m1){
                int rr = std::lower_bound(buf.begin(), buf.end(), r - d) - buf.begin() + 1;
                int ll = std::lower_bound(buf.begin(), buf.end(), l + d - 1) - buf.begin() + 1;
                L += {p, 2, rr}, L += {p, 3, ll};
            }
            std::sort(L.begin(), L.end());
            for(auto [x, o, i] : L){
                if(o == -1) T.update(i, -1);
                else if(o == 1) T.update(i, 1);
                else if(o == 2) ans += T.query(i);
                else ans -= T.query(i);
            }
        }
        cout << ans << '\n';
    }
    void clear(){}
}

signed main(){
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    int T = 1; // cin >> T;
    while(T--) solution::solve(), solution::clear();
    return 0;
}

// [A_田字格.cpp] Think twice, code once.
// Written by xiezheyuan

[Medium-Hard] P3266 [JLOI2015] 骗我呢

给定一个 $n \times m$ 的矩阵 $x_{i,j}$，其中每个元素满足 $0 \le x_{i,j} \le m$。求满足以下条件的矩阵的数量，结果对 $10^9+7$ 取模：

1. 对于所有 $1 \le i \le n, 1 \le j < m$，有 $x_{i,j} < x_{i,j+1}$。
2. 对于所有 $1 < i \le n, 1 \le j < m$，有 $x_{i,j} < x_{i-1,j+1}$。

$1\leq n,m\leq 10^6$。

先来了解一个经典的问题：

在平面直角坐标系上有两条直线 $y=x+a$ 和 $y=x+b$，你从 $(0,0)$ 出发，每次只能向右走或向上走一个单位，但是不能与这两条直线相交（允许相切），求走到 $(n,m)$ 的方案数。

在解决这个经典问题之前，不妨先考虑解决一个更加经典的问题：

在平面直角坐标系上有一条直线 $y=x+a$，你从 $(0,0)$ 出发，每次只能向右或者向上走一个单位，但是不能与这条直线相交（允许相切），求走到 $(n,m)$ 的方案数。