题解：P8501 [NOI2022] 二次整数规划问题

· · 题解

首先，约束条件可能和固有上下界产生冲突，导致区间内有些数不可能被取到，我们直接类似于 SPFA 跑 n 轮松弛一下这些约束就行了，这样子每个数在取值区间内的数都可能被取到。然后固定一些 l_i=r_i 的数，还有一些数 l_i\neq r_i，但是 1\in [l_i,r_i]，5\in [l_i,r_i]，这个时候就把 1,5 从值域中踢出去，原因会在下一段提及。

最大化的式子中可以发现，首先 v_1 和 v_k 无贡献，其次取 1,k 的时候对于 G 的贡献没有那么优。所以通过调整法，我们可以发现选 1 或选 k 不如选 [2,k-1]。除了那些必须选 1 和 k 的位置外，值域可以缩小到 [2,k-1]。

对于部分分，k=3,4 都可以运用基础调整法迅速得到最优解的形式。对于 k=3 的时候，能选 2 就选 2。对于 k=4 的时候，贡献函数关于未知数为线性形式，所以那些未确定数要么全选 2，要么全选 3。

对于 k=5 的时候。建议再读一遍第一段的第一步变量取值处理，然后再看这段文字。未确定数的取值有 \{2,3,4\}。首先思考约束关系，那些必定取某个值也就是 l_i=r_i 的位置如果对于未确定数产生约束的时候肯定在松弛的时候已经产生过了。比较难限制的是动区间的约束，比如两个取值为 [2,4] 和 [2,4] 的变量，要求距离 \le 1，也就是不能一个 2，另一个 4，在松弛那一步约束不到这里。可以发现在 \{2,3,4\} 集合内，约束为 1 的条件传递两次没有意义，因为传递两次之后就变成了距离约束为 2，而集合内部的最大距离也就是 2，所以没有用了。但是距离为 0 的约束传递多次还是紧的，可以用并查集处理。

于是，总结一下变量取值范围的约束大部分在松弛那一步已经处理完了，还有两个需要处理，一个是两个变量要求距离 \le 1（这个两个变量只可能是 m 条约束中的某个），还有一个是一些变量要求相等 （不一定为 m 条约束之一，可能是通过并查集传递的）。

处理完变量取值范围之后，我们回到最优化问题上来。可以发现贡献函数形式必定如下所示，其中 V_i=v_i+?，这个 ? 代表一些距离大于 1 的变量贡献需要扣掉，反正是关于未知数是一个线性关系。

W=10^6(n^2-2c_2c_4)+\sum\limits_{i=2}^4c_iV_i+\lambda

化简式子，就是就是形如 -M(c_2-A_1)(c_4-A_2)+A_3。

令 (c_2,c_4)\to (x,y)，我们需要最小化 (x-A_1)(y-A_4)。

其中可以发现目前已知三个点肯定能取到，分别是 (x_{\min},y_{\min})，(x_{\min},y_{\max}) 和 (x_{\max},y_{\min})。x_{\min} 就是只有那些松弛过后 l_i=r_i 的点产生贡献，x_{\max} 则是把所有 2\in[l_i,r_i] 的点都压到了 2，保证能取到是因为这些动区间点聚集在一起之后距离肯定为 0 也就肯定满足约束，而剩下点还是那一句话定点对于动区间点的约束已经在松弛一步约束过了，所以这三个点肯定能取到。

我们现在需要处理出一个可能取到最优答案的点集，遇到哪个询问，就直接对于点集内所有点代入算一遍看看哪个最优就行了。你直接把所有可能的合法点放入这个点集，显然复杂度是要爆的。因此我们需要挑出一些有用的点放入，一些 100\% 不可能产生贡献的点就不放进去了。

题目要处理的就是平移之后的 x'y' 乘积最小值。可以发现如果平移之后，还存在点处于 1,2,4 象限，那么肯定包含上面三个定点之一，选择它们是最好的（这个可以自己画画图就能看出来）。因此先把这三个点纳入点集。

问题是如果所有点都被移动到了第三象限怎么办？乘法之后负负得正，我们可以作出判断，只有上凸壳上面的点是有效的（注意如果是 x,y>0 那么是下凸壳，但是这里 x,y<0，所以是上凸壳）。

根据引理，

在一个 V\times V 的矩形上，凸包点的个数是 O(V^{\frac{2}{3}}) 的。

于是我们就把有用点点集大小缩小到了凸包大小量级，也就是 O(V^{\frac{2}{3}})。

