Added line projection/reflection, circle line intersection, and a debugging utility (#133)

Author: Chillee

content/geometry/CircleLine.h

@@ -0,0 +1,25 @@
+/**
+ * Author: Victor Lecomte, chilli
+ * Date: 2019-10-29
+ * License: CC0
+ * Source: https://vlecomte.github.io/cp-geo.pdf
+ * Description: Finds the intersection between a circle and a line. Returns a
+ * vector of either 0, 1, or 2 intersection points. P is intended to be
+ * Point<double>
+ * Status: unit tested
+ */
+
+#pragma once
+
+#include "Point.h"
+#include "lineDistance.h"
+#include "LineProjectionReflection.h"
+
+template<class P>
+vector<P> circleLine(P c, double r, P a, P b) {
+	double h2 = r*r - a.cross(b,c)*a.cross(b,c)/(b-a).dist2();
+	if (h2 < 0) return {};
+	P p = lineProj(a, b, c), h = (b-a).unit() * sqrt(h2);
+	if (h2 == 0) return {p};
+	return {p - h, p + h};
+}

content/geometry/LineProjectionReflection.h

@@ -0,0 +1,21 @@
+/**
+ * Author: Victor Lecomte, chilli
+ * Date: 2019-10-29
+ * License: CC0
+ * Source: https://vlecomte.github.io/cp-geo.pdf
+ * Description: Projects point p onto line ab. Set refl=true to get reflection
+ * of point p across line ab insted. The wrong point will be returned if P is
+ * an integer point and the desired point doesn't have integer coordinates.
+ * Products of three coordinates are used in intermediate steps so watch out
+ * for overflow.
+ * Status: fuzz-tested
+ */
+#pragma once
+
+#include "Point.h"
+
+template<class P>
+P lineProj(P a, P b, P p, bool refl=false) {
+	P v = b - a;
+	return p - v.perp()*(1+refl)*v.cross(p-a)/v.dist2();
+}

content/geometry/Point.h

@@ -34,4 +34,6 @@ struct Point {
 	// returns point rotated 'a' radians ccw around the origin
 	P rotate(double a) const {
 		return P(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a)); }
+	friend ostream& operator<<(ostream& os, P p) {
+		return os << "(" << p.x << "," << p.y << ")"; }
 };

fuzz-tests/geometry/CircleLine.cpp

@@ -0,0 +1,34 @@
+#include <bits/stdc++.h>
+using namespace std;
+
+#define rep(i, a, b) for (int i = a; i < (b); ++i)
+#define trav(a, x) for (auto &a : x)
+#define all(x) x.begin(), x.end()
+#define sz(x) (int)(x).size()
+typedef long long ll;
+typedef pair<int, int> pii;
+typedef vector<int> vi;
+
+#include "../../content/geometry/CircleLine.h"
+
+typedef Point<double> P;
+int main() {
+    cin.sync_with_stdio(0);
+    cin.tie(0);
+    cin.exceptions(cin.failbit);
+    {
+        auto res = circleLine(P(0, 0), 1, P(-1, -1), P(1, 1));
+        assert(res.size() == 2);
+        assert((res[1]-P(sqrt(2)/2, sqrt(2)/2)).dist() < 1e-8);
+    }
+    {
+        auto res = circleLine(P(0, 0), 1, P(-5,  1), P(5, 1));
+        assert(res.size() == 1);
+        assert((res[0]-P(0,1)).dist() < 1e-8);
+    }
+    {
+        auto res = circleLine(P(4, 4), 1, P(0,  0), P(5, 0));
+        assert(res.size() == 0);
+    }
+    cout<<"Tests passed!"<<endl;
+}

fuzz-tests/geometry/LineProjectionReflection.cpp

@@ -0,0 +1,31 @@
+#include <bits/stdc++.h>
+using namespace std;
+
+#define rep(i, a, b) for (int i = a; i < (b); ++i)
+#define trav(a, x) for (auto &a : x)
+#define all(x) x.begin(), x.end()
+#define sz(x) (int)(x).size()
+typedef long long ll;
+typedef pair<int, int> pii;
+typedef vector<int> vi;
+
+#include "../../content/geometry/LineProjectionReflection.h"
+
+typedef Point<double> P;
+int main() {
+    cin.sync_with_stdio(0);
+    cin.tie(0);
+    const int lim = 5;
+    for (int i = 0; i < 100000; i++) {
+        P p = P(rand() % lim, rand() % lim);
+        P a = P(rand() % lim, rand() % lim);
+        P b = P(rand() % lim, rand() % lim);
+        while (a == b)
+            b = P(rand() % lim, rand() % lim);
+        auto proj = lineProj(a, b, p, false);
+        auto refl = lineProj(a, b, p, true);
+        assert(lineDist(a, b, proj) < 1e-8);
+        auto manProj = (refl + p) / 2;
+        assert((proj-manProj).dist() < 1e-8);
+    }
+}