CodeArchive Geometry
From Ubcacm
Geometry
Includes/typedefs.
#include <algorithm> #include <iostream> #include <vector> #include <complex> #include <cmath> #include <cstdlib> #include <cstdio> #include <set> #include <map> using namespace std; typedef long long ll; typedef long double ld; typedef unsigned long long ull; typedef complex < ll > point; typedef complex < ld > fpoint; typedef vector < point > pol;
Integer geometry. Please punch me as hard as possible if I try to use floating point while the problem can be done in integers.
inline ll sign(const ll n){
return !n ? 0 : (n < 0 ? -1 : 1);
}
inline ll cp(const point &v1, const point &v2){
return v1.real()*v2.imag() - v1.imag()*v2.real();
}
inline ll dp(const point &v1, const point &v2){
return v1.real()*v2.real() + v1.imag()*v2.imag();
}
pair < ull, ull > get_absr(const point &a, const point &b, const point &c, const point &d){
return make_pair(abs(cp(a - c, d - c)), abs(cp(d - c, b - a)));
}
inline bool ray_crosses_point(const point &r1, const point &r2, const point &p){
return dp(r2 - r1, p - r1) >= 0 && cp(r2 - r1, p - r1) == 0; // in front and on the same line
}
// Prereqs: non-parallel, (non-zero). Excluding end points of the segment!!!
inline bool line_crosses_segment(const point &l1, const point &l2, const point &s1, const point &s2){
return sign(cp(l2 - l1, s1 - l1)) == -sign(cp(l2 - l1, s2 - l1)); // on opposite sides
}
// Prereqs: non-parallel, non-zero. Excluding end points of the segment, but including the beginning of the ray
inline bool ray_crosses_segment(const point &r1, const point &r2, const point &s1, const point &s2){
return line_crosses_segment(r1, r2, s1, s2) && sign(cp(r1 - s1, s2 - s1))*sign(cp(s2 - s1, r2 - r1)) >= 0;
}
// Assumes counter-clockwise orientation. Did ray r1 -> r2 enter/leave the polygon while going through p2 ?
inline ll in_out_point(const point &r1, const point &r2, const point &p1, const point &p2, const point &p3){
ll s1 = sign(cp(r2 - r1, p1 - r1)), s3 = sign(cp(r2 - r1, p3 - r1));
if(s1 == s3) return 0; // no change
if(s1 == 1 && s3 == -1) return 1; // entering
if(s1 == -1 && s3 == 1) return -1; // leaving
ll d1 = dp(r2 - r1, p1 - r1), d2 = dp(r2 - r1, p2 - r1), d3 = dp(r2 - r1, p3 - r1);
if(s3 == 0 && d3 > d2 && s1 > 0 || s1 == 0 && d2 < d1 && s3 < 0) return 1;
if(s3 == 0 && d3 < d2 && s1 < 0 || s1 == 0 && d2 > d1 && s3 > 0) return -1;
return 0;
}
// Assumes counter-clockwise orientation, non-parallel, intersection takes place somewhere in the middle.
inline ll in_out_segment(const point &r1, const point &r2, const point &s1, const point &s2){
return cp(r2 - r1, s1 - r1) > 0 ? 1 : -1;
}
ull mb = (1ll << 31); // mult. base, 1*2 bits to catch overflow
inline pair < ull, ull > mult(ull a, ull b){
ull low = (a & (mb - 1))*(b & (mb - 1));
ull middle = (a >> 31)*(b & (mb - 1)) + (a & (mb - 1))*(b >> 31);
low += ((middle & (mb - 1)) << 31);
return make_pair((a >> 31)*(b >> 31) + (middle >> 31) + (low >> 62), low & ((1ll << 62) - 1));
}
Floating point geometry - use it carefully
#define EPS 1e-9
const ld PI = acos((ld)-1);
inline ld cp(const point &v1, const point &v2){
return v1.real()*v2.imag() - v1.imag()*v2.real();
}
inline ld dp(const point &v1, const point &v2){
return v1.real()*v2.real() + v1.imag()*v2.imag();
}
// floating point sign ;) *TESTED ON 10902*
inline ld sign(const ld n){
return n > EPS ? 1 : (n < -EPS ? -1 : 0);
}
// includes end points *TESTED ON 10902*
bool segment_crosses_segment(const fpoint &a1, const fpoint &a2, const fpoint &b1, const fpoint &b2){
return sign(cp(a2 - a1, b1 - a1))*sign(cp(a2 - a1, b2 - a1)) <= 0 &&
sign(cp(b2 - b1, a1 - b1))*sign(cp(b2 - b1, a2 - b1)) <= 0;
}
// *TESTED ON 303*
inline point line_inter(const point &a, const point &b, const point &c, const point &d){
return a + cp(a - c, d - c)/cp(d - c, b - a) * (b - a);
}
inline ld get_r(const point &a, const point &b, const point &c, const point &d){
return cp(a - c, d - c) / cp(d - c, b - a);
}
// SIGNED line to point distance. Point on the right of vector (a -> b) will be NEGATIVE.
inline ld lp_dist(const point &a, const point &b, const point &p){
return cp(b - a, p - a) / abs(b - a);
}
// Projection of line segment from a to p to vector (a -> b), SIGNED - positive in front
inline ld proj_dist(const point &a, const point &b, const point &p){
return dp(b - a, p - a) / abs(b - a);
}
// Line segment (a, b) to point p distance.
inline ld lsp_dist(const point &a, const point &b, const point &p){
return proj_dist(a, b, p) > 0 && proj_dist(b, a, p) > 0 ?
abs(lp_dist(a, b, p)) : min(abs(a - p), abs(b - p));
}
// Closest point on line segment (a, b) to point p.
inline point cp_segment(const point &a, const point &b, const point &p){
if(proj_dist(a, b, p) > 0 && proj_dist(b, a, p) > 0)
return abs(lp_dist(a, b, p)) < EPS ? p : line_inter(a, b, p, p + perp(a - b));
return abs(a - p) < abs(b - p) ? a : b;
}
// Left to the vector (a -> b) will be cut off. Erases duplicate points.
pol cut_polygon(const pol &v, const point &a, const point &b){
pol out;
for(int j = 0, i = v.size() - 1; j < (int)v.size(); i = j++){
if( lp_dist(a, b, v[i]) < EPS ) out.push_back(v[i]);
if( lp_dist(a, b, v[i]) < -EPS && lp_dist(a, b, v[j]) > EPS ||
lp_dist(a, b, v[i]) > EPS && lp_dist(a, b, v[j]) < -EPS ){
point p = line_inter(a, b, v[i], v[j]);
if(!out.size() || abs(out.back() - p) > EPS) out.push_back(p);
}
}
while(out.size() && abs(out[0] - out.back()) < EPS) out.pop_back();
return out;
}
// Area of a polygon (convex or concave). Always non-negative.
ld area(pol v){
ld out = 0;
for(int j = 0, i = v.size() - 1; j < (int)v.size(); i = j++) out += cp(v[i], v[j]);
return abs(out)/2;
}
// Lexicographical comparison for points.
bool cmp_lex(const point &a, const point &b){
if(a.real() == b.real()) return a.imag() < b.imag();
return a.real() < b.real();
}
// Coparison by angle relative to the base (primary) and distance (secondary)
// WARNING!!! NORM is the square of length, watch for overflow!!!!!!!!!!!!!!!
point base;
bool cmp_base(const point &a, const point &b){
if(cp(a - base, b - base) == 0) return norm(a - base) < norm(b - base);
return cp(a - base, b - base) < 0;
}
// Finds CONVEX HULL of a polygon. Might not (will not) work with floating point.
pol convex_hull(pol v){
sort(v.begin(), v.end(), cmp_lex), base = v[0];
sort(v.begin() + 1, v.end(), cmp_base);
pol out;
for(int i = 0; i < (int)v.size(); i++){
while(out.size() > 1 && cp(out.back() - out[out.size() - 2], v[i] - out[out.size() - 2]) >= 0) out.pop_back();
out.push_back(v[i]);
}
return out;
}
// Point on line segment. Checks if p lies on line segment (a, b).
inline bool on_segment(const point &a, const point &b, const point &p){
return proj_dist(a, b, p) > -EPS && proj_dist(b, a, p) > -EPS && abs(lp_dist(a, b, p)) < EPS;
}
// Point on the boundry of a polygon
bool on_boundry(const pol &v, const point &p){
for(int i = v.size() - 1, j = 0; j < (int)v.size(); i = j++)
if(on_segment(v[i], v[j], p)) return true;
return false;
}
// I assume that the point is not on any of the line segments...
// the following two functions were tested on 881 and 10867 !!!!!!!!!!!!!
ld get_chains_angle(const pol &v, const point &p){
ld out = 0;
for(int i = v.size() - 1, j = 0; j < (int)v.size(); i = j++)
out += atan2(cp(v[i] - p, v[j] - p), dp(v[i] - p, v[j] - p));
return out;
}
// orientation does not matter !!! :)
// tested on 881 *!!!NO BOUNDRY CHECKING!!!!!*
bool point_in_polygon(const pol &v, const point &p){
return abs(get_chains_angle(v, p)) > PI; // will be either 2*PI or 0
}
// intersection of two (could be concave) polygons
// !!! not tested on any UVA problems !!! but looks to work
vector < pol > pol_inter(pol v1, pol v2){
vector < pair < point, point > > segm;
for(int k = 0; k < 2; k++){
for(int i = v1.size() - 1, j = 0; j < (int)v1.size(); i = j++){
pol p; p.push_back(v1[i]), p.push_back(v1[j]);
for(int ii = v2.size() - 1, jj = 0; jj < (int)v2.size(); ii = jj++)
if(abs(cp(v1[i] - v1[j], v2[ii] - v2[jj])) > EPS){
ld a = get_r(v1[i], v1[j], v2[ii], v2[jj]);
ld b = get_r(v2[ii], v2[jj], v1[i], v1[j]);
if(a > -EPS && a < 1 + EPS && b > -EPS && b < 1 + EPS) p.push_back(line_inter(v1[i], v1[j], v2[ii], v2[jj]));
}
sort(p.begin(), p.end(), cmp_lex);
for(int i = 0; i + 1 < (int)p.size(); i++){
point a = (p[i] + p[i + 1])/ld(2);
if(on_boundry(v2, a) || point_in_polygon(v2, a)) segm.push_back(make_pair(p[i], p[i + 1]));
}
}
v1.swap(v2);
}
vector < pol > out;
vector < bool > vis(segm.size(), false);
for(int i = 0; i < (int)segm.size(); i++)
if(!vis[i]){
out.resize(out.size() + 1), out.back().push_back(segm[i].first), out.back().push_back(segm[i].second);
point now = segm[i].second;
while(true){
bool done = true;
for(int j = i + 1; j < (int)segm.size(); j++)
if(!vis[j]){
if(abs(segm[j].first - now) < EPS){
done = false, now = segm[j].second, vis[j] = 1;
break;
}else if(abs(segm[j].second - now) < EPS){
done = false, now = segm[j].first, vis[j] = 1;
break;
}
}
if(done) break;
out.back().push_back(now);
}
if(abs(out.back().back() - out.back()[0]) < EPS) out.back().pop_back();
}
return out;
}
// Closest pair specific comparison - imag has priority
// Tested on 10245, 10750
struct cmp_cp_st{
bool operator()(const point &p1, const point &p2) const{
return p1.imag() == p2.imag() ? p1.real() < p2.real() : p1.imag() < p2.imag();
}
};
// Closest pair in n*log(n). Returns the distance.
// Modify apropriately if want it to work in ints.
// Tested on 10245, 10750
ld closest_pair(vector < point > v){
sort(v.begin(), v.end(), cmp_lex);
ld best = 1e99; int low = 0;
set < point, cmp_cp_st > help;
for(int i = 0; i < (int)v.size(); i++){
while(low < i && (v[i] - v[low]).real() > best)
help.erase(v[low++]);
for(__typeof(help.end()) it = help.lower_bound(v[i] - point(1e99, best));
it != help.end() && (*it - v[i]).imag() < best; it++) best = min(best, abs(*it - v[i]));
help.insert(v[i]);
}
return best;
}
// given two circles, find all tangent angles.
// returns four angles for the first circle.
inline vector < ld > get_tangents(const point &p1, const ld &r1, const point &p2, const ld &r2){
point v = p2 - p1; // vector from p1 to p2
ld d = abs(v); // distance between centers
ld a = acos((r1 - r2)/d); // outer tangent's angle relative to circle1 -> circle2
ld b = acos((r1 + r2)/d); // inner .......
ld s = atan2(v.imag(), v.real()); // basis angle relative to the circle1
vector < ld > out;
out.push_back(s + a), out.push_back(s - a);
out.push_back(s + b), out.push_back(s - b);
return out;
}
ld good_angle(ld a){
return fmod(fmod(a, 2*PI) + 2*PI, 2*PI);
}
ld ang_dist(ld a1, ld a2){
ld d = fmod(abs(a1 - a2), 2*PI);
return min(d, 2*PI - d);
}
// looks to work in N^3
vector < vector < pair < int, ld > > > build_graph(const vector < pair < point, ld > > &data){
vector < ld > angle(4*data.size()*data.size(), 0), temp;
for(int i = 0; i < (int)data.size(); i++)
for(int j = 0; j < (int)data.size(); j++)
if(i != j){
temp = get_tangents(data[i].first, data[i].second, data[j].first, data[j].second);
for(int k = 0; k < 4; k++) angle[4*i*data.size() + 4*j + k] = temp[k]; // i -> j
}
vector < point > pnts;
for(int i = 0; i < (int)angle.size(); i++){
point p = data[i/(4*data.size())].first;
ld r = data[i/(4*data.size())].second;
ld a = angle[i];
pnts.push_back(p + r*point(cos(a), sin(a)));
}
vector < vector < pair < int, ld > > > out(pnts.size());
// do the distances between points on same circle
for(int i = 0; i < (int)angle.size(); i += 4*data.size()){ // for each circle
int id = i/(4*data.size()); // circle #
vector < pair < ld, int > > help;
for(int j = 0; j < (int)data.size()*4; j++)
help.push_back(make_pair(good_angle(angle[i + j]), i + j));
sort(help.begin(), help.end());
for(int a = help.size() - 1, b = 0; b < (int)help.size(); a = b++)
out[help[a].second].push_back(make_pair(help[b].second, data[id].second*ang_dist(help[a].first, help[b].first))),
out[help[b].second].push_back(make_pair(help[a].second, data[id].second*ang_dist(help[a].first, help[b].first)));
}
// do the distances between different circles, whenever possible...
for(int i = 0; i < (int)data.size(); i++)
for(int j = i + 1; j < (int)data.size(); j++)
for(int k = 0; k < 4; k++){
int ii = 4*i*data.size() + 4*j + k;
int jj = 4*j*data.size() + 4*i + (k^1^((k&2) >> 1));
bool good = true; // check is line if far enough from each circle
for(int kk = 0; kk < (int)data.size() && good; kk++)
good &= (lsp_dist(pnts[ii], pnts[jj], data[kk].first) + EPS >= data[kk].second);
if(good){
out[ii].push_back(make_pair(jj, abs(pnts[ii] - pnts[jj])));
out[jj].push_back(make_pair(ii, abs(pnts[ii] - pnts[jj])));
}
}
return out;
}
-- Main.YuryKholondyrev - 29 Jan 2006
