P4642 题解

· · 题解

发现一道没有题解的题。

怎么辉石呢？是啊怎么辉石呢？

题意

给定 A,B,C，求：

& \max z = \sum_{i = 1}^{n} C_i x_i\\ &s.t \begin{cases} \sum_{i = 1}^{n} A_i x_i = S\\ \sum_{i = 1}^{n} B_i x_i = T \end{cases} \end{aligned}

解法

虽然已经是一个线性规划的标准形式，但是可以发现这个玩意就非常的屑，把它直接当成线性规划做肯定会 TLE，想要当成带约束多元函数极值做拉格朗日乘数结果求出偏导数没法反解 x，做高斯消元也会 TLE。

但是有一个重要的条件是线性约束只有两个，考虑转化问题，使用线性规划的对偶原理将问题转化为变量只有两个。

对偶得到如下线性规划：

& \min z = Sx + Ty\\ &s.t \begin{cases} \forall i \in [1,n] A_i x + B_i y \ge C_i\\ x,y \in \mathbb{R} \end{cases} \end{aligned}

线性规划除了单纯形法外还有一个常见解法：图解法，即通过几何作图求得可行域，然后求解。

对于对偶问题，只有两个变量的情况，我们可以把 x,y 当作平面上的两维，每个限制 A_i x + B_i y - C_i \ge 0 可以看作直线的一般式方程，分割平面得到的可行域为一个 半平面，问题转化为在这些半平面的交集中找到一点使 Sx + Ty 最小。

我们可以首先离线询问，然后按照 S,T 排序最后统一处理即可。

时间复杂度 \mathcal{O} (n\log n)

代码

去掉了读写优化的部分。

typedef std::pair<f80,f80> Point;

// 半平面的每个点
Point p[N];

struct Line {
    i64 a, b, c;
}l[N],q[M];

inline bool cmp(const Line &x, const Line &y) {
    i64 tmp = x.a * y.b - x.b * y.a;
    if(tmp != 0) return tmp > 0;
    tmp = x.c * y.b - x.b * y.c;
    return tmp > 0;
}

int stk[N];

// 直线一般式方程求交点
inline Point Intersect(const Line &x, const Line &y) {
    i64 d = x.a * y.b - x.b * y.a;
    f80 px = x.c * y.b - x.b * y.c;
    f80 py = x.a * y.c - x.c * y.a;
    return std::make_pair(px / d,py / d);
}

inline bool PointOnLine(const Line &pl,const Point &pp) {
    f80 tmp = pl.a * pp.first + pl.b * pp.second - pl.c;
    if(fabs(tmp) < eps) return 1;
    return !(tmp < 0);
}

f64 ans[M];

void SolveSpec(int m) {
    // 只有一个限制条件的时候需要特判，不然会出现 nan
    while (m--) {
        int S = readI(),T = readI();
        // 平行时无解
        if(S * l[1].b != T * l[1].a)
            PutStr("IMPOSSIBLE");
        else    
            writeF((f64)l[1].c * S / l[1].a,5),enter;
    }
}

void Solve(int n,int m) {
    // 半平面交
    int tp = 2;
    stk[1] = 1,stk[2] = 2;
    p[1].first = 0,p[1].second = (f80)l[1].c / l[1].b;
    p[2] = Intersect(l[1],l[2]);

    rep(i,3,n) {
        while(tp > 1) {
            int pre = stk[tp - 1],top = stk[tp];
            if(PointOnLine(l[top],Intersect(l[i],l[pre]))) --tp;
            else break;
        }
        p[tp + 1] = Intersect(l[i],l[stk[tp]]);
        stk[++tp] = i;
    }

    rep(i,1,m) {
        q[i].a = readI();
        q[i].b = readI();
        q[i].c = i;
        ans[i] = -2;
    }
    std::sort(q + 1,q + m + 1,cmp);

    rep(i,1,m) {
        if(q[i].a * l[1].b > q[i].b * l[1].a)
            ans[q[i].c] = -1;
        else break;
    }
    int cur = 1;
    rep(i,1,m) {
        if(ans[q[i].c] == -1) continue;
        int pl = stk[tp];
        if(q[i].a * l[pl].b < q[i].b * l[pl].a)
            break;
        while(cur <= tp) {
            pl = stk[cur];
            i64 d = q[i].a * l[pl].b - q[i].b * l[pl].a;
            if(d >= 0) {
                ans[q[i].c] = q[i].a * p[cur].first + q[i].b * p[cur].second;
                break;
            }
            else {
                ++cur;
                continue;
            }
        }
        if(cur > tp) break;
    }

    rep(i,1,m) if(ans[i] > 0)
        writeF(ans[i],5),enter;
    else
        PutStr("IMPOSSIBLE");
}

int main() {
    InitIO();
    int n = readI(),m = readI();
    rep(i,1,n) {
        l[i].a = readI();
        l[i].b = readI();
        l[i].c = readI();
    }
    std::sort(l + 1,l + n + 1,cmp);
    int siz = 1;
    rep(i,2,n) {
        if(l[i].a * l[i - 1].b == l[i - 1].a * l[i].b)
            continue;
        l[++siz] = l[i];
    }
    if(siz == 1) SolveSpec(m);
    else Solve(siz,m);
    EndIO();
    return 0;
}

参考

董克凡 浅谈线性规划与对偶问题