Started fixing the bugs

Author: Chillee

content/geometry/HalfPlane.h

@@ -36,7 +36,13 @@ vector<P> halfPlaneIntersection(vector<Line> vs) {
 		if (angDiff(vs[i], vs[i - 1]) == 0) continue;
 		while (ah<at && sideOf(sp(vs[i]),ans[at-1]) < 0) at--;
 		while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah])<0) ah++;
-		ans[at++] = lineInter(sp(vs[i]), sp(deq[at])).second;
+		auto res = lineInter(sp(vs[i]), sp(deq[at]));
+		auto res2 = lineInter2(sp(vs[i]), sp(deq[at]));
+		// if (isnan(res2.x) || isinf(res2.x)) continue;
+		bool b = (isnan(res2.x) || isinf(res2.y)) == (res.first != 1);
+		assert(b);
+		if (res.first != 1) continue;
+		ans[at++] = res.second;
 		deq[at] = vs[i];
 	}
 	if (at - ah <= 2) return {};

content/geometry/lineIntersection.h

@@ -28,3 +28,8 @@ pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
 	else
 		return {1, s2-(e2-s2)*s1.cross(e1, s2)/d};
 }
+template<class P> P lineInter2(P p1, P p2, P q1, P q2) {
+    double a1 = q1.cross(q2, p1), a2 = -q1.cross(q2, p2);
+    auto res = (p1 * a2 + p2 * a1) / (a1 + a2);
+	return res;
+}

fuzz-tests/geometry/HalfPlane.cpp

@@ -10,10 +10,11 @@ typedef long long ll;
 typedef pair<int, int> pii;
 typedef vector<int> vi;
 
-#include "../../content/geometry/HalfPlane.h"
 #include "../../content/geometry/PolygonArea.h"
-
+typedef Point<double> P;
 ostream &operator<<(ostream &os, P p) { return cout << "(" << p.x << "," << p.y << ")"; }
+#include "../../content/geometry/HalfPlane.h"
+
 
 namespace sjtu {
 typedef double T;
@@ -203,7 +204,7 @@ void testEmpty() {
 }
 void testRandom() {
     int lim = 1e1;
-    for (int i = 0; i < 10000; i++) {
+    for (int i = 0; i < 10000000; i++) {
         vector<Line> t;
         for (int i = 0; i < 6; i++) {
             Line cand{P(0, 0), P(0, 0)};
@@ -214,9 +215,11 @@ void testRandom() {
         addInf(t, lim);
         auto res = halfPlaneIntersection(t);
         double area = sjtu::halfPlaneIntersection(t);
-        double diff = abs(2 * area - polygonArea2(res));
+        double resArea = polygonArea2(res)/2;
+        if (isnan(resArea)) resArea = 0;
+        double diff = abs(area - resArea);
         if (diff > EPS) {
-            cout << diff << ' ' << area << ' ' << endl;
+            cout << diff << ' ' << area << ' ' <<resArea<< endl;
             for (auto i : t)
                 cout << i[0] << "->" << i[1] << ' ';
             cout << endl;
@@ -237,6 +240,10 @@ int main() {
     testRandom();
     // Failure case for MIT's half plane cod
     // vector<Line> t({{P(9, 8), P(9, 1)}, {P(3, 3), P(3, 5)}, {P(7, 6), P(0, 8)}});
+    // Failure case for old code.
+    // vector<Line> t({{P(3,0),P(3,3)},{P(5,3), P(5,0)}, {P(4,-2), P(0,1)}, {P(3,-1), P(0,-1)}});
+    // addInf(t);
+    // cout << polygonArea2(halfPlaneIntersection(t)) << endl;
     // addInf(t);
     // cout << polygonArea2(halfPlaneIntersection(t)) << endl;
     // cout << MIT::Intersection_Area(convert(t)) << endl;