[Counting Sort, Radix Sort] 計數排序, 基數排序

 Counting Sort是一種效率很高的排序方式，複雜度為$\ord{n+k}$，其中$k$是Bucket的大小，由此可知僅限於整數且數字範圍不能太大。根據觀察在很多應用中會有對物件以編號大小進行排序的行為，在這方面應該能做到很大的加速。

另外一個問題是Counting Sort雖然簡單，很多人甚至可以自己想到實作方法，但這也導致了標準的作法常常被人忽略。因此這裡就來給大家展示標準的Counting sort：

	#include <algorithm>
	#include <numeric>
	template <typename InputIter, typename OutputIter, typename BucketIter,
	typename KeyFunction>
	void counting_sort(InputIter First, InputIter Last, BucketIter BucketFirst,
	BucketIter BucketLast, OutputIter OutputFirst,
	KeyFunction Key) {
	std::fill(BucketFirst, BucketLast, 0);
	for (auto Iter = First; Iter != Last; ++Iter)
	++(*(BucketFirst + Key(*Iter)));
	std::partial_sum(BucketFirst, BucketLast, BucketFirst);
	do {
	--Last;
	*(OutputFirst + --(*(BucketFirst + Key(*Last)))) = std::move(*Last);
	} while (Last != First);
	}

view raw CountingSort.hpp hosted with ❤ by GitHub

參數的解釋如下：

• First, Last:
  和std::sort的定義一樣，需要排序的範圍，注意不一定要是random access iterator。
• BucketFirst, BucketLast:
  Counting Sort正統的實作方式會有一個正整數陣列作為Bucket，考量到各種應用所以這裡接傳Bucket的範圍進來能做的優化會比較多，必須要是random access iterator。
• OutputFirst:
  Counting Sort的output是直接將input存到另一個陣列中，因此OutputFirst指的是Output陣列的開頭，必須要是random access iteator，且要注意output的空間是足夠的。這邊將input存進output時是用std::move的方式，如果想要保留原本input的話可以將其拿掉。
• Key:
  這是一個函數，Key(x)必須要回傳一個0~(BucketLast-BucketFirst-1)的正整數作為被排序物件x的關鍵值。

有了Counting sort，Radix Sort就比較好解釋了。首先正統的Counting sort是stable sort，所以Key值相同的東西排序前後的先後順序是不變的。因此可以透過多次的Counting Sort來完成一些原本Counting Sort無法完成的事情。

以整數(int, -2147483648~2147483647)排序為例，可以先針對第十億位做為Key進行排序，接著再對第一億位做為Key、第一千萬位做為Key...直到十位數、個位數作為Key，最後再以正負號最為Key進行排序，這樣就可以完成一般範圍的整數排序。

實際上一般不會這樣用，通常是用在有多個Key值的情況，以下面的程式碼來說，可以自行執行看看花費的時間有多少：

	#include "CountingSort.hpp"
	#include <cassert>
	#include <chrono>
	#include <iostream>
	#include <random>
	#include <vector>
	using namespace std;
	int main() {
	vector<pair<int, int>> V1;
	mt19937 MT(7122);
	const int N = 10000000, A = 1000, B = 10000;
	for (int i = 0; i < N; ++i) {
	V1.emplace_back(MT() % A, MT() % B);
	}
	decltype(V1) V2 = V1;
	auto Begin = chrono::high_resolution_clock::now();
	sort(V2.begin(), V2.end());
	auto End = chrono::high_resolution_clock::now();
	cout << "std::sort: "
	<< chrono::duration_cast<chrono::milliseconds>(End - Begin).count()
	<< "ms\n";
	Begin = chrono::high_resolution_clock::now();
	decltype(V2) Tmp(V1.size());
	vector<int> Bucket(max(A, B), 0);
	counting_sort(V1.begin(), V1.end(), Bucket.begin(), Bucket.begin() + B,
	Tmp.begin(), [&](const auto &P) -> int { return P.second; });
	counting_sort(Tmp.begin(), Tmp.end(), Bucket.begin(), Bucket.begin() + A,
	V1.begin(), [&](const auto &P) -> int { return P.first; });
	End = chrono::high_resolution_clock::now();
	cout << "radix_sort: "
	<< chrono::duration_cast<chrono::milliseconds>(End - Begin).count()
	<< "ms\n";
	assert(V1 == V2);
	return 0;
	}