可惜我们无法直接显式建立这个凸包。不过根据 P5540 的 trick，我们可以解决这个凸包问题。我们知道凸包上的两点 P(x_{\min},y_{\max}) 和 Q(x_{\max},y_{\min})。我们可以通过已有的两点去构建这个凸包，直接寻找所有点还是不容易，考虑哪个点是一定在凸包上的，那么就是一个点 R 满足在直线 PQ 上方且离 PQ 距离最大。找到 R 点之后更新答案，并且分治为 (P,R) 和 (R,Q) 去解决问题，这样子我们每次找到一个肯定在凸包上的点，总的寻找复杂度也就是凸包上点的个数乘以单次寻找时间复杂度。

现在考虑如何找到 R 点，上述条件等价于三角形 PQR 的面积最大，也就是 \vec{RP}\times \vec{RQ} 最大。代入坐标，利用叉乘公式计算之后可以得到，希望最大化 (X_Q-X_P)y_R+(Y_P-Y_Q)x_R，也就是最小化 (X_P-X_Q)y_R+(Y_Q-Y_P)x_R。

也就是说每个数如果选 2 就会有 Y_Q-Y_P 的代价，如果选 4 就会有 X_P-X_Q 的贡献，选 3 没有贡献，同时一些数之间会有在上文中的两个约束。这显然是一个切糕模型的形式。每次利用切糕模型建图，代价可能是负的，平移 n 即可。跑最小割即可。利用最小割的流量信息求出一组方案，然后找到凸壳上的点，继续分治就行了。注意找到点之后要验证其是否在凸包上，以防复杂度退化。

时间复杂度 O(n^{\frac{2}{3}}(q+\rm Dinic()))。

