2018-Sp08-lyk/j

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
#include<bits/stdc++.h>
using namespace std;

const double eps = 1E-9;
int sgn(double x){
    return x < -eps ? -1 : x > eps ? 1 : 0;
} 

struct vec{
    double x, y;

    vec(): x(0), y(0){}; 
    vec(double x, double y): x(x), y(y){};
    vec(FILE *in){fscanf(in, "%lf %lf", &x, &y);}

    vec operator + (vec v){return vec(x + v.x, y + v.y);}
    vec operator - (vec v){return vec(x - v.x, y - v.y);}
    vec operator * (double k){return vec(x * k, y * k);}
};

double cross(vec a, vec b){return a.x * b.y - a.y * b.x;}
bool line_intersection(vec p, vec q, vec a, vec b, vec &o, double &t){
    //p + t(q - p) = a + u(b - a)
    //t(q - p) X (b - a) = (a - p) X (b - a)
    double U = cross(a - p, b - a), D = cross(q - p, b - a);
    if(sgn(D) == 0) return 0;
    t = U / D;
    o = p + (q - p) * t;
    return 1;
}

typedef vector<vec> polygon;
vector<polygon> tri, stri[2];

double area(polygon p){
    int n = p.size();
    if(n < 3) return 0;
    double res = 0;
    for(int i = 1; i + 1 < n; i += 1) res += fabs(cross(p[i] - p[0], p[i + 1] - p[0])) / 2;
    return res;
}

void polygon_cut(polygon p, vec a, vec b, polygon & res){
    int n = p.size();
    for(int i = 0; i < n; i += 1){
        if(sgn(cross(p[i] - a, b - a)) <= 0) res.push_back(p[i]);
        vec c;
        double t = 0;
        if(line_intersection(p[i], p[(i + 1) % n], a, b, c, t) && t > 0 && t < 1)
            res.push_back(c);
    }
}

int main(){
    int n;
    double w, h;
    scanf("%d %lf %lf", &n, &w, &h);
    for(int i = 0; i < n; i += 1)
        tri.push_back({vec(stdin), vec(stdin), vec(stdin)});
    double xL = 0, xR = w, yD = 0, yU = h;
    for(bool dir = true; xR - xL > eps || yU - yD > eps; dir ^= 1){
//      printf("%lf %lf %lf %lf\n", xL, xR, yD, yU);
//      for(polygon p: tri){
//          for(vec v: p) printf("(%lf, %lf), ", v.x, v.y);
//          printf("\n");
//      }
        vec a, b;
        if(dir){
            a = vec((xL + xR) / 2, yD);
            b = vec((xL + xR) / 2, yU);
        }
        else{
            a = vec(xL, (yD + yU) / 2);
            b = vec(xR, (yD + yU) / 2);
        }
        double sarea[2] = {0, 0};
        for(polygon p: tri){
            polygon sp[2];
            for(int i = 0; i < 2; i += 1){
                if(i == 0) polygon_cut(p, a, b, sp[i]);
                else polygon_cut(p, b, a, sp[i]);
                if(sp[i].size()){
                    sarea[i] += area(sp[i]);
                    stri[i].push_back(sp[i]);
                }
            }
        }
        tri.clear();
        if(sgn(sarea[0] - sarea[1]) == - 1 || sgn(sarea[0] - sarea[1]) == 0 && stri[0].size() < stri[1].size()){
            tri.resize(stri[0].size());
            int n = tri.size();
            for(int i = 0; i < n; i += 1) tri[i].assign(stri[0][i].begin(), stri[0][i].end());
            if(dir) xR = (xL + xR) / 2;
            else yD = (yD + yU) / 2;
        }
        else{
            tri.resize(stri[1].size());
            int n = tri.size();
            for(int i = 0; i < n; i += 1) tri[i].assign(stri[1][i].begin(), stri[1][i].end());
            if(dir) xL = (xL + xR) / 2;
            else yU = (yD + yU) / 2;
        }
        stri[0].clear();
        stri[1].clear();
    }
    printf("%lf %lf", (xL + xR) / 2, (yD + yU) / 2);
}
}}}