view raw main.cpp hosted with ❤ by GitHub

[Discretize Relabeling] 離散化器

離散化是演算法競賽常用的操作，在各種實際問題上也能看到其應用。最基本的情況，是對於n個可排序的元素，製造一個map使得它們可以和自己的名次一一對應，但通常的應用中這n個元素確定之後就不太會有增減的動作，因此可以存到vector中排序去除重複的部分，搜索的部分就用二分搜尋來取代。

	#include <algorithm>
	#include <stdexcept>
	#include <vector>
	template <typename T, typename Alloc = std::allocator<T>>
	class Discretizer : private std::vector<T, Alloc> {
	void build() {
	std::sort(std::vector<T, Alloc>::begin(), std::vector<T, Alloc>::end());
	std::vector<T, Alloc>::erase(std::unique(std::vector<T, Alloc>::begin(),
	std::vector<T, Alloc>::end()),
	std::vector<T, Alloc>::end());
	}
	public:
	using const_iterator = typename std::vector<T, Alloc>::const_iterator;
	using const_reverse_iterator =
	typename std::vector<T, Alloc>::const_reverse_iterator;
	using const_reference = typename std::vector<T, Alloc>::const_reference;
	using size_type = typename std::vector<T, Alloc>::size_type;
	Discretizer() = default;
	template <class InputIterator>
	Discretizer(InputIterator first, InputIterator last,
	const Alloc &alloc = Alloc())
	: std::vector<T, Alloc>(first, last, alloc) {
	build();
	}
	Discretizer(const typename std::vector<T, Alloc> &x)
	: std::vector<T, Alloc>(x) {
	build();
	}
	Discretizer(const typename std::vector<T, Alloc> &x, const Alloc &alloc)
	: std::vector<T, Alloc>(x, alloc) {
	build();
	}
	Discretizer(typename std::vector<T, Alloc> &&x)
	: std::vector<T, Alloc>(std::move(x)) {
	build();
	}
	Discretizer(typename std::vector<T, Alloc> &&x, const Alloc &alloc)
	: std::vector<T, Alloc>(std::move(x), alloc) {
	build();
	}
	Discretizer(std::initializer_list<T> il, const Alloc &alloc = Alloc())
	: std::vector<T, Alloc>(il, alloc) {
	build();
	}
	Discretizer(const Discretizer &x) : std::vector<T, Alloc>(x) {}
	Discretizer(Discretizer &&x) : std::vector<T, Alloc>(std::move(x)) {}
	Discretizer(Discretizer &&x, const Alloc &alloc)
	: std::vector<T, Alloc>(std::move(x), alloc) {}
	Discretizer &operator=(const Discretizer &x) {
	std::vector<T, Alloc>::operator=(x);
	return *this;
	}
	Discretizer &operator=(Discretizer &&x) {
	std::vector<T, Alloc>::operator=(std::move(x));
	return *this;
	}
	const_iterator begin() const noexcept {
	return std::vector<T, Alloc>::begin();
	}
	const_iterator end() const noexcept { return std::vector<T, Alloc>::end(); }
	const_reverse_iterator rbegin() const noexcept {
	return std::vector<T, Alloc>::rbegin();
	}
	const_reverse_iterator rend() const noexcept {
	return std::vector<T, Alloc>::rend();
	}
	size_type size() const noexcept { return std::vector<T, Alloc>::size(); }
	size_type capacity() const noexcept {
	return std::vector<T, Alloc>::capacity();
	}
	bool empty() const noexcept { return std::vector<T, Alloc>::empty(); }
	void shrink_to_fit() { std::vector<T, Alloc>::shrink_to_fit(); }
	const_reference operator[](size_type n) const {
	return std::vector<T, Alloc>::operator[](n);
	}
	const_reference at(size_type n) const { return std::vector<T, Alloc>::at(n); }
	const_reference front() const { return std::vector<T, Alloc>::front(); }
	const_reference back() const { return std::vector<T, Alloc>::back(); }
	void pop_back() { std::vector<T, Alloc>::pop_back(); }
	void clear() noexcept { std::vector<T, Alloc>::clear(); }
	void swap(Discretizer<T, Alloc> &x) { std::vector<T, Alloc>::swap(x); }
	const_iterator lower_bound(const_reference x) const {
	return std::lower_bound(begin(), end(), x);
	}
	const_iterator upper_bound(const_reference x) const {
	return std::upper_bound(begin(), end(), x);
	}
	size_type getIndex(const_reference x) const {
	auto Iter = lower_bound(x);
	if (Iter == end() || *Iter != x)
	throw std::out_of_range("value not exist.");
	return Iter - begin();
	}
	};

