SOJ 1059. Exocenter of a Trian

题目大意：给定一个三角形ABC，求三角形的Exocenter。Exocenter定义：对三角形ABC三边作正方形，得到正方形ADEB,BJHC,CFGA,连接DG,EJ,FH，取DG、EJ、FH的中点L、M、N，直线LA、BM、CN的交点即为Exocenter.

解题思路：求出直线LA和MB，用高斯消元法求交点O。先用垂直条件、等长条件和钝角条件求D、G、E、J的坐标，再求L、M的坐标。

注意：浮点误差，两垂直直线和三种情况。

代码如下：

#include <cstdio>
#include <cmath>
#include <utility>
#include <climits>
#include <vector>
using namespace std;

//typedef pair<double, double> Cor;

struct Cor{
    double x;
    double y;
    Cor(double x_ = 0, double y_ = 0) : x(x_), y(y_) {}
    Cor operator-(const Cor &other) const {
        return Cor(x - other.x, y - other.y);
    }

    Cor operator+(const Cor &other) const {
        return Cor(x + other.x, y + other.y);
    }
    Cor operator/(const double &num) const {
        return Cor(x / num, y / num);
    }
    static Cor midCor(const Cor &cor1, const Cor &cor2) {
        return (cor1 + cor2) / 2;
    }
};

const double error_range = 0.0000001;
bool doubleEqual(const double & a, const double &b) {
    return fabs(a - b) <= error_range;
}

bool doubleLarger(const double &a, const double &b) {
    return a - b > error_range;
}

bool doubleSmaller(const double &a, const double &b) {
    return b - a > error_range;
}

double dot(const Cor &A, const Cor &B) {
    return A.x * B.x + A.y * B.y;
}

Cor getExtrianglePoint(const Cor &B, const Cor &A, const Cor &C) {
    Cor one, two;
    // 如果原直线平行于x轴
    if (A.y == B.y) {
        one.x = A.x;
        one.y = A.y + fabs(B.x - A.x);
        two.x = A.x;
        two.y = A.y - fabs(B.x - A.x);
    } else if (A.x == B.x) {
        one.y = A.y;
        one.x = A.x + fabs(B.y - A.y);
        two.y = A.y;
        two.x = A.x - fabs(B.y - A.y);
    } else {
        one.x = A.x + fabs(B.y - A.y);
        one.y = A.y - (B.x - A.x) * (one.x - A.x) / (B.y - A.y);
        two.x = A.x - fabs(B.y - A.y);
        two.y = A.y - (B.x - A.x) * (two.x - A.x) / (B.y - A.y);
    }

    if (doubleSmaller(dot(one - A, C - A), 0 )) {
        return one;
    } else if (doubleSmaller(dot(two - A, C - A), 0)) {
        return two;
    } else {
        printf("error
");
        return Cor();
    }

}

vector<double> Gauss(double (*matrix)[3], const int &n) {
    vector<int> select(n - 1);
    vector<bool> visited(n - 1, false);
    bool notSingle = false;
    for (int i = 0; i < n - 1; ++i) {
        int maximum_num;
        double maximum = 0;
        for (int j = 0; j < n - 1; ++j) {
            if (!visited[j] && doubleLarger(fabs(matrix[j][i]), fabs(maximum))) {
                maximum_num = j;
                maximum = matrix[j][i];
            }
        }
        if (doubleEqual(maximum, 0)) {
            notSingle = true;
            break;
        }
        visited[maximum_num] = true;
        select[i] = maximum_num;
        for (int j = 0; j < n - 1; ++j) {
            if (j == maximum_num) continue;
            if (doubleEqual(matrix[j][i], 0)) continue;
            double coefficient = -1 * matrix[j][i] / maximum;
            for (int k = i; k < n; ++k) {
                matrix[j][k] = matrix[maximum_num][k] * coefficient + matrix[j][k];
            }
        }
    }
    vector<double> ans;
    if (notSingle) return ans;
    for (int i = 0; i < select.size(); ++i) {
        ans.push_back(matrix[select[i]][n - 1] / matrix[select[i]][i]);
    }
    return ans;
}

double GM[2][3];



int main() {
    int n;
    Cor A, B, C;
    scanf("%d", &n);
    while (n--) {
        scanf("%lf%lf%lf%lf%lf%lf", &A.x, &A.y, &B.x, &B.y, &C.x, &C.y);
        Cor D = getExtrianglePoint(B, A, C);
        Cor G = getExtrianglePoint(C, A, B);
        Cor L = (D + G) / 2;
        Cor E = getExtrianglePoint(A, B, C);
        Cor J = getExtrianglePoint(C, B, A);
        Cor M = (E + J) / 2;
        
        if (doubleEqual(L.x, A.x)) {
            GM[0][0] = 1;
            GM[0][1] = 0;
            GM[0][2] = L.x;
        } else {
            GM[0][0] = (L.y - A.y) / (L.x - A.x);
            GM[0][1] = -1;
            GM[0][2] = (L.y - A.y) / (L.x - A.x) * A.x - A.y;
        }

        if (doubleEqual(M.x, B.x)) {
            GM[1][0] = 1;
            GM[1][1] = 0;
            GM[1][2] = M.x;
        } else {
            GM[1][0] = (M.y - B.y) / (M.x - B.x);
            GM[1][1] = -1;
            GM[1][2] = (M.y - B.y) / (M.x - B.x) * B.x - B.y;
        }

        vector<double> ans = Gauss(GM, 3);
        for (int i = 0; i < ans.size(); ++i) {
            printf("%.4lf%c", doubleEqual(ans[i], 0) ? 0 : ans[i], i == ans.size() - 1 ? '
' : ' ');
        }
    }
    return 0;
}