#include<bits/stdc++.h>
using namespace std;
const double eps = 1E-9;
int sgn(double x){
    return x < -eps ? -1 : x > eps ? 1 : 0;
} 
struct vec{
    double x, y;
    vec(): x(0), y(0){}; 
    vec(double x, double y): x(x), y(y){};
    vec(FILE *in){fscanf(in, "%lf %lf", &x, &y);}
    vec operator + (vec v){return vec(x + v.x, y + v.y);}
    vec operator - (vec v){return vec(x - v.x, y - v.y);}
    vec operator * (double k){return vec(x * k, y * k);}
};
double cross(vec a, vec b){return a.x * b.y - a.y * b.x;}
bool line_intersection(vec p, vec q, vec a, vec b, vec &o, double &t){
    //p + t(q - p) = a + u(b - a)
    //t(q - p) X (b - a) = (a - p) X (b - a)
    double U = cross(a - p, b - a), D = cross(q - p, b - a);
    if(sgn(D) == 0) return 0;
    t = U / D;
    o = p + (q - p) * t;
    return 1;
}
typedef vector<vec> polygon;
vector<polygon> tri, stri[2];
double area(polygon p){
    int n = p.size();
    if(n < 3) return 0;
    double res = 0;
    for(int i = 1; i + 1 < n; i += 1) res += fabs(cross(p[i] - p[0], p[i + 1] - p[0])) / 2;
    return res;
}
void polygon_cut(polygon p, vec a, vec b, polygon & res){
    int n = p.size();
    for(int i = 0; i < n; i += 1){
        if(sgn(cross(p[i] - a, b - a)) <= 0) res.push_back(p[i]);
        vec c;
        double t = 0;
        if(line_intersection(p[i], p[(i + 1) % n], a, b, c, t) && t > 0 && t < 1)
            res.push_back(c);
    }
}
int main(){
    int n;
    double w, h;
    scanf("%d %lf %lf", &n, &w, &h);
    for(int i = 0; i < n; i += 1)
        tri.push_back({vec(stdin), vec(stdin), vec(stdin)});
    double xL = 0, xR = w, yD = 0, yU = h;
    for(bool dir = true; xR - xL > eps || yU - yD > eps; dir ^= 1){
//      printf("%lf %lf %lf %lf\n", xL, xR, yD, yU);
//      for(polygon p: tri){
//          for(vec v: p) printf("(%lf, %lf), ", v.x, v.y);
//          printf("\n");
//      }
        vec a, b;
        if(dir){
            a = vec((xL + xR) / 2, yD);
            b = vec((xL + xR) / 2, yU);
        }
        else{
            a = vec(xL, (yD + yU) / 2);
            b = vec(xR, (yD + yU) / 2);
        }
        double sarea[2] = {0, 0};
        for(polygon p: tri){
            polygon sp[2];
            for(int i = 0; i < 2; i += 1){
                if(i == 0) polygon_cut(p, a, b, sp[i]);
                else polygon_cut(p, b, a, sp[i]);
                if(sp[i].size()){
                    sarea[i] += area(sp[i]);
                    stri[i].push_back(sp[i]);
                }
            }
        }
        tri.clear();
        if(sgn(sarea[0] - sarea[1]) == - 1 || sgn(sarea[0] - sarea[1]) == 0 && stri[0].size() < stri[1].size()){
            tri.resize(stri[0].size());
            int n = tri.size();
            for(int i = 0; i < n; i += 1) tri[i].assign(stri[0][i].begin(), stri[0][i].end());
            if(dir) xR = (xL + xR) / 2;
            else yD = (yD + yU) / 2;
        }
        else{
            tri.resize(stri[1].size());
            int n = tri.size();
            for(int i = 0; i < n; i += 1) tri[i].assign(stri[1][i].begin(), stri[1][i].end());
            if(dir) xL = (xL + xR) / 2;
            else yU = (yD + yU) / 2;
        }
        stri[0].clear();
        stri[1].clear();
    }
    printf("%lf %lf", (xL + xR) / 2, (yD + yU) / 2);
}