++w++ ICPCライブラリ
++w++(amylase, y3eadgbe, atetubou(is2012))で使用する ICPC ライブラリ置き場。 (last updated: 7/13 2:20)
- 基本的には練習で書いたコードのうち汎用性の高いものを載せている感じです。
- 自由に使ってもらって構いませんが、正当性は保証しません。
- つい最近まで点と線分の距離が間違ってたりしました。そういう感じです。
- 自由に使ってもらって構いませんが、正当性は保証しません。
- TODO
- 線型計画を二段階単体法にする。
- 幾何(線分アレンジメント)
- 文字列(KMP method, Aho-Corasick algorithm)
- (ここには載せないが、蟻本から最大流・最小費用流・拡張ユークリッド・連立線型合同式を持っていく)
- 余裕があれば有名問題を増やす。
- グラフ彩色
目次
テンプレート
#include <iostream> #include <vector> #include <map> #include <algorithm> #include <queue> #include <cstring> #include <complex> using namespace std; #define rep(i,n) for(int i=0;i<n;++i) #define pb(a) push_back(a) #define mp(a,b) make_pair(a,b) #define F first #define S second #define SZ(a) (int)((a).size())
- だいたいコードの先頭にあると考えておいてね。
幾何
基本要素
// primitives
typedef long double Real;
typedef complex<Real> Point;
typedef pair<Point, Point> Line;
typedef pair<Point, Real> Circle;
typedef vector<Point> Poly;
Real eps = 1e-9;
inline istream& operator>>(istream& s, Point& p) {return s >> p.real() >> p.imag();}
inline bool near(const Point& p, const Point& q){return abs(p - q) < eps;}- Real : long double, double の切り替え用に。
- Point : 点もしくはベクトル。abs(x) で長さ。arg(x) で偏角[-pi, +pi]。a + bi = pt(a, b)
- Line : 線分・直線。通る2点（線分では両端）。
- Circle : 円。中心点と半径。
- Poly : 多角形。反時計回りに頂点を持つ。
- eps : 許容誤差。問題によって値を修正する。
operator>> : (x座標) (y座標) と来た時に一発で読み取る。
- near : 誤差を許容して同じ点かどうか見るよ。
内積・外積の大きさ
inline Real dot(const Point& p, const Point& q){return p.real() * q.real() + p.imag() * q.imag();}
inline Real cross(const Point& p, const Point& q){return p.real() * q.imag() - p.imag() * q.real();}- dot : a, b の内積。射影を取ったりするときに使う。
- cross : a, b が張る平行四辺形の面積。三角形の面積とか。a, b を通る直線上の点 x が cross(x, b) = cross(a, b) を満たすことから、abの法線を軸としたパラメータとの解釈もある。
3点の位置関係
// ccw_b : old version of ccw.
inline Real ccw_b(const Point& p, const Point& q, const Point& r) {return cross(q - p, r - p);}
/* ccw :
CD : counter direction
CW : clock wise
OS : on segment
CCW : counter clock wise
D : direction
*/
enum LPposit { P_CD = -2, P_CW = -1, P_OS = 0, P_CCW = 1, P_D = 2};
LPposit ccw(const Point& p, const Point& q, const Point& r) {
Real c = ccw_b(p, q, r);
if (c < -eps) return P_CW;
if (c > eps) return P_CCW;
if (dot(q - p, r - p) < -eps) return P_CD;
if (dot(p - q, r - q) < -eps) return P_D;
return P_OS;
}- ccw_b : 上で述べた cross のパラメータ解釈での値を計算。
- ccw : コメント参照!
点と直線の距離
inline Real Sabs(const Line& l) {return abs(l.first - l.second); }
inline Real LPdist(const Line& l, const Point& p) {return abs(ccw_b(l.first, l.second, p)) / Sabs(l); }
inline Real SPdist(Line l, Point p) {
Real a = abs(l.first - p);
Real b = abs(l.second - p);
Real c = Labs(l);
if (b * b + c * c > a * a && a * a + c * c > b * b){
return LPdist(l, p);
}
return min(a, b);
}- Sabs : 線分の長さ。
- LPdist : 点と直線の距離。
- SPdist : 点と線分の距離。
- L : line, S : segment, P : point
線分の交差判定
bool crossS(const Line& p, const Line& q){
return
ccw(p.first, p.second, q.first) * ccw(p.first, p.second, q.second) <= 0 &&
ccw(q.first, q.second, p.first) * ccw(q.first, q.second, p.second) <= 0;
}- 線分 p と q が交差するかどうか。
直線の交点
Point intersect(const Line& p, const Line& q) {
Point vp = p.second - p.first;
Point vq = q.second - q.first;
Point c(cross(vp, p.first), cross(vq, q.first));
return Point(cross(c, Point(vp.real(), vq.real())), cross(c, Point(vp.imag(), vq.imag()))) / cross(vp, vq);
}
多角形の面積
Real area(const Poly& p) {
Real ret = 0.0;
for(Poly::iterator it = p.begin(); it != p.end(); it++) ret += cross(*it, *(it + 1));
return ret / 2;
}- 多角形の面積。凸性は仮定しない。O(n)。
- 時計回りの時は面積が -1 倍される。abs 取ろう。
円と直線の交点
vector<Point> circle_line_intersect(line l,circle c){
vector<Point> ret;
Real di = dist(l,c.first);
Real r=c.second;
if(di+EPS > r) return ret;
Point v=(l.second-l.first);
v/=abs(v);
Point rv=v*Point(0,1);
rv*=di;
if(dist(l,c.first+rv) > di+EPS) rv = -rv;
v*=sqrt(r*r-di*di);
ret.push_back(c.first+rv-v);
ret.push_back(c.first+rv+v);
return ret;
}- 円と直線の交点を求める。
- 円が直線と接しない場合はemptyを返し、それ以外は取り敢えず二つの要素を持つvectorを返す。
円の共通接線
vector<Line> common_tangent (Circle p, Circle q) {
Real pr = p.second, qr = q.second;
Point pc = p.first, qc = q.first;
Real d = abs(pc - qc), dr = abs(pr - qr), sr = abs(pr + qr);
vector<Line> ret;
if (d > sr) {
// cross pair
Point cp = (pc * qr + qc * pr) / sr;
vector<Point> pts = tangent(cp, p), qts = tangent(cp, q);
ret.push_back(Line(pts[0], qts[0]));
ret.push_back(Line(pts[1], qts[1]));
} else if (abs(d - sr) < eps) {
// cross pair coinside
Point cp = (pc * qr + qc * pr) / sr;
ret.push_back(Line(cp, cp));
}
if (d > dr) {
// outer pair
if (abs(pr - qr) < eps) {
Point v = (qc - pc) / d;
v *= Point(0, 1);
ret.push_back(Line(pc + v, qc + v));
ret.push_back(Line(pc - v, qc - v));
} else {
Point cp = pc + (qc - pc) * pr / (pr - qr);
vector<Point> pts = tangent(cp, p), qts = tangent(cp, q);
ret.push_back(Line(pts[0], qts[0]));
ret.push_back(Line(pts[1], qts[1]));
}
} else if (abs(d - dr) < eps) {
// outer pair coinside
Point cp = (qc - pc) * pr / (pr - qr);
ret.push_back(Line(cp, cp));
}
return ret;
}vector<line> の各要素が接線。
凸包
// convex-hull
// return minimum convex polygon that contains every point in vs.
// depend on ccw
namespace std{
bool operator<(const Point& a, const Point& b) {
return a.real() != b.real() ? a.real() < b.real() : a.imag() < b.imag();
}
}
Poly convexHull(vector<Point> ps) {
int n = ps.size();
sort(ps.begin(), ps.end());
Poly ret(2 * n);
int m = 0;
for (int i = 0; i < n; i++) {
while (m >= 2 && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
ret[m++] = ps[i];
}
int t = m;
for (int i = n-2; i >= 0; i--) {
while (m >= t && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
ret[m++] = ps[i];
}
ret.resize(m - 1);
return ret;
}- Point の座標軸ソートもついてる。
- convexHull : ps に含まれる点をすべて含む最小の凸多角形を返す。spaghetthi source を参考にしたよ。
円の交点
vector<Point> CCintersect (Circle c, Circle d) {
vector<Point> ret;
const Real dist = abs(c.first - d.first);
const Real cr = c.second;
const Real dr = d.second;
if (dist > cr + dr) return ret;
if (dist < abs(cr - dr)) return ret;
const Real s = (cr + dr + dist) / 2.;
const Real area = sqrt(s * (s - cr) * (s - dr) * (s - dist));
const Real h = 2 * area / dist;
Point v = d.first - c.first; v /= abs(v);
const Point m = c.first + sqrt(cr * cr - h * h) * v;
const Point n = v * Point(0, 1);
ret.push_back(m + n * h);
ret.push_back(m - n * h);
return ret;
}
数論
ミラー・ラビン素数判定法
long long powmod(long long x, long long p, long long m){
if(p == 0) return 1;
long long rt = powmod(x, p/2, m);
if(p % 2 == 0){
return rt * rt % m;
} else {
return (rt * rt % m) * x % m;
}
}
int miller_test[9] = {2,3,5,7,11,13,17,19,23};
bool isprime(long long n){
if(n <= ISPR_MAX) return !isnpr[n];
long long d = n-1, s = 0;
while(d%2 == 0){
d /= 2;
s++;
}
for(int i=0; i<9; i++){
bool iscomp = true;
li x = powmod(miller_test[i], d, n);
iscomp = iscomp && (x % n != 1);
for(int r=0; r<s; r++){
iscomp = iscomp && (x % n != n-1);
x = x * x % n;
}
if(iscomp) return false;
}
return true;
}
ポラード・ロー素因数分解法
long long myrand(long long c, long long n, long long x){
return (x * x + c) % n;
}
long long pollard(long long n){
long long x = 2, y = 2, d = 1, c = rand();
while(d == 1){
x = myrand(c, n, x);
y = myrand(c, n, myrand(c, n, y));
d = __gcd(abs(x-y), n);
}
return d;
}
素数列挙
- 改良エラトステネスの篩
- 普通に篩うより2倍ぐらい早い
vector<int> getPrime(const int n) {
const int ub = (n - 1) / 2;
const int sqrtub = (sqrt(n) - 1) / 2;
vector<int> res;
if (n <= 1) return res;
res.push_back(2);
bool *isNotPrime = new bool[ub + 1];
for (int i = 0; i <= ub; i++)isNotPrime[i] = (i % 3 == 1);
isNotPrime[1] = false;
for (int i = 2; i <= sqrtub; i++) {
if (!isNotPrime[i]) {
int d = i * 2 + 1;
for (int j = 3 * i + 1; j <= ub; j += d) {
isNotPrime[j] = true;
}
}
}
for (int i = 1; i <= ub; i++) {
if (!isNotPrime[i]) {
res.push_back(i * 2 + 1);
}
}
delete[] isNotPrime;
return res;
}
数学
行列乗算
typedef vector<vector<li> > matrix;
const li mod = 1000000007;
matrix ident(const int& n) {
matrix ret(n, vi(n, 0));
rep(i, n) ret[i][i] = 1;
return ret;
}
matrix matadd(const matrix& p, const matrix& q) {
int n = p.size(), m = p[0].size();
matrix ret = p;
rep(i, n) rep(j, m) ret[i][j] = (p[i][j] + q[i][j]) % mod;
return ret;
}
matrix matmul(const matrix& p, const matrix& q) {
int n = p.size(), m = q[0].size(), l = p[0].size();
matrix ret = matrix(n, vi(m, 0));
rep(i, n) rep(j, m) rep(k, l) ret[i][j] = (ret[i][j] + p[i][k] * q[k][j] % mod) % mod;
return ret;
}
matrix matpow(const matrix& p, const li& x) {
if (x == 0) return ident(p.size());
if (x == 1) return p;
matrix ret = matpow(p, x / 2);
ret = matmul(ret, ret);
if (x % 2) ret = matmul(ret, p);
return ret;
}- 2次元配列で行列を表現。エラーチェックなし。mod とってます。
mod 2 上の行列演算
- bit-wise parallelized version. 爆速です。
畳み込み
typedef double Real;
typedef complex<Real> Complex;
vector<Complex> fastFourierTransform(const vector<Complex>& x) {
if (__builtin_popcount(x.size()) > 1) {
cerr << "fastFourierTransform: x.size() should be power of 2." << endl;
return x;
}
const int n = x.size();
if (n == 1) {
return x;
}
vector<Complex> even(n / 2), odd(n / 2), ret(n);
for (int i = 0; i < n; ++i) {
if (i % 2 == 0) {
even[i / 2] = x[i];
} else {
odd[i / 2] = x[i];
}
}
even = fastFourierTransform(even);
odd = fastFourierTransform(odd);
for (int i = 0; i < n; ++i) {
int j = i % (n / 2);
ret[i] = even[j] + odd[j] * exp(Complex(0, -2 * M_PI * i / n));
}
return ret;
}
vector<Complex> inverseFastFourierTransform(vector<Complex> x) {
const int n = x.size();
for (int i = 0; i < n; ++i) {
x[i] = conj(x[i]);
}
x = fastFourierTransform(x);
for (int i = 0; i < n; ++i) {
x[i] = conj(x[i]) / Complex(n);
}
return x;
}
int conv_minpow2(int x) {
int ret = 1;
while (ret < x) ret <<= 1;
return ret;
}
vector<Complex> convolution(vector<Complex> x,
vector<Complex> y) {
int l = 2 * max(conv_minpow2(x.size()), conv_minpow2(y.size()));
x.resize(l);
y.resize(l);
x = fastFourierTransform(x);
y = fastFourierTransform(y);
vector<Complex> ret(l);
for (int i = 0; i < l; i++) {
ret[i] = x[i] * y[i];
}
ret = inverseFastFourierTransform(ret);
while (not ret.empty() && ret.back() == Complex(0)) ret.pop_back();
return ret;
}- 高速フーリエ変換による畳み込み演算。O(n log n)。
- 多項式の掛け算がナイーブにやるとO(n^2)だが畳込みを使うと計算量が落ちたり。
有名問題
線型計画法
// linear programming by simplex method
// solve maximize(cx) s.t. Ax <= a
// if not bounded, return -1.
double simplex(vector<double>& c,
vector<vector<double> >& A,
vector<double>& a) {
// introduce slack variable
int n = c.size(), m = a.size(), dim = n + m;
vector<int> N(n), B(m);
for (int i = 0; i < n; ++i) N[i] = i;
for (int i = 0; i < m; ++i) {
B[i] = i + n;
c.push_back(0);
for (int j = 0; j < n; j++) A[i][j] *= -1;
for (int j = 0; j < m; j++) A[i].push_back(0);
}
double ans = 0;
while (true) {
// check optimized or not
int s = -1;
for (int i = 0; i < n; i++) if (c[N[i]] > 0) s = N[i];
if (s < 0) break;
// check bounded or not
double bound = 1e300;
int r = -1;
for (int i = 0; i < m; i++) {
if (A[i][N[s]] < 0) {
double nbound = -a[i] / A[i][N[s]];
if (nbound < bound) {
bound = nbound;
r = i;
}
}
}
if (r < 0) return -1;
// pivotting
for (int i = 0; i < dim; i++) if (i != N[s]) A[r][i] /= -A[r][N[s]];
a[r] /= -A[r][N[s]];
A[r][B[r]] = 1.0 / A[r][N[s]];
A[r][N[s]] = 0;
for (int i = 0; i < m; i++) {
if (i == r) continue;
for (int j = 0; j < dim; j++) if (j != N[s]) A[i][j] += A[i][N[s]] * A[r][j];
a[i] += A[i][N[s]] * a[r];
A[i][N[s]] = 0;
}
ans += c[N[s]] * a[r];
for (int i = 0; i < dim; i++) if (i != N[s]) c[i] += c[N[s]] * A[r][i];
c[N[s]] = 0;
swap(N[s], B[r]);
}
return ans;
}simplex(c, A, a) : maximize cx, where Ax <= a, x >= 0 なる線型計画問題を単体法で解く。
- x = 0 が実行可能解であることを仮定している（つまり、TODOとして二段階単体法に改善することがあげられる）。
- 発散する場合には -1 を返す。有限な最適解は必ず正になっているはずなので、正負を見れば発散したかどうかがわかる。
リンク
参考文献
- プログラミングコンテストチャレンジブック