Added some experimentation code

Author: Chillee

content/geometry/HalfPlane.h

@@ -25,7 +25,9 @@ bool cmp(Line a, Line b) {
 	auto s = angDiff(a, b);
 	return s == 0 ? sideOf(sp(b), a[0]) >= 0 : s < 0;
 }
-vector<P> halfPlaneIntersection(vector<Line> vs) {
+std::random_device rd;
+std::mt19937 e2(rd());
+vector<P> halfPlaneIntersection(vector<Line> vs, const double EPS=1e-8) {
 	sort(all(vs), cmp);
 	// for (auto l: vs) {
 	// 	cout<<l[0]<<" -> "<<l[1]<<endl;
@@ -39,10 +41,12 @@ 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], 1e-8) < 0) at--;
-		while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah], 1e-8)<0) ah++;
-		auto res = lineInter(sp(vs[i]), sp(deq[at]));
+const double mult = 1.414;
+std::uniform_real_distribution<> dist(-EPS, EPS);
+while (ah<at && sideOf(sp(vs[i]),ans[at-1], mult*EPS) < 0) at--;
+while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah], mult*EPS)<0) ah++;
+auto res = lineInter(sp(vs[i]), sp(deq[at]));
+res.second = res.second + P(dist(e2), dist(e2));
 		if (res.first != 1) {
 			// cout<<"46"<<endl;
 			continue;

content/geometry/lineIntersection.h

@@ -23,20 +23,20 @@ Products of three coordinates are used in intermediate steps so watch out for ov
 
 #include "Point.h"
 
-template<class P>
-pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
-	auto d = (e1-s1).cross(e2-s2);
-	if (d == 0)  //if parallel
-		return {-((e1-s1).cross(s2-s1)==0 || s2==e2), P(0,0)};
-	else
-		return {1, s2-(e2-s2)*(e1-s1).cross(s2-s1)/d};
-}
-
 // template<class P>
 // pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
-// 	auto d = (e1 - s1).cross(e2 - s2);
-// 	if (d == 0) // if parallel
-// 		return {-(s1.cross(e1, s2) == 0), P(0, 0)};
-// 	auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
-// 	return {1, (s1 * p + e1 * q) / d};
+// 	auto d = (e1-s1).cross(e2-s2);
+// 	if (d == 0)  //if parallel
+// 		return {-((e1-s1).cross(s2-s1)==0 || s2==e2), P(0,0)};
+// 	else
+// 		return {1, s2-(e2-s2)*(e1-s1).cross(s2-s1)/d};
 // }
+
+template<class P>
+pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
+	auto d = (e1 - s1).cross(e2 - s2);
+	if (d == 0) // if parallel
+		return {-(s1.cross(e1, s2) == 0), P(0, 0)};
+	auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
+	return {1, (s1 * p + e1 * q) / d};
+}

fuzz-tests/geometry/HalfPlane.cpp

@@ -301,8 +301,9 @@ void testEmpty() {
     auto res = halfPlaneIntersection(t);
     assert(sz(res) == 0);
 }
-void testRandom() {
-    int lim = 5;
+void testRandom(const double EPS) {
+    int lim = 3;
+    double mxDiff = 0;
     for (int i = 0; i < 10000000; i++) {
         vector<Line> t;
         for (int i = 0; i < 4; i++) {
@@ -327,12 +328,14 @@ void testRandom() {
             t.push_back(cand);
         }
         addInf(t, lim);
-        auto res = halfPlaneIntersection(t);
+        auto res = halfPlaneIntersection(t, EPS);
         double area = sjtu::halfPlaneIntersection(t);
         double resArea = abs(polygonArea2(res) / 2);
         // double resArea = mit::Intersection_Area(convert(t));
         double diff = abs(area - resArea);
-        if (diff > EPS || isnan(diff)) {
+        mxDiff = max(diff, mxDiff);
+        continue;
+        if (diff > .1 || isnan(diff)) {
             cout << diff << ' ' << area << ' ' << resArea << endl;
             for (auto i : t)
                 cout << i[0] << "->" << i[1] << ' ';
@@ -344,9 +347,11 @@ void testRandom() {
             for (auto i : res)
                 cout << i << ',';
             cout << endl;
+            assert(false);
         }
-        assert(diff <= EPS);
     }
+    cout<<"eps: "<<EPS<<endl;
+    cout<<"mxDiff: "<<mxDiff<<endl;
 }
 vector<P> convexHull(vector<P> pts) {
 	if (sz(pts) <= 1) return pts;
@@ -367,8 +372,9 @@ int main() {
     // testInf();
     // testLine();
     // testPoint();
-    // testRandom();
-    // return 0;
+    for (double e : {1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2})
+        testRandom(e);
+    return 0;
     // testEmpty();
     // Case that messes with precision
     // vector<Line> cases({{P(1,0),P(0,2)},{P(0,1),P(2,0)},{P(0,0),P(1,1)},{P(2,2),P(1,1)},{P(3,3),P(-3,3)},{P(-3,3),P(-3,-3)},{P(-3,-3),P(3,-3)},{P(3,-3),P(3,3)}});