view raw Discretizer.h hosted with ❤ by GitHub

	#include <iostream>
	#include "Discretizer.h"
	int main() {
	std::vector<int> V{3, 1, 9, 2, 2, 7, 3, 100};
	Discretizer<int> D(std::move(V));
	for (auto e : D)
	std::cout << e << ' ';
	std::cout << '\n';
	std::cout << D.getIndex(7) << '\n';
	std::cout << D.getIndex(100) << '\n';
	std::cout << D.upper_bound(9) - D.lower_bound(2) << '\n';
	return 0;
	}

view raw main.cpp hosted with ❤ by GitHub

[Multiple line segment intersection] Bentley–Ottmann 演算法

 基本上這個問題就是給你一些線段(格式通常為兩個端點)，你要找出這些線段的交點。直觀的做法兩兩進行計算會花上$\ord{n^2}$的時間，但大多數的情況下交點不會很多。為了解決這個問題，修改自Shamos–Hoey演算法的Bentley–Ottmann演算法可以在$\ord{(n+k)\log n}$的時間內找出所有交點，其中$k$是交點數量。

這裡附上實作時需要用到的基本資料結構：

	template <typename T> struct point {
	T x, y;
	point() {}
	point(const T &x, const T &y) : x(x), y(y) {}
	point operator+(const point &b) const { return point(x + b.x, y + b.y); }
	point operator-(const point &b) const { return point(x - b.x, y - b.y); }
	point operator*(const T &b) const { return point(x * b, y * b); }
	bool operator==(const point &b) const { return x == b.x && y == b.y; }
	T dot(const point &b) const { return x * b.x + y * b.y; }
	T cross(const point &b) const { return x * b.y - y * b.x; }
	};
	template <typename T> struct line {
	line() {}
	point<T> p1, p2;
	T a, b, c; // ax+by+c=0
	line(const point<T> &p1, const point<T> &p2) : p1(p1), p2(p2) {}
	void pton() {
	a = p1.y - p2.y;
	b = p2.x - p1.x;
	c = -a * p1.x - b * p1.y;
	}
	T ori(const point<T> &p) const { return (p2 - p1).cross(p - p1); }
	T btw(const point<T> &p) const { return (p1 - p).dot(p2 - p); }
	int seg_intersect(const line &l) const {
	// -1: Infinitude of intersections
	// 0: No intersection
	// 1: One intersection
	// 2: Collinear and intersect at p1
	// 3: Collinear and intersect at p2
	T c1 = ori(l.p1), c2 = ori(l.p2);
	T c3 = l.ori(p1), c4 = l.ori(p2);
	if (c1 == 0 && c2 == 0) {
	bool b1 = btw(l.p1) >= 0, b2 = btw(l.p2) >= 0;
	T a3 = l.btw(p1), a4 = l.btw(p2);
	if (b1 && b2 && a3 == 0 && a4 >= 0)
	return 2;
	if (b1 && b2 && a3 >= 0 && a4 == 0)
	return 3;
	if (b1 && b2 && a3 >= 0 && a4 >= 0)
	return 0;
	return -1;
	} else if (c1 * c2 <= 0 && c3 * c4 <= 0)
	return 1;
	return 0;
	}
	point<T> line_intersection(const line &l) const {
	point<T> a = p2 - p1, b = l.p2 - l.p1, s = l.p1 - p1;
	// if(a.cross(b)==0) return INF;
	return p1 + a * (s.cross(b) / a.cross(b));
	}
	};
	template <typename T> using segment = line<T>;

view raw Basic.cpp hosted with ❤ by GitHub

演算法使用掃描線進行。掃描線是一條垂直線從左邊掃到右邊(有些實作是水平線從上面掃到下面)，並且在遇到事件點的時候進行相關處理。

線段的兩端點以及交點都作為事件點被紀錄在最終結果中。對於每個事件點$P$，我們會計算三個集合：

