显然$S=frac{1}{2}ah$是不可用的。(压根感觉不可优化)
考虑向量的做法:$S = frac{1}{2}(A imes B + B imes C + C imes A)$。(相当于把一个三角形拆成了三个以原点作为其中一个顶点的"有向"三角形)
于是考虑利用向量叉积对向量加法的分配律进行优化。
枚举第一条直线,剩下的直线按照极角序加入,不断计算交点和有向面积。
对于直线的方向,我是根据原点在直线的哪一侧决定的。(比如定向后,原点在它左侧)
然后画三条直线相交,讨论原点在哪,然后再讨论怎么计算有向面积。
画一张图仅供参考,1和2表示是交点被计算的顺序。
开心地发现原点在三角形内部的时候计算叉积的时候需要取相反数计入答案。
这很烦。所以取一个超级远的地方的点作为原点就可以成功避开了这个问题。
Code
1 /** 2 * Codeforces 3 * Problem#528E 4 * Accepted 5 * Time: 78ms 6 * Memory: 100k 7 */ 8 #include <bits/stdc++.h> 9 using namespace std; 10 typedef bool boolean; 11 12 //Comparison of floating point constants 13 const double eps = 1e-10; 14 //π 15 const double pi = acos((double)-1); 16 17 //Define the points and the vectors 18 typedef class Point { 19 public: 20 double x; 21 double y; 22 Point(const double x = 0.0, const double y = 0.0):x(x), y(y) { } 23 }Point, Vector; 24 25 const Point O(1e7, 1e7); 26 27 Vector operator + (Vector a, Vector b) { 28 return Vector(a.x + b.x, a.y + b.y); 29 } 30 31 Vector operator - (Vector a, Vector b) { 32 return Vector(a.x - b.x, a.y - b.y); 33 } 34 35 Vector operator * (Vector a, double b) { 36 return Vector(a.x * b, a.y * b); 37 } 38 39 Vector operator * (double b, Vector a) { 40 return Vector(a.x * b, a.y * b); 41 } 42 43 Vector operator / (Vector a, double b) { 44 return Vector(a.x / b, a.y / b); 45 } 46 47 Vector operator - (Vector a) { 48 return Vector(-a.x, -a.y); 49 } 50 51 int dcmp(double x) { 52 if(fabs(x) < eps) return 0; 53 return (x < 0) ? (-1) : (1); 54 } 55 56 double Dot(Vector a, Vector b) { 57 return a.x * b.x + a.y * b.y; 58 } 59 60 double Cross(Vector a, Vector b) { 61 return a.x * b.y - a.y * b.x; 62 } 63 64 double Area(Point a, Point b, Point c) { 65 return fabs(Cross(b - a, c - a) / 2); 66 } 67 68 Point getLineIntersection(Point A, Vector v, Point B, Vector u) { 69 Vector w = B - A; 70 double t = (Cross(w, u) / Cross(v, u)); 71 return A + t * v; 72 } 73 74 typedef class Line { 75 public: 76 Point p; 77 Vector v; 78 double ang; 79 int id; 80 81 Line() { } 82 Line(int a, int b, int c, int id):id(id) { 83 if (!b) { 84 p = Point(c * 1.0 / a, 0); 85 v = Point(0, 1); 86 } else { 87 p = Point(0, c * 1.0 / b); 88 v = Point(-b, a); 89 } 90 if (Cross(O - p, v) > 0) 91 v = -v; 92 ang = atan2(v.y, v.x); 93 } 94 95 boolean operator < (Line b) const { 96 return ang < b.ang; 97 } 98 }Line; 99 100 ostream& operator << (ostream& os, Point p) { 101 os << "(" << p.x << " " << p.y << ")"; 102 return os; 103 } 104 105 int n; 106 Line *ls; 107 108 inline void init() { 109 scanf("%d", &n); 110 ls = new Line[(n + 1)]; 111 for (int i = 1, a, b, c; i <= n; i++) { 112 scanf("%d%d%d", &a, &b, &c); 113 ls[i] = Line(a, b, c, i); 114 } 115 } 116 117 double res = 0.0; 118 inline void solve() { 119 sort(ls + 1, ls + n + 1); 120 // for (int i = 1; i <= n; i++) 121 // cerr << ls[i].p << " " << ls[i].v << " " << ls[i].ang << endl; 122 for (int i = 1; i <= n; i++) { 123 Point sP(0, 0), P; 124 for (int j = i + 1; j <= n; j++) { 125 P = getLineIntersection(ls[i].p, ls[i].v, ls[j].p, ls[j].v) - O; 126 // int d = dcmp(Cross(ls[i].v, ls[j].v)); 127 res += Cross(sP, P); 128 sP = sP + P; 129 } 130 for (int j = 1; j < i; j++) { 131 P = getLineIntersection(ls[i].p, ls[i].v, ls[j].p, ls[j].v) - O; 132 // int d = dcmp(Cross(ls[i].v, ls[j].v)); 133 res += Cross(sP, P); 134 sP = sP + P; 135 } 136 } 137 printf("%.9lf", res * 3 / n / (n - 1) / (n - 2)); 138 } 139 140 int main() { 141 init(); 142 solve(); 143 return 0; 144 }