计算几何
记录一些东西。
凸包
采用求一次上凸包然后求一次下凸包在拼起来实现。
\(type = 0\) 是下凸包， 反之是上凸包。
inline void graham(vector<vec2> &vec, int type) {
if(vec.size() == 0) { cerr << "Hull is empty !\n" << endl; exit(0); }
sort(vec.begin(), vec.end());
vector<vec2> rec; rec.pb(* vec.begin()); int sz = 0;
for(int i = 1; i < vec.size(); i ++) {
if(type == 0) while(sz >= 1 && le(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
else while(sz >= 1 && ge(cross(rec[sz - 1], rec[sz], vec[i]), 0)) rec.pop_back(), sz --;
rec.push_back(vec[i]); sz ++;
}
swap(vec, rec);
}
inline void graham_full(vector<vec2> &vec) {
vector<vec2> v1 = vec, v2 = vec;
graham(v1, 0); graham(v2, 1);
v1.pop_back(); for(int i = v2.size() - 1; i >= 1; i --) v1.push_back(v2[i]); swap(vec, v1);
}
凸包直径
采用旋转卡壳实现。每次就看移动指针以后面积会不会变大就好了。
inline double convDiameter(vector<vec2> vec) {
if(vec.size() == 2) { return (vec[0] - vec[1]).norm2(); }
vec.pb(vec[0]);
int j = 2, n = vec.size() - 1;
double res = 0;
for(int i = 0; i < vec.size() - 1; i ++) {
while(abs(cross(vec[i], vec[i + 1], vec[j])) < abs(cross(vec[i], vec[i + 1], vec[j + 1]))) j = (j + 1) % n;
Max(res, max((vec[i] - vec[j]).norm2(), (vec[i + 1] - vec[j]).norm2()));
}
return res;
}
直线交点
用三角形面积然后定比分点实现。
inline vec2 intersection(const line &l1, const line &l2) {
double ls = cross(l1.p1, l1.p2, l2.p1);
double rs = cross(l1.p1, l1.p2, l2.p2);
return l2.p1 + (l2.p2 - l2.p1) * ls / (ls - rs);
}
半平面交
把直线按照极角排序，平行的保留内侧在前， 注意排序判断是否在内侧的时候使用严格大于，从头开始加， 如果队头或者队尾的两条直线交点在当前直线右侧则弹出， 如果最后队尾两条直线交点在队头直线右侧则弹出队尾。
直线是逆时针存的。
inline void HPI(vector<line> &lv) {
vector<pair<line, double> > sorted(lv.size());
for(int i = 0; i < lv.size(); i ++) sorted[i].fi = lv[i], sorted[i].se = atan2(lv[i].direct().y, lv[i].direct().x);
sort(sorted.begin(), sorted.end(), [] (auto a, auto b) -> bool {
if(eq(a.se, b.se)) {
if(lt(cross(a.fi.p1, a.fi.p2, b.fi.p2), 0)) return 1;
else return 0;
}
else return a.se < b.se;
} );
for(int i = 0; i < lv.size(); i ++) lv[i] = sorted[i].fi;
deque<line> q;
q.push_back(lv[0]);
for(int i = 1; i < lv.size(); i ++) if(! parallel(lv[i], lv[i - 1])) {
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_back();
else break;
}
while(q.size() > 1) {
vec2 p = intersection(* q.begin(), * ++ q.begin());
if(lt(cross(lv[i].p1, lv[i].p2, p), 0)) q.pop_front();
else break;
}
q.push_back(lv[i]);
}
while(q.size() > 1) {
vec2 p = intersection(* --q.end(), * -- -- q.end());
if(lt(cross(q.begin() -> p1, q.begin() -> p2, p), 0)) q.pop_back();
else break;
}
lv = vector<line> (q.size());
for(int i = 0; i < q.size(); i ++) lv[i] = q[i];
}
闵可夫斯基和
要求 \(v1, v2\) 是凸包且极角排序完毕。
把所有边拿下来归并即可。注意最后弹出一个初始点。
inline vector<vec2> Minkowski(vector<vec2> v1, vector<vec2> v2) {
// v1, v2 is sorted
vector<vec2> s1(v1.size()), s2(v2.size());
for(int i = 1; i < s1.size(); i ++) s1[i - 1] = v1[i] - v1[i - 1]; s1[s1.size() - 1] = v1[0] - v1[s1.size() - 1];
for(int i = 1; i < s2.size(); i ++) s2[i - 1] = v2[i] - v2[i - 1]; s2[s2.size() - 1] = v2[0] - v2[s2.size() - 1];
vector<vec2> hull(v1.size() + v2.size() + 1);
int p1 = 0, p2 = 0, cnt = 0;
hull[cnt ++] = v1[0] + v2[0];
while(p1 < s1.size() && p2 < s2.size()) {
hull[cnt] = hull[cnt - 1] + (ge(s1[p1] ^ s2[p2], 0) ? s1[p1 ++] : s2[p2 ++]);
cnt ++;
}
while(p1 < s1.size()) hull[cnt] = hull[cnt - 1] + s1[p1 ++], cnt ++;
while(p2 < s2.size()) hull[cnt] = hull[cnt - 1] + s2[p2 ++], cnt ++;
hull.pop_back();
return hull;
}
辛普森积分
积分公式 ： \((R - L) \times (f(L) + f(R) + 4f(M)) \times \frac{1}{6}\)
db calc(db L, db R) {
db mid = (L + R) / 2;
return (f(L) + f(R) + 4 * f(mid)) / 6 * (R - L);
}
db simpson(db L, db R, int dep) {
if(abs(L - R) <= 1e-9 && dep >= 4) return 0;
db mid = (L + R) / 2;
if(abs(calc(L, R) - calc(L, mid) - calc(mid, R)) <= 1e-9 && dep >= 4) return calc(L, R);
return simpson(L, mid, dep + 1) + simpson(mid, R, dep + 1);
}
圆圆交点
把一个圆放到圆心， 另外一个旋转到水平， 解方程以后转回去。
平面图转对偶图
使用最小左转法。把所有边拆为两条有向边， 然后在每个点上极角排序， 每次一个边在终点找反向边然后看极角序在反向边前面的一条边记为 \(nxt\) ， 然后跳 \(nxt\) 编号， 根据面积正负判断是否是无界域。
void build() {
for(int i = 1; i <= n; i ++) sort(e[i].begin(), e[i].end());
for(int i = 2; i <= tot; i ++) {
int y = E[i].y;
auto it = lower_bound(e[y].begin(), e[y].end(), E[i ^ 1]);
if(it == e[y].begin()) it = -- e[y].end(); else -- it;
nxt[i] = it -> id;
}
for(int i = 2; i <= tot; i ++) if(! pos[i]) {
pos[i] = pos[nxt[i]] = ++ cnt;
int u = E[i].x;
for(int j = nxt[i]; E[j].y != E[i].x; j = nxt[j], pos[j] = cnt)
s[cnt] += ( (nod[E[j].x] - nod[u]) ^ (nod[E[j].y] - nod[u]) );
if(s[cnt] < 1e-9) rt = cnt;
}
for(int i = 2; i <= tot; i ++) {
hv[pos[i]].pb(i);
g[pos[i]].pb( {E[i].w, pos[i ^ 1]} ) ;
}
}
处理圆和多边形包含关系
使用 \(set\) 直接维护， 具体见 [SCOI2012]Blinker 的噩梦 。