• U集合：所有以$P$為起始點的線段集合
• C集合：所有包含$P$的線段集合
• L集合：所有以$P$為結束點的線段集合

當然要先保證每條線段的起始點移動會在結束點的左方，只要得到線段後稍微判斷一下就可以做到了。每個事件點找出這三個集合後就可以很容易的判斷相交資訊，但要注意的是會有以下的退化情形：

• 線段退化成點：這種情況該點的U和L都會包含該線段。
• 兩線段重合：只有重合處的兩端點會被紀錄為事件點，可以根據UCL判斷出是否線段重合
• 垂直線段：排序點和線段時如果x一樣就按照y來比較
  

最後是掃描線的資料結構，需要一棵平衡的BST根據當前掃描線和各個線段切點的y值進行排序，但這件事是可以用STL做到的！我們把當前事件點傳進比較函數裡面進行計算，因為在任何一個時刻BST中的資料都是根據當前的比較函數由小排到大的，應該不算undefined behavior。另外該演算法的浮點數誤差很大，建議使用時套上處理誤差的模板或是直接用分數計算：

	#include <algorithm>
	#include <map>
	#include <set>
	#include <vector>
	template <typename T> struct SegmentIntersection {
	struct PointInfo {
	std::vector<const segment<T> *> U, C, L;
	// U: the set of segments start at the point.
	// C: the set of segments contain the point.
	// L: the set of segments end at the point.
	};
	private:
	struct PointCmp {
	bool operator()(const point<T> &a, const point<T> &b) const {
	return (a.x < b.x) || (a.x == b.x && a.y < b.y);
	}
	};
	struct LineCmp {
	const point<T> &P;
	LineCmp(const point<T> &P) : P(P) {}
	T getY(const segment<T> *s) const {
	if (s->b == 0)
	return P.y;
	return (s->a * P.x + s->c) / -s->b;
	}
	bool operator()(const line<T> *a, const line<T> *b) const {
	return getY(a) < getY(b);
	}
	};
	const std::vector<segment<T>> Segs;
	std::map<point<T>, PointInfo, PointCmp> PointInfoMap;
	std::set<point<T>, PointCmp> Queue;
	point<T> Current;
	std::multiset<const segment<T> *, LineCmp> BST;
	std::vector<segment<T>> initSegs(std::vector<segment<T>> &&Segs) {
	for (auto &S : Segs) {
	if (!PointCmp()(S.p1, S.p2))
	std::swap(S.p1, S.p2);
	S.pton();
	}
	return Segs;
	}
	void init() {
	for (auto &S : Segs) {
	PointInfoMap[S.p1].U.emplace_back(&S);
	PointInfoMap[S.p2].L.emplace_back(&S);
	Queue.emplace(S.p1);
	Queue.emplace(S.p2);
	}
	}
	void FindNewEvent(const segment<T> *A, const segment<T> *B) {
	auto Type = A->seg_intersect(*B);
	if (Type <= 0)
	return;
	point<T> P;
	if (Type == 2)
	P = A->p1;
	else if (Type == 3)
	P = A->p2;
	else
	P = A->line_intersection(*B);
	if (PointCmp()(Current, P))
	Queue.emplace(P);
	}
	void HandleEventPoint() {
	auto &Info = PointInfoMap[Current];
	segment<T> Tmp(Current, Current);
	auto LBound = BST.lower_bound(&Tmp);
	auto UBound = BST.upper_bound(&Tmp);
	std::copy_if(
	LBound, UBound, std::back_inserter(Info.C),
	[&](const segment<T> *S) -> bool { return !(S->p2 == Current); });
	BST.erase(LBound, UBound);
	auto UC = Info.U;
	UC.insert(UC.end(), Info.C.begin(), Info.C.end());
	UC.erase(std::remove_if(
	UC.begin(), UC.end(),
	[&](const segment<T> *S) -> bool { return S->p1 == S->p2; }),
	UC.end());
	std::sort(UC.begin(), UC.end(),
	[&](const segment<T> *A, const segment<T> *B) -> bool {
	return (A->p2 - Current).cross(B->p2 - Current) > 0;
	});
	if (UC.empty()) {
	if (UBound != BST.end() && UBound != BST.begin()) {
	auto Sr = *UBound;
	auto Sl = *(--UBound);
	FindNewEvent(Sl, Sr);
	}
	} else {
	if (UBound != BST.end()) {
	auto Sr = *UBound;
	FindNewEvent(UC.back(), Sr);
	}
	if (UBound != BST.begin()) {
	auto Sl = *(--UBound);
	FindNewEvent(Sl, UC.front());
	}
	}
	for (auto S : UC)
	BST.emplace(S);
	}
	public:
	SegmentIntersection(std::vector<segment<T>> Segs)
	: Segs(initSegs(std::move(Segs))), BST(Current) {
	init();
	while (Queue.size()) {
	auto It = Queue.begin();
	Current = *It;
	Queue.erase(It);
	HandleEventPoint();
	}
	}
	const std::vector<segment<T>> getSegments() const { return Segs; }
	const std::map<point<T>, PointInfo, PointCmp> &getPointInfoMap() const {
	return PointInfoMap;
	}
	};

