其他的咕咕咕了,不喜歡計算幾何 qaq
基本定義
向量
- \(\overrightarrow{PQ}=Q-P\)
叉積
二維 a.x * b.y-a.y * b.x
三維
sum += a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y
sum -= a.z * b.y * c.x + a.y * b.x * c.z + a.x * b.z * c.y
基本操作
將一個向量旋轉一定角度
Vec rotate(Vec a, LD b){
LD s = sin(b), c = cos(b);
return {a.x * c - a.y * s, a.x * s + a.y * c};
}
判斷一個點在直線的哪邊
有直線上的一點 \(P\),直線的方向向量 \(v\),想知道 \(Q\) 在直線哪邊:
利用叉積的性質,若 \(\overrightarrow{PQ}\times v >0\),則 \(Q\) 在直線逆時針方向,否則在順時針方向。
bool Left(Vec a, Line b){ return cross(a - b.p , b.v) > 0; }
判斷兩圓之間的關係
返回的是公切線個數
int Pos(Circle a, Circle b){
LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b);
if(dis > a.r + b.r) return 4;
if(dis == a.r + b.r) return 3;
if(dis > a.r - b.r) return 2;
if(dis == a.r - b.r) return 1;
return 0;
}
求兩條直線的交點
有直線 \(AB,CD\),求交點 \(E\):
首先確定是否只有一個交點,然後因為記錄的是直線上的一個點和直線的方向向量,所以只需要知道這個點與交點的距離 \(l\) ,再將這個點沿方向向量平移 \(l\) 個單位長度即可。利用正弦定理求解
由上圖可知,\(|\mathbf{a\times b|=|a||b|\sin\beta,|u\times b|=|u||b|\sin\theta}\)。
作商得:
交點即點 \(B+T\mathbf{a}\)。
Vec Inter(Line a, Line b){ return a.p + cross(b.v, (b.p - a.p)) / cross(b.v, a.v) * a.v;}
求兩圓的交點
注意要先判斷有沒有交點
pair <Vec, Vec> Inter(Circle a, Circle b){
LD x = a.r, y = b.r, z = Dis(a.O, b.O);
LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z));
Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x;
LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z));
Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y;
return {i1, i2};
}
求多邊形周長
double PloygonDis(Vec *a, int n){
double sum = 0;
lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]);
return sum;
}
求多邊形面積
double PloygonArea(Vec *a, int n){
double sum = 0;
lfor(i, 3, n) sum += cross(a[i - 1] - a[1], a[i] - a[1]);
return sum / 2;
}
極角排序
實現
- 使用
atan2(y, x)
函式,返回值的範圍是 \([-\pi,\pi]\); - 使用叉積大於 \(0\) 的性質,注意因為叉積無法判 180 度以上,所以可能要結合象限排序;
凸包
二維凸包
特殊結論
-
二維平面四個點求凸包面積 \(\rightarrow\) 任選三個點面積之和 / 2
二維平面三個點面積 \(\rightarrow\) 二個二維向量行列式值的絕對值 / 2
-
三維空間五個點求凸包體積 \(\rightarrow\) 任選四個點體積之和 / 2
三維空間四個點體積 \(\rightarrow\) 三個三維向量行列式值的絕對值 / 6
半平面交
定義
有若干條有向直線,要求保留每條直線其中一側,求最後保留的範圍。
離線 \(O(n\log n)\)
用有向直線(一個點和一個方向向量)表示半平面,以下預設半平面在有向直線的左側。
對有向直線按方向向量的極角排序,維護一個雙端佇列,儲存當前構成半平面的直線以及相鄰兩直線的交點。
每次加入一條有向直線,如果隊首 / 隊尾的交點在直線右側(用叉積判)則彈掉隊首 / 隊尾的直線。
需要注意的細節:
- 加入直線時,先彈隊尾,再彈隊首。
- 特判平行直線,在右側的要彈掉。
- 最後還要檢查隊尾交點是否在隊首直線的右側,如果是也要彈掉。
- 如果題目給出的半平面不一定有限制邊界,則應該手動加入一個 INF 邊界。
double HPI(Line *a, int n){
static deque <Vec> I; while(!I.empty()) I.pop_back();
static deque <Line> Q; while(!Q.empty()) Q.pop_back();
sort(a + 1, a + n + 1);
Q.push_back(a[1]);
lfor(i, 2, n){
while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back();
while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front();
if(a[i].at2 != Q.back().at2) Q.push_back(a[i]);
else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){
Q.back() = a[i]; if(!I.empty()) I.pop_back();
}else continue;
auto qwq = Q.rbegin();
if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back()));
}
while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back();
if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back()));
int cnt = 0; static Vec *b = new Vec[I.size() + 1];
for(auto x : I) b[++cnt] = x;
return PloygonDis(b, cnt);
}
積分
對積分的感性理解
-
有一個函式 \(f(x)\),求其在區間 \([a,b]\) 與 \(x\) 軸圍成的面積,\(x\) 軸上為正,\(x\) 軸下為負。
-
那麼 \(\int\) 類比與 \(\sum\) 符號,同時用 \(\mathrm{d}x\) 表示將 \(x\) 分成很多很多很小的份。
-
同時也不難意識到,這個空間是封閉的才可以求面積。
如何積分
大部分的函式都是無法精確積分的,於是採用一些公式來逼近。
自適應辛普森積分
公式:
- 二次函式的積分可以精確計算,辛普森積分即是一種拿二次函式來擬合的方式。
- 在二次函式的情況下,該式求出的即準確積分,推導過程。
自適應:
因為 \(f(x)\) 不是二次函式,那麼當然不能直接積分,於是就有了根據誤差調整的自適應做法。
LD Ars(LD l, LD r, LD eps, LD val){
LD mid = (l + r) / 2;
LD L = simpson(l, mid), R = simpson(mid, r);
if(fabs(L + R - val) <= eps) return L + R;
return Ars(l, mid, eps / 2, L) + Ars(mid, r, eps / 2, R);
}
- 非常直接的想法,如果不滿足 eps,就繼續細分割槽間。
- 注意因為要合併兩個區間,所以對誤差的要求在提高。
- 因為是自適應的,所以複雜度玄學。
- 對於擬合易錯的函式,可強制迭代一定層數。
積分在 OI 中的應用
一般用來求各種面積。
沒封裝的簡陋玩意
const LD Pi = acos(-1);
struct Vec{ LD x, y; };
void Out(Vec a){ cerr << a.x << ' ' << a.y << endl; }
void In(Vec &a){ scanf("%Lf%Lf", &a.x, &a.y); }
bool operator ==(Vec a, Vec b){ return a.x == b.x && a.y == b.y; }
bool operator !=(Vec a, Vec b){ return a.x != b.x || a.y != b.y; }
bool operator <(Vec a, Vec b){ return atan2(a.y, a.x) < atan2(b.y, b.x); }
LD atan2(Vec a){ return atan2(a.y, a.x); }
LD Cross(Vec a, Vec b){ return a.x * b.y - a.y * b.x; }
LD Dis(Vec a, Vec b){ return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); }
Vec operator -(Vec a, Vec b){ return {a.x - b.x, a.y - b.y}; }
Vec operator +(Vec a, Vec b){ return {a.x + b.x, a.y + b.y}; }
Vec operator /(Vec a, LD b){ return {a.x / b, a.y / b}; }
Vec operator *(Vec a, LD b){ return {a.x * b, a.y * b}; }
Vec operator *(LD b, Vec a){ return {a.x * b, a.y * b}; }
Vec rotate(Vec a, LD b){
LD s = sin(b), c = cos(b);
return {a.x * c - a.y * s, a.x * s + a.y * c};
}
struct Line{
Vec p, v; LD at2;
Line(){}
Line(Vec a, Vec b, LD c){ p = a, v = b, at2 = c; }
Line(Vec a, Vec b){ p = a, v = b - a, at2 = atan2(v.y, v.x); }
};
bool operator <(Line a, Line b){ return a.at2 < b.at2; }
bool Left(Vec a, Line b){ return Cross(a - b.p , b.v) > 0; }
Vec Inter(Line a, Line b){ return a.p + Cross(b.v, (b.p - a.p)) / Cross(b.v, a.v) * a.v;}
Vec operator +(Vec a, Line b){ return {a.x + b.v.x, a.y + b.v.y}; }
Line rotate(Line a, LD b){ return {a.p, rotate(a.v, b), a.at2}; }
Line operator *(Line a, LD b){ return (Line){a.p, (Vec){a.v.x * b, a.v.y * b}, a.at2}; }
Line operator /(Line a, LD b){ return (Line){a.p, (Vec){a.v.x / b, a.v.y / b}, a.at2}; }
struct Circle{ Vec O; LD r; };
void In(Circle &a){ scanf("%Lf", &a.r), In(a.O); }
int Pos(Circle a, Circle b){
LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b);
if(dis > a.r + b.r) return 4;
if(dis == a.r + b.r) return 3;
if(dis > a.r - b.r) return 2;
if(dis == a.r - b.r) return 1;
return 0;
}
pair <Vec, Vec> Inter(Circle a, Circle b){
LD x = a.r, y = b.r, z = Dis(a.O, b.O);
LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z));
Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x;
LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z));
Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y;
return {i1, i2};
}
double PloygonArea(Vec *a, int n){
double sum = 0;
lfor(i, 3, n) sum += Cross(a[i - 1] - a[1], a[i] - a[1]);
return sum / 2;
}
double PloygonDis(Vec *a, int n){
double sum = 0;
lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]);
return sum;
}
double HPI(Line *a, int n){
static deque <Vec> I; while(!I.empty()) I.pop_back();
static deque <Line> Q; while(!Q.empty()) Q.pop_back();
sort(a + 1, a + n + 1);
Q.push_back(a[1]);
lfor(i, 2, n){
while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back();
while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front();
if(a[i].at2 != Q.back().at2) Q.push_back(a[i]);
else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){
Q.back() = a[i]; if(!I.empty()) I.pop_back();
}else continue;
auto qwq = Q.rbegin();
if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back()));
}
while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back();
if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back()));
int cnt = 0; static Vec *b = new Vec[I.size() + 1];
for(auto x : I) b[++cnt] = x;
return PloygonDis(b, cnt);
}