Hihocoder 1429 A New Ground Heating Device （几何 圆并）

Problem

A brand new photosensitive ground heating device is under developing.

This time, to test these devices and help local farmers, engineers placed several devices on the ground in a greenhouse. (You can assume the ground of the greenhouse as a plane and the height of all devices are zero)

There is only a light source in the greenhouse. Its coordinate at the plane must be (0,0),but you can put it at any height above the ground(or on the ground) as you want.

The effective warming radius of the $i$th device is $\frac{W}{L_i\times Z_i}$ , among which $w$ is the power of the light source, $L_i$ is the distance between light source and the $i$th device, while owing to different degree of wear and tear, photoresistance factor of each device, $Z_i$, may be different.

The winter is so harsh that a piece of land in this greenhouse is suitable for planting only when it is heated by at least $K$ devices at the same time.

Considering the efficiency of production, farmers require that the area of arable land should be greater than $S$.

Engineers wonder, to satisfy the heating demand raised by farmers, what is the maximum height of the light source.

主要意思就是 n 个机器，每个机器制热的覆盖半径是 $\frac{W}{L_i\times Z_i}$ ，其中 W 和 $Z_i$ 均为定值，$L_i$ 指机器与 light source 的三维距离（机器的坐标已固定，且高度一定为 0 ，light source 的横纵坐标均为 0， 高度可变）。求 light source 可行的最高高度，使得满足至少有 S 的面积能被至少 K 个机器制热半径覆盖。

解题思路

二分枚举 light source 的高度，对每次产生的高度判断被 $\ge K$ 个机器覆盖的面积是否大于等于 S。

对于求 K 机器覆盖的面积，套用 K 次圆并的模板。

:cry: 模板题。

代码

#include<bits/stdc++.h>
#define sqr(x) ((x)*(x))
using namespace std;
const int N = 200 + 10;
const double eps = 1e-8;
const double pi = acos(-1.0);
double area[N];
int T, n, W, K, S, x[N], y[N];
double z[N];
int dcmp(double x) {
    if (x < -eps) return -1; else return x > eps;
}
struct cp {
    double x, y, r, angle;
    int d;
    cp(){}
    cp(double xx, double yy, double ang = 0, int t = 0) {
        x = xx;  y = yy;  angle = ang;  d = t;
    }
    void get() {
        scanf("%lf%lf%lf", &x, &y, &r);
        d = 1;
    }
}cir[N], tp[N * 2];
double dis(cp a, cp b) {
    return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
double cross(cp p0, cp p1, cp p2) {
    return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);
}
int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &cp2) {
    double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx;
    double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my;
    double sq = mx2 + my2, d = -(sq - sqr(r1 - r2)) * (sq - sqr(r1 + r2));
    if (d + eps < 0) return 0; if (d < eps) d = 0; else d = sqrt(d);
    double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2;
    double y = my * ((r1 + r2) * (r1 - r2) + my * sy) + sy * mx2;
    double dx = mx * d, dy = my * d; sq *= 2;
    cp1.x = (x - dy) / sq; cp1.y = (y + dx) / sq;
    cp2.x = (x + dy) / sq; cp2.y = (y - dx) / sq;
    if (d > eps) return 2; else return 1;
}
bool circmp(const cp& u, const cp& v) {
    return dcmp(u.r - v.r) < 0;
}
bool cmp(const cp& u, const cp& v) {
    if (dcmp(u.angle - v.angle)) return u.angle < v.angle;
    return u.d > v.d;
}
double calc(cp cir, cp cp1, cp cp2) {
    double ans = (cp2.angle - cp1.angle) * sqr(cir.r) 
        - cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2);
    return ans / 2;
}
void CirUnion(cp cir[], int n) {
    cp cp1, cp2;
    sort(cir, cir + n, circmp);
    for (int i = 0; i < n; ++i)
        for (int j = i + 1; j < n; ++j)
            if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0)
                cir[i].d++;
    for (int i = 0; i < n; ++i) {
        int tn = 0, cnt = 0;
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r,
                cp2, cp1) < 2) continue;
            cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x);
            cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x);
            cp1.d = 1;    tp[tn++] = cp1;
            cp2.d = -1;   tp[tn++] = cp2;
            if (dcmp(cp1.angle - cp2.angle) > 0) cnt++;
        }
        tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, pi, -cnt);
        tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -pi, cnt);
        sort(tp, tp + tn, cmp);
        int p, s = cir[i].d + tp[0].d;
        for (int j = 1; j < tn; ++j) {
            p = s;  s += tp[j].d;
            area[p] += calc(cir[i], tp[j - 1], tp[j]);
        }
    }
}
double getR(double X, double Y, double h, double Z) {
    double l = sqrt(X*X + Y*Y + h*h);
    return W * 1.0 / l / Z;
}
bool jug(double h) {
    memset(cir, 0, sizeof(cir));
    for(int i=0;i<n;i++)
    {
        cir[i].x = x[i];
        cir[i].y = y[i];
        cir[i].r = getR(x[i], y[i], h, z[i]);
        cir[i].d = 1;
    }
    memset(area, 0, sizeof(area));
    CirUnion(cir, n);
    if(dcmp(area[K]-S) >= 0)    return true;
    return false;
}
int main() {
    scanf("%d", &T);
    while(T-- && scanf("%d %d %d %d", &n, &W, &K, &S))
    {
        for(int i=0;i<n;i++)
            scanf("%d %d %lf", &x[i], &y[i], &z[i]);
        if(jug(500)) {
            printf("Oops!\n");
        } else if(jug(0) == false) {
            printf("No solution!\n");   
        } else {
            double l = 0,   r = 500, mid, height = -1;
            for(int i=1;i<100;i++)
            {
                mid = (l+r) / 2;
                if(jug(mid))    l = mid,    height = mid;
                else    r = mid;
            }
            printf("%.4lf\n", height);
        }
    }
}