view raw SegmentIntersection.cpp hosted with ❤ by GitHub

最後是測試的部分，以下圖做為測試範例：

將該圖轉換成我們接受的input如下：

10

-2 7 2 0

-2 7 -2 0

-2 6 2 5

-2 6 2 2

-2 4 2 7

-2 4 2 2

-2 4 4 1

-2 0 2 2

0 1 0 1

0 3 4 1

最後附上測試程式碼，需要的話可以自己執行看看：

	#include <cmath>
	#include <iostream>
	const double EPS = 1e-9;
	struct Double {
	double d;
	Double(double d = 0) : d(d) {}
	Double operator-() const { return -d; }
	Double operator+(const Double &b) const { return d + b.d; }
	Double operator-(const Double &b) const { return d - b.d; }
	Double operator*(const Double &b) const { return d * b.d; }
	Double operator/(const Double &b) const { return d / b.d; }
	Double operator+=(const Double &b) { return d += b.d; }
	Double operator-=(const Double &b) { return d -= b.d; }
	Double operator*=(const Double &b) { return d *= b.d; }
	Double operator/=(const Double &b) { return d /= b.d; }
	bool operator<(const Double &b) const { return d - b.d < -EPS; }
	bool operator>(const Double &b) const { return d - b.d > EPS; }
	bool operator==(const Double &b) const { return fabs(d - b.d) <= EPS; }
	bool operator!=(const Double &b) const { return fabs(d - b.d) > EPS; }
	bool operator<=(const Double &b) const { return d - b.d <= EPS; }
	bool operator>=(const Double &b) const { return d - b.d >= -EPS; }
	friend std::ostream &operator<<(std::ostream &os, const Double &db) {
	return os << db.d;
	}
	friend std::istream &operator>>(std::istream &is, Double &db) {
	return is >> db.d;
	}
	};
	void getResult(const SegmentIntersection<Double> &SI) {
	std::cout << R"(-----------------Info---------------------
	U: the set of segments start at the point.
	C: the set of segments contain the point.
	L: the set of segments end at the point.
	------------------------------------------
	)";
	for (auto p : SI.getPointInfoMap()) {
	std::cout << "(" << p.first.x << ',' << p.first.y << "):\n";
	std::cout << " U:\n";
	for (auto s : p.second.U) {
	std::cout << " (" << s->p1.x << ',' << s->p1.y << ") (" << s->p2.x
	<< ',' << s->p2.y << ")\n";
	}
	std::cout << " C:\n";
	for (auto s : p.second.C) {
	std::cout << " (" << s->p1.x << ',' << s->p1.y << ") (" << s->p2.x
	<< ',' << s->p2.y << ")\n";
	}
	std::cout << " L:\n";
	for (auto s : p.second.L) {
	std::cout << " (" << s->p1.x << ',' << s->p1.y << ") (" << s->p2.x
	<< ',' << s->p2.y << ")\n";
	}
	std::cout << '\n';
	}
	}
	int main() {
	int n;
	std::cin >> n;
	std::vector<segment<Double>> segments;
	while (n--) {
	point<Double> a, b;
	std::cin >> a.x >> a.y >> b.x >> b.y;
	segments.emplace_back(a, b);
	}
	SegmentIntersection<Double> SI(segments);
	getResult(SI);
	return 0;
	}

view raw main.cpp hosted with ❤ by GitHub