#include<bits/stdc++.h>
#define pb emplace_back
#define fi first
#define se second
#define mp make_pair
using namespace std;
typedef long long ll;
typedef pair<int,int> pii; 
void cmax(int &x,int y){ x=x>y?x:y; }
void cmin(int &x,int y){ x=x<y?x:y; }
const int maxn=610;
const int maxm=1e6+10;
const ll Vg=1e6;
const ll inf=1e13;
int L[maxn],R[maxn],n,K,m,q; ll val[maxn];
int cnt[6],U[maxm],V[maxm],W[maxm];
vector<pii> vec; int id[maxn][2],c=0;
struct DSU{
    int fa[maxn],sz[maxn];
    void init(){ for(int i=1;i<=n;i++) fa[i]=i,sz[i]=1; }
    int find(int u){ return fa[u]==u?u:fa[u]=find(fa[u]); }
    bool merge(int u,int v){
        u=find(u); v=find(v);
        if(u==v) return 0;
        if(sz[u]<sz[v]) swap(u,v);
        fa[v]=u; sz[u]+=sz[v]; return 1;
    }
}dsu;
struct Dinic{
    struct edge{
        int to,Next; ll f;
    }edges[maxm]; queue<int> q;
    int head[maxm],now[maxm],d[maxm],tot;
    void init(){ tot=1; for(int i=1;i<=c+2;i++) head[i]=0; }//别忘记调用了 
    void add(int u,int v,ll f){
        edges[++tot]=(edge){v,head[u],f}; head[u]=tot;
        edges[++tot]=(edge){u,head[v],0}; head[v]=tot;
    }
    bool bfs(int s,int t){
        while(q.size()) q.pop(); for(int i=1;i<=c+2;i++) d[i]=0; 
        d[s]=1; q.push(s); now[s]=head[s];//清空 now(s) 
        while(!q.empty()){
            int u=q.front(); q.pop();
            for(int i=head[u];i;i=edges[i].Next){
                int v=edges[i].to; if(d[v]||!edges[i].f) continue;//注意两个条件 
                now[v]=head[v]; d[v]=d[u]+1; q.push(v);//清空 now(v) 
                if(v==t) return 1; 
            }
        }
        return 0;
    }
    ll dinic(int s,int t,ll flow){
        if(s==t) return flow;
        ll rest=flow,k;
        for(int i=now[s];i&&rest;i=edges[i].Next){//判断 rest 是否有剩余 
            now[s]=i; int v=edges[i].to;//动态更新 now 
            if(!edges[i].f||d[v]!=d[s]+1) continue;// 注意 d 的条件 
            k=dinic(v,t,min(rest,edges[i].f));
            if(!k) d[v]=0; //注意 
            rest-=k; edges[i].f-=k; edges[i^1].f+=k;
        }
        return flow-rest;
    }
    pii solve(int s,int t){
        ll maxflow=0,flow;
        while(bfs(s,t))
            while(flow=dinic(s,t,inf)) maxflow+=flow;
        int c2=0,c4=0; bfs(s,t);
        for(int i=1;i<=n;i++){
            if(dsu.find(i)!=i||L[i]==1||R[i]==5) continue;
            if(!d[id[i][0]]) c2+=dsu.sz[i];
            else if(d[id[i][1]]) c4+=dsu.sz[i];
        }
        return mp(c2,c4);
    }
}DC;
void init(){
    dsu.init(); vec.clear(); c=0;
    for(int i=1;i<=5;i++) cnt[i]=0;
}
namespace K3{
    void solve(){
        ll ans; 
        while(q--){
            cin>>val[2];
            int c2=n-cnt[1]-cnt[3];
            ans=Vg*(1ll*n*n-2ll*cnt[1]*cnt[3])+c2*val[2];
            cout<<ans<<"\n";
        }
    }
}
namespace K4{
    void solve(){
        ll ans; 
        while(q--){
            cin>>val[2]>>val[3];
            int r=n-cnt[1]-cnt[2]-cnt[3]-cnt[4];
            int c2=cnt[2]+r;
            ans=Vg*(1ll*n*n-2ll*cnt[1]*(cnt[3]+cnt[4])-2ll*c2*cnt[4])+c2*val[2]+cnt[3]*val[3];
            int c3=cnt[3]+r;
            ans=max(ans,Vg*(1ll*n*n-2ll*cnt[1]*(c3+cnt[4])-2ll*cnt[2]*cnt[4])+cnt[2]*val[2]+c3*val[3]);
            cout<<ans<<"\n";
        }
    }
}
bool chk(pii A,pii B,pii C){
    pii D=mp(B.fi-A.fi,B.se-A.se);
    pii E=mp(C.fi-A.fi,C.se-A.se);
    if(1ll*D.fi*E.se-1ll*D.se*E.fi>0) return 1;
    else return 0;
}
void solve(pii P,pii Q){
    DC.init(); int S=c+1,T=c+2;
    int X=Q.se-P.se+n,Y=P.fi-Q.fi+n;
    for(int i=1;i<=n;i++){
        if(dsu.find(i)!=i||R[i]==1||L[i]==5) continue;
        int sz=dsu.sz[i];
        DC.add(S,id[i][0],(L[i]==2?sz*X:inf));
        DC.add(id[i][0],id[i][1],(L[i]<=3&&R[i]>=3)?sz*n:inf);
        DC.add(id[i][1],T,(R[i]==4?sz*Y:inf));
    }
    for(int i=1;i<=m;i++){
        int u=dsu.find(U[i]),v=dsu.find(V[i]);
        if(W[i]!=1||u==v) continue;
        DC.add(id[u][1],id[v][0],inf); 
        DC.add(id[v][1],id[u][0],inf);
    }
    pii R=DC.solve(S,T);
    if(!chk(P,Q,R)) return ;
    vec.pb(R); solve(P,R); solve(R,Q); 
}
ll calc(){
    ll res=0;
    for(int i=2;i<=K-1;i++) res+=cnt[i]*val[i];
    for(int i=1;i<=K;i++){
        res+=Vg*cnt[i]*cnt[i];
        if(i+1<=K) res+=Vg*cnt[i]*cnt[i+1];
        if(i) res+=Vg*cnt[i-1]*cnt[i];
    }
    return res;
}
void solve(){
    cin>>K>>n>>m>>q; init();
    for(int i=1;i<=n;i++)
        for(auto j:{0,1})
            id[i][j]=++c;
    for(int i=1;i<=n;i++) cin>>L[i]>>R[i];
    for(int i=1;i<=m;i++) cin>>U[i]>>V[i]>>W[i];
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++){
            int x=U[j],y=V[j],w=W[j];
            cmax(L[x],L[y]-w); cmin(R[x],R[y]+w);
            cmax(L[y],L[x]-w); cmin(R[y],R[x]+w);
        }
    }
    int mx2=0,mx4=0;
    for(int i=1;i<=n;i++){
        if(L[i]<K&&R[i]>1) cmax(L[i],2),cmin(R[i],K-1);
        if(L[i]==R[i]) cnt[L[i]]++;
        mx2+=(L[i]==2); mx4+=(R[i]==4);
    }
    if(K==3){ K3::solve(); return ; } 
    if(K==4){ K4::solve(); return ; }
    for(int i=1;i<=m;i++)
        if(!W[i]) dsu.merge(U[i],V[i]); 
    vec.pb(mp(cnt[2],cnt[4])); vec.pb(mp(mx2,cnt[4])); vec.pb(mp(cnt[2],mx4));
    solve(vec[2],vec[1]);
    while(q--){
        for(auto i:{2,3,4}) cin>>val[i];
        ll ans=0;
        for(auto z:vec){
            cnt[2]=z.fi,cnt[4]=z.se;
            cnt[3]=n-cnt[1]-cnt[2]-cnt[4]-cnt[5];
            ans=max(ans,calc());
        }
        cout<<ans<<"\n";
    }
}
int main(){
    ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
    int c,T; cin>>c>>T;
    while(T--) solve();
    return 0;
}