merge mzuenni changes

8 files changed, 675 insertions, 13 deletions

diff --git a/test/geometry/delaunay.cpp b/test/geometry/delaunay.cpp
index 7f8ec30..5740b95 100644
--- a/test/geometry/delaunay.cpp
+++ b/test/geometry/delaunay.cpp
@@ -16,11 +16,11 @@ vector<pt> convexHull(vector<pt> pts){
 	vector<pt> h(2 * sz(pts));
 	auto half = [&](auto begin, auto end, int t) {
 		for (auto it = begin; it != end; it++) {
-			while (k > t && cross(h[k-2], h[k-1], *it) < 0) k--;//allow collinear points!
+			while (k > t && cross(h[k-2], h[k-1], *it) < 0) k--; //allow collinear points!
 			h[k++] = *it;
 	}};
-	half(all(pts), 1);// Untere Hülle.
-	half(next(pts.rbegin()), pts.rend(), k);// Obere Hülle.
+	half(all(pts), 1); // Untere Hülle.
+	half(next(pts.rbegin()), pts.rend(), k); // Obere Hülle.
 	h.resize(k);
 	return h;
 }
diff --git a/test/geometry/formulas3d.cpp b/test/geometry/formulas3d.cpp
new file mode 100644
index 0000000..f6c04b7
--- /dev/null
+++ b/test/geometry/formulas3d.cpp
@@ -0,0 +1,236 @@
+#include "../util.h"
+constexpr double EPS = 1e-6;
+struct pt3 {
+	double x, y, z;
+	pt3 operator-(const pt3 o) const {
+		return {x-o.x, y-o.y, z-o.z};
+	}
+	pt3 operator+(const pt3 o) const {
+		return {x+o.x, y+o.y, z+o.z};
+	}
+
+	pt3 operator*(double m) const {
+		return {x*m, y*m, z*m};
+	}
+
+	pt3 operator/(double d) const {
+		return {x/d, y/d, z/d};
+	}
+
+	friend ostream& operator<<(ostream& os, pt3 p) {
+		return os << "(" << p.x << ", " << p.y << ", " << p.z << ")";
+	}
+};
+
+pt3 cross(pt3 a, pt3 b);
+double abs(pt3 p);
+double distToLine(pt3 a, pt3 b, pt3 p) {
+	return abs(cross(p - a, b - a)) / abs(b - a);
+}
+
+#include <geometry/formulas3d.cpp>
+
+namespace Random {
+	pt3 point3d(ll l, ll r) {
+		return {
+			(double)integer<ll>(l, r),
+			(double)integer<ll>(l, r),
+			(double)integer<ll>(l, r),
+		};
+	}
+
+	template<size_t C>
+	std::array<pt3, C> distinct(ll l, ll r) {
+		std::array<pt3, C> res = {};
+		for (size_t i = 0; i < C; i++) {
+			bool unqiue;
+			do {
+				res[i] = point3d(l, r);
+				unqiue = true;
+				for (size_t j = 0; j < i; j++) {
+					unqiue &= res[j].x != res[i].x ||
+					          res[j].y != res[i].y ||
+					          res[j].z != res[i].z;
+				}
+			} while (!unqiue);
+		}
+		return res;
+	}
+}
+
+void test_dot(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto p = Random::point3d(-range, range);
+		auto q = Random::point3d(-range, range);
+
+		auto expected = p.x*q.x + p.y*q.y + p.z*q.z;
+		auto got = dot(p, q);
+
+		if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested dot: " << queries << endl;
+}
+
+void test_cross(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto p = Random::point3d(-range, range);
+		auto q = Random::point3d(-range, range);
+
+		auto expected = pt3{p.y*q.z - p.z*q.y,
+		                    p.z*q.x - p.x*q.z,
+		                    p.x*q.y - p.y*q.x};
+		auto got = cross(p, q);
+
+		if (got.x != expected.x) cerr << "error: x" << FAIL;
+		if (got.y != expected.y) cerr << "error: y" << FAIL;
+		if (got.z != expected.z) cerr << "error: z" << FAIL;
+		queries++;
+	}
+	cerr << "tested cross: " << queries << endl;
+}
+
+void test_abs(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto p = Random::point3d(-range, range);
+
+		auto expected = sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
+		auto got = abs(p);
+
+		if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested abs: " << queries << endl;
+}
+
+void test_mixed(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto a = Random::point3d(-range, range);
+		auto b = Random::point3d(-range, range);
+		auto c = Random::point3d(-range, range);
+
+		auto expected = dot(cross(a, b), c);
+		auto got = mixed(a, b, c);
+
+		if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested mixed: " << queries << endl;
+}
+
+void test_ccw(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto a = Random::point3d(-range, range);
+		auto b = Random::point3d(-range, range);
+		auto c = Random::point3d(-range, range);
+		auto d = Random::point3d(-range, range);
+
+		ll expected = mixed(b - a, c - a, d - a);
+		if (expected < 0) expected = -1;
+		if (expected > 0) expected = 1;
+		auto got = ccw(a, b, c, d);
+
+		if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested ccw: " << queries << endl;
+}
+
+void test_distToPlane(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [p, q] = Random::distinct<2>(-range, range);
+		auto [a, b, c] = Random::distinct<3>(-range, range);
+
+		auto norm = cross(a - c, b - c);
+		norm = norm / abs(norm);
+
+		auto gotA = distToPlane(a, b, c, p);
+		auto gotB = distToPlane(a, b, c, q);
+		auto got = ccw(a, b, c, p) * ccw(a, b, c, q) < 0 ? gotA + gotB : abs(gotA - gotB);
+
+		auto expected = abs(dot(norm, p - q));
+
+		if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested distToPlane: " << queries << endl;
+}
+
+void stress_pointOnPlane(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto p = Random::point3d(-range, range);
+		auto [a, b, c] = Random::distinct<3>(-range, range);
+
+		bool expected = ccw(a, b, c, p) == 0;
+		bool got = pointOnPlane(a, b, c, p);
+
+		if (got != expected) cerr << "error" << FAIL;
+		queries++;
+	}
+	cerr << "tested pointOnPlane: " << queries << endl;
+}
+
+void test_linePlaneIntersection(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [p, q] = Random::distinct<2>(-range, range);
+		auto [a, b, c] = Random::distinct<3>(-range, range);
+
+		if (abs(dot(p - q, cross(a-c, b-c))) < 1e-6) continue; 
+
+		auto got = linePlaneIntersection(p, q, a, b, c);
+
+		if (distToLine(p, q, got) > 1e-6) cerr << "error: 1" << FAIL;
+		if (distToPlane(a, b, c, got) > 1e-6) cerr << "error: 2" << FAIL;
+		queries++;
+	}
+	cerr << "tested linePlaneIntersection: " << queries << endl;
+}
+
+void test_lineLineDist(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [p, q] = Random::distinct<2>(-range, range);
+		auto [a, b] = Random::distinct<2>(-range, range);
+
+		double expected = 0;
+		if (ccw(a, b, p, q) != 0) {
+			pt3 c = a + p - q;
+			expected = distToPlane(a, b, c, p);
+		}
+		auto got = lineLineDist(p, q, a, b);
+		
+		if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL;
+		queries++;
+	}
+	cerr << "tested lineLineDist: " << queries << endl;
+}
+
+int main() {
+	test_dot(100);
+	test_dot(1'000'000);
+	test_cross(100);
+	test_cross(1'000'000);
+	test_abs(100);
+	test_abs(1'000'000);
+	test_mixed(100);
+	test_mixed(1'000'000);
+	test_ccw(100);
+	test_ccw(1'000'000);
+
+	test_distToPlane(100);
+	test_distToPlane(1'000'000);
+	stress_pointOnPlane(100);
+	stress_pointOnPlane(1'000'000);
+	test_linePlaneIntersection(100);
+	test_linePlaneIntersection(1'000);//can be far away
+	test_lineLineDist(100);
+	test_lineLineDist(1'000'000);
+}
diff --git a/test/geometry/hpi.cpp b/test/geometry/hpi.cpp
new file mode 100644
index 0000000..a2326bc
--- /dev/null
+++ b/test/geometry/hpi.cpp
@@ -0,0 +1,291 @@
+#include "../util.h"
+constexpr ll EPS = 0;
+#define double ll
+#define polar polar<ll>
+#include <geometry/formulas.cpp>
+#undef polar
+#undef double
+#include "../geometry.h"
+ll sgn(ll x) {
+	return (x > 0) - (x < 0);
+}
+#include <geometry/hpi.cpp>
+
+//https://cp-algorithms.com/geometry/halfplane-intersection.html
+namespace cpalgo {
+	// Redefine epsilon and infinity as necessary. Be mindful of precision errors.
+	const long double eps = 1e-9, inf = 1e9; 
+
+	// Basic point/vector struct.
+	struct Point { 
+
+		long double x, y;
+		explicit Point(long double x_ = 0, long double y_ = 0) : x(x_), y(y_) {}
+		Point(pt p) : x(real(p)), y(imag(p)) {}
+
+		// Addition, substraction, multiply by constant, dot product, cross product.
+
+		friend Point operator + (const Point& p, const Point& q) {
+			return Point(p.x + q.x, p.y + q.y); 
+		}
+
+		friend Point operator - (const Point& p, const Point& q) { 
+			return Point(p.x - q.x, p.y - q.y); 
+		}
+
+		friend Point operator * (const Point& p, const long double& k) { 
+			return Point(p.x * k, p.y * k); 
+		} 
+
+		friend long double dot(const Point& p, const Point& q) {
+			return p.x * q.x + p.y * q.y;
+		}
+
+		friend long double cross(const Point& p, const Point& q) { 
+			return p.x * q.y - p.y * q.x; 
+		}
+
+		friend std::ostream& operator<<(std::ostream& os, const Point& p) {
+			return os << "(" << p.x << ", " << p.y << ")";
+		}
+
+
+	};
+
+	// Basic half-plane struct.
+	struct Halfplane { 
+
+		// 'p' is a passing point of the line and 'pq' is the direction vector of the line.
+		Point p, pq; 
+		long double angle;
+
+		Halfplane() {}
+		Halfplane(const Point& a, const Point& b) : p(a), pq(b - a) {
+			angle = atan2l(pq.x, -pq.y);//patched to get same sort
+		}
+		Halfplane(array<pt, 2> ps) : Halfplane(ps[0], ps[1]) {}
+		Halfplane(hp h) : Halfplane(h.from, h.to) {}
+
+		// Check if point 'r' is outside this half-plane. 
+		// Every half-plane allows the region to the LEFT of its line.
+		bool out(const Point& r) {
+			return cross(pq, r - p) < -eps; 
+		}
+
+		// Comparator for sorting. 
+		bool operator < (const Halfplane& e) const { 
+			return angle < e.angle;
+		} 
+
+		// Intersection point of the lines of two half-planes. It is assumed they're never parallel.
+		friend Point inter(const Halfplane& s, const Halfplane& t) {
+			long double alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq);
+			return s.p + (s.pq * alpha);
+		}
+
+		friend std::ostream& operator<<(std::ostream& os, const Halfplane& hp) {
+			return os << hp.p << " " << hp.p+hp.pq;
+		}
+	};
+
+	// Actual algorithm
+	vector<Point> hp_intersect(vector<Halfplane>& H) { 
+
+		/*Point box[4] = {  // Bounding box in CCW order
+			Point(inf, inf), 
+			Point(-inf, inf), 
+			Point(-inf, -inf), 
+			Point(inf, -inf) 
+		};
+
+		for(int i = 0; i<4; i++) { // Add bounding box half-planes.
+			Halfplane aux(box[i], box[(i+1) % 4]);
+			H.push_back(aux);
+		}*/
+		// Add bounding box half-planes, fixed:
+		H.push_back({pt( 1e9,  1e9), pt( 1e9-1,  1e9  )});
+		H.push_back({pt(-1e9,  1e9), pt(-1e9  ,  1e9-1)});
+		H.push_back({pt(-1e9, -1e9), pt(-1e9+1, -1e9  )});
+		H.push_back({pt( 1e9, -1e9), pt( 1e9  , -1e9+1)});
+
+		// Sort by angle and start algorithm
+		sort(H.begin(), H.end());
+		deque<Halfplane> dq;
+		int len = 0;
+		for(int i = 0; i < int(H.size()); i++) {
+
+			// Remove from the back of the deque while last half-plane is redundant
+			while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) {
+				dq.pop_back();
+				--len;
+			}
+
+			// Remove from the front of the deque while first half-plane is redundant
+			while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
+				dq.pop_front();
+				--len;
+			}
+
+			// Special case check: Parallel half-planes
+			if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) {
+				// Opposite parallel half-planes that ended up checked against each other.
+				if (dot(H[i].pq, dq[len-1].pq) < 0.0)
+					return vector<Point>();
+
+				// Same direction half-plane: keep only the leftmost half-plane.
+				if (H[i].out(dq[len-1].p)) {
+					dq.pop_back();
+					--len;
+				}
+				else continue;
+			}
+
+			// Add new half-plane
+			dq.push_back(H[i]);
+			++len;
+		}
+
+		// Final cleanup: Check half-planes at the front against the back and vice-versa
+		while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) {
+			dq.pop_back();
+			--len;
+		}
+
+		while (len > 2 && dq[len-1].out(inter(dq[0], dq[1]))) {
+			dq.pop_front();
+			--len;
+		}
+
+		// Report empty intersection if necessary
+		if (len < 3) return vector<Point>();
+
+		// Reconstruct the convex polygon from the remaining half-planes.
+		vector<Point> ret(len);
+		for(int i = 0; i+1 < len; i++) {
+			ret[i] = inter(dq[i], dq[i+1]);
+		}
+		ret.back() = inter(dq[len-1], dq[0]);
+		return ret;
+	}
+}
+
+//check if a is dominated by b and c
+bool naiveCheck(cpalgo::Halfplane a, cpalgo::Halfplane b, cpalgo::Halfplane c) {
+	return a.out(inter(b, c));
+}
+
+void test_check(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto a = Random::line(range);
+		auto b = Random::line(range);
+		auto c = b;
+		while (cross(b[0] - b[1], c[0] - c[1]) == 0) c = Random::line(range);
+		
+		bool got = hp(a[0], a[1]).check(hp(b[0], b[1]), hp(c[0], c[1]));
+		bool expected = naiveCheck(a, b, c);
+
+		if (got != expected) {
+			cout << tries << endl;
+			cout << a[0] << " " << a[1] << endl;
+			cout << b[0] << " " << b[1] << endl;
+			cout << c[0] << " " << c[1] << endl;
+			cout << boolalpha << got << " " << expected << endl;
+			cerr << "error" << FAIL;
+		}
+		queries++;
+	}
+	cerr << "tested random queries: " << queries << endl;
+}
+
+void stress_test(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto centre = Random::point<ll>(-range / 2, range / 2);
+		int n = Random::integer<int>(3, 30);
+
+		vector<hp> hps1 = {//+cpalgo bounding box
+			{pt( 1e9,  1e9), pt(-1e9,  1e9)},
+			{pt(-1e9,  1e9), pt(-1e9, -1e9)},
+			{pt(-1e9, -1e9), pt( 1e9, -1e9)},
+			{pt( 1e9, -1e9), pt( 1e9,  1e9)},
+		};
+
+		vector<cpalgo::Halfplane> hps2;
+		for (int i = 0; i < n; i++) {
+			auto [a, b] = Random::line(range);
+			if (cross(a, b, centre) < 0) swap(a, b);
+			hps1.emplace_back(a, b);
+			hps2.emplace_back(a, b);
+		}
+
+		auto expected = cpalgo::hp_intersect(hps2);
+		auto gotHp = intersect(hps1);
+		if (sz(gotHp) == 1) cerr << "WHAT?" << FAIL;
+		for (hp h : gotHp) {
+			if (h.dummy()) cerr << "DUMMY!" << FAIL;//we added a bounding box...
+		}
+		vector<cpalgo::Point> got;
+		if (!gotHp.empty()) {
+			gotHp.push_back(gotHp[0]);
+			for (int i = 0; i + 1 < sz(gotHp); i++) {
+				got.push_back(inter(cpalgo::Halfplane(gotHp[i]), cpalgo::Halfplane(gotHp[i + 1])));
+			}
+		}
+
+		bool eq = sz(got) == sz(expected);
+		for (ll i = 0; eq && i < sz(got); i++) {
+			eq &= float_error(got[i].x, expected[i].x) < 1e-6;
+			eq &= float_error(got[i].y, expected[i].y) < 1e-6;
+		}
+
+		if (!eq) {
+			cerr << tries << endl;
+			cerr << setprecision(20) << endl;
+			for (auto tmp : hps2) cerr << tmp << endl;
+			cerr << endl;
+			for (auto tmp : expected) cerr << tmp << endl;
+			cerr << endl;
+			for (auto tmp : got) cerr << tmp << endl;
+			cerr << FAIL;
+		}
+		queries += n;
+	}
+	cerr << "tested random queries: " << queries << endl;
+}
+
+/*constexpr int N = 1'000'000;
+void performance_test() {
+	timer t;
+	hash_t hash = 0;
+	double maxTime = 0;
+
+	vector<pt> ps;
+	for (int i = 0; i*i <= N; i++) {
+		for (int j = 0; j*j <= N; j++) {
+			ps.emplace_back(i, j);
+		}
+	}
+	t.start();
+	hash = shortestDist(ps);
+	t.stop();
+	maxTime = max(maxTime, t.time);
+
+	ps = Random::points<ll>(N, -1'000'000'000, 1'000'000'000);
+	t.reset();
+	t.start();
+	hash += shortestDist(ps);
+	t.stop();
+	maxTime = max(maxTime, t.time);
+
+	if (maxTime > 500) cerr << "too slow: " << maxTime << FAIL;
+	cerr << "tested performance: " << maxTime << "ms (hash: " << hash << ")" << endl;
+}*/
+
+int main() {
+	test_check(10);
+	test_check(100);
+	stress_test(10);
+	stress_test(100);
+	//performance_test();
+}
diff --git a/test/geometry/lines.cpp b/test/geometry/lines.cpp
new file mode 100644
index 0000000..7b1b99a
--- /dev/null
+++ b/test/geometry/lines.cpp
@@ -0,0 +1,107 @@
+#include "../util.h"
+constexpr double EPS = 1e-6;
+#define ll double
+double gcd(double x, double /**/) {return x;} //hacky
+#include <geometry/formulas.cpp>
+#undef ll
+#include <geometry/linesAndSegments.cpp>
+#include <geometry/lines.cpp>
+
+#include "../geometry.h"
+
+void stress_pointsToLine(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [a, b] = Random::line(range);
+
+		line gotA(a, b);
+		if (real(a) != real(b)) {
+			gotA.a /= gotA.b;
+			gotA.c /= gotA.b;
+			gotA.b /= gotA.b;
+		} else {
+			gotA.c /= gotA.a;
+			gotA.a /= gotA.a;
+		}
+		line gotB = pointsToLine(a, b);
+
+		if (!same(gotA, gotB)) cerr << "error" << FAIL;
+		queries++;
+	}
+	cerr << "tested pointsToLine: " << queries << endl;
+}
+
+void stress_same(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [a, b] = Random::line(range);
+		auto [c, d] = Random::line(range);
+
+		line lAB = pointsToLine(a, b);
+		line lCD = pointsToLine(c, d);
+
+		auto got = same(lAB, lCD);
+		auto expected = pointOnLine(a, b, c) && pointOnLine(a, b, d);
+
+		if (got != expected) cerr << "error" << FAIL;
+		queries++;
+	}
+	cerr << "tested same: " << queries << endl;
+}
+
+void stress_parallel(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [a, b] = Random::line(range);
+		auto [c, d] = Random::line(range);
+
+		line lAB = pointsToLine(a, b);
+		line lCD = pointsToLine(c, d);
+
+		auto got = parallel(lAB, lCD);
+		auto expected = cross(b-a, d-c) == 0;
+
+		if (got != expected) cerr << "error" << FAIL;
+		queries++;
+	}
+	cerr << "tested parallel: " << queries << endl;
+}
+
+void stress_intersect(ll range) {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto [a, b] = Random::line(range);
+		auto [c, d] = Random::line(range);
+
+		line lAB = pointsToLine(a, b);
+		line lCD = pointsToLine(c, d);
+
+		if (same(lAB, lCD)) continue;
+
+		pt gotPT;
+		auto got = intersect(lAB, lCD, gotPT);
+		auto expected = lineIntersection(a, b, c, d);
+
+		if (got != expected) cerr << "error: 1" << FAIL;
+		if (got) {
+			pt expectedPt = lineIntersection2(a, b, c, d);
+			if (float_error(real(gotPT), real(expectedPt)) > 1e-6) cerr << "error: 2" << FAIL;
+			if (float_error(imag(gotPT), imag(expectedPt)) > 1e-6) cerr << "error: 2" << FAIL;
+		}
+		queries++;
+	}
+	cerr << "tested intersect: " << queries << endl;
+}
+
+int main() {
+	stress_pointsToLine(100);
+	stress_pointsToLine(100'000);
+	stress_same(10);
+	stress_same(100);
+	stress_same(1'000'000'000);//no precision issue since this will alwas be false...
+	stress_parallel(10);
+	stress_parallel(100);
+	stress_parallel(1'000'000'000);
+	stress_intersect(100);
+	stress_intersect(1'000'000'000);
+}
diff --git a/test/geometry/linesAndSegments.cpp b/test/geometry/linesAndSegments.cpp
index 2943a67..a2da3ba 100644
--- a/test/geometry/linesAndSegments.cpp
+++ b/test/geometry/linesAndSegments.cpp
@@ -30,7 +30,7 @@ void stress_lineIntersection(ll range) {
 		auto [c, d] = Random::line(range);
 		if (ccw(a, b, c) == 0 && ccw(a, b, d) == 0) continue;
 
-		bool expected = ccw(0, a-b, c-d) == 0;
+		bool expected = ccw(0, a-b, c-d) != 0;
 		bool got = lineIntersection(a, b, c, d);
 
 		if (got != expected) cerr << "error" << FAIL;
@@ -49,7 +49,7 @@ void stress_lineIntersection2(ll range) {
 		auto got = lineIntersection2(a, b, c, d);
 
 		if (distToLine(a, b, got) > 1e-6) cerr << "error: 1" << FAIL;
-		if (distToLine(a, b, got) > 1e-6) cerr << "error: 2" << FAIL;
+		if (distToLine(c, d, got) > 1e-6) cerr << "error: 2" << FAIL;
 		queries++;
 	}
 	cerr << "tested lineIntersection2: " << queries << endl;
@@ -172,7 +172,7 @@ void stress_segmentIntersection2(ll range) {
 		if (!got.empty() != tmp) cerr << "error: 1" << FAIL;
 		for (pt p : got) {
 			if (distToSegment(a, b, p) > 1e-6) cerr << "error: 2" << FAIL;
-			if (distToSegment(a, b, p) > 1e-6) cerr << "error: 3" << FAIL;
+			if (distToSegment(c, d, p) > 1e-6) cerr << "error: 3" << FAIL;
 		}
 		if (tmp) {
 			double gotDist = abs(got.front() - got.back());
@@ -220,7 +220,7 @@ int main() {
 	stress_lineIntersection(100);
 	stress_lineIntersection(1'000'000'000);
 	stress_lineIntersection2(100);
-	stress_lineIntersection2(1'000'000);
+	stress_lineIntersection2(1'000);//intersection can bet at N^3...
 	stress_distToLine(100);
 	stress_distToLine(1'000'000'000);
 	stress_projectToLine(100);
diff --git a/test/geometry/polygon.cpp b/test/geometry/polygon.cpp
index 1dd46ca..643ea70 100644
--- a/test/geometry/polygon.cpp
+++ b/test/geometry/polygon.cpp
@@ -10,7 +10,7 @@ double abs(pt p) {
 	return hypot(real(p), imag(p));
 }
 // Liegt p auf der Strecke a-b?
-bool pointOnLineSegment(pt a, pt b, pt p) {
+bool pointOnSegment(pt a, pt b, pt p) {
 	if (cross(a, b, p) != 0) return false;
 	auto dist = norm(a - b);
 	return norm(a - p) <= dist && norm(b - p) <= dist;
@@ -65,7 +65,7 @@ void test_windingNumber(ll range) {
 				int cur = details::lineSegmentIntersection(p, p + pt(1, 2'000'000'007), ps[j], ps[j+1]);
 				if (ptLess(ps[j], ps[j+1])) expected -= cur;
 				else expected += cur;
-				onBorder |= pointOnLineSegment(ps[j], ps[j+1], p);
+				onBorder |= pointOnSegment(ps[j], ps[j+1], p);
 			}
 			if (onBorder) continue;
 			if (area(ps) < 0) expected = -expected;
@@ -93,7 +93,7 @@ void test_inside(ll range) {
 			bool onBorder = false;
 			for (int j = 0; j < n; j++) {
 				count += details::lineSegmentIntersection(p, p + pt(1, 2'000'000'007), ps[j], ps[j+1]);
-				onBorder |= pointOnLineSegment(ps[j], ps[j+1], p);
+				onBorder |= pointOnSegment(ps[j], ps[j+1], p);
 			}
 			bool expected = (count % 2) && !onBorder;
 			bool got = inside(p, ps);
diff --git a/test/geometry/segmentIntersection.cpp b/test/geometry/segmentIntersection.cpp
index 9862be5..6d3ddd6 100644
--- a/test/geometry/segmentIntersection.cpp
+++ b/test/geometry/segmentIntersection.cpp
@@ -14,7 +14,7 @@ bool pointOnLineSegment(pt a, pt b, pt p) {
 }
 
 // Test auf Streckenschnitt zwischen a-b und c-d.
-bool lineSegmentIntersection(pt a, pt b, pt c, pt d) {
+bool segmentIntersection(pt a, pt b, pt c, pt d) {
 	if (ccw(a, b, c) == 0 && ccw(a, b, d) == 0)
 			return pointOnLineSegment(a,b,c) ||
 			       pointOnLineSegment(a,b,d) ||
@@ -42,7 +42,7 @@ vector<seg> randomSegs(int n, ll range) {
 bool naive(vector<seg>& segs) {
 	for (ll i = 0; i < sz(segs); i++) {
 		for (ll j = 0; j < i; j++) {
-			if (lineSegmentIntersection(segs[i].a, segs[i].b, segs[j].a, segs[j].b)) return true;
+			if (segmentIntersection(segs[i].a, segs[i].b, segs[j].a, segs[j].b)) return true;
 		}
 	}
 	return false;
@@ -60,7 +60,7 @@ void stress_test(ll range) {
 		if (got != (b >= 0)) cerr << "error: invalid ans" << FAIL;
 		auto expected = naive(segs);
 		if (got != expected) cerr << "error: intersection not found" << FAIL;
-		if (got && !lineSegmentIntersection(segs[a].a, segs[a].b, segs[b].a, segs[b].b)) cerr << "error: no intersection" << FAIL;
+		if (got && !segmentIntersection(segs[a].a, segs[a].b, segs[b].a, segs[b].b)) cerr << "error: no intersection" << FAIL;
 		queries += n;
 		intersection += got;
 		notIntersection += !got;
diff --git a/test/geometry/spheres.cpp b/test/geometry/spheres.cpp
new file mode 100644
index 0000000..16e0ebd
--- /dev/null
+++ b/test/geometry/spheres.cpp
@@ -0,0 +1,28 @@
+#include "../util.h"
+constexpr double PI = acos(-1.0);
+#pragma GCC diagnostic ignored "-Wshadow"
+#include <geometry/spheres.cpp>
+
+void test_consistent() {
+	ll queries = 0;
+	for (int tries = 0; tries < 100'000; tries++) {
+		auto pLat = Random::real<double>(-180, 180);
+		auto pLon = Random::real<double>(0, 360);
+		auto qLat = Random::real<double>(-180, 180);
+		auto qLon = Random::real<double>(0, 360);
+
+		point p(pLat, pLon);
+		point q(qLat, qLon);
+
+		auto gotA = gcDist(pLat, pLon, qLat, qLon, 1);
+		auto gotB = gcDist(p, q);
+
+		if (abs(gotA - gotB) > 1e-6) cerr << "gotA: " << gotA << ", gotB: " << gotB << FAIL;
+		queries++;
+	}
+	cerr << "tested random: " << queries << endl;
+}
+
+int main() {
+	test_consistent();
+}