/// online batch
#define inf 1e18
struct Line {
mutable long double a, b, x;
Line(long double _a, long double _b) : a(_a), b(_b), x(inf) { }
Line(long double x) : x(x), a(inf), b(inf) { }
bool operator< (const Line & q) const {
if (q.a == inf || a == inf)
return x < q.x;
if (a == q.a)
return b < q.b;
return a > q.a;
}
long double intersect(const Line & q, bool change = 0) const {
if (change)
x = (q.b - b) / (a - q.a);
return (q.b - b) / (a - q.a);
}
};
struct Batch : set <Line>
{
long double Query(long double x) {
iterator it = lower_bound(x);
return x * it->a + it->b;
}
void Add(Line x) {
if (find(x) != end())
return;
iterator it = insert(x).first;
if (bad(it) || (it != begin() && it->a == prev(it)->a)) {
erase(it);
return;
}
if (next(it) != end() && it->a == next(it)->a)
erase(next(it));
while (next(it) != end() && bad(next(it)))
erase(next(it));
while (it != begin() && bad(prev(it)))
erase(prev(it));
if (it != begin())
prev(it)->intersect(*it, 1);
if (next(it) != end())
it->intersect(*next(it), 1);
else
it->x = inf;
}
bool bad(iterator it) {
if (it == begin() || next(it) == end())
return false;
return (prev(it)->intersect(*(next(it))) >= it->intersect(*next(it)));
}
};
/** Order statistic set
* Author: Simon Lindholm
* Date: 2016-03-22
* License: CC0
* Source: hacKIT, NWERC 2015
* Description: A set (not multiset!) with support for finding the n'th
* element, and finding the index of an element.
* Time: O(\log N)
*/
#pragma once
#include <bits/extc++.h> /** keep-include */
using namespace __gnu_pbds;
template <class T>
using Tree = tree<T, null_type, less<T>, rb_tree_tag,
tree_order_statistics_node_update>;
void example() {
Tree<int> t, t2; t.insert(8);
auto it = t.insert(10).first;
assert(it == t.lower_bound(9));
assert(t.order_of_key(10) == 1);
assert(t.order_of_key(11) == 2);
assert(*t.find_by_order(0) == 8);
t.join(t2); // assuming T < T2 or T > T2, merge t2 into t
}
/** Min rotation
* Author: Stjepan Glavina
* Source: https://github.com/stjepang/snippets/blob/master/min_rotation.cpp
* Description: Finds the lexicographically smallest rotation of a string.
* Time: O(N)
* Status: Fuzz-tested
* Usage:
* rotate(v.begin(), v.begin()+MinRotation(v), v.end());
*/
#pragma once
int MinRotation(string s) {
int a = 0, n = s.size(); s += s;
for (int b = 0; b < n; ++b)
for (int i = 0; i < n; ++i) {
if (a + i == b || s[a + i] < s[b + i]) {
b += max(0, i - 1); break;
}
if (s[a + i] > s[b + i]) { a = b; break; }
}
return a;
}
/** Suffix Array + LCP
* Author: Lukas Polacek
* Date: 2009-10-27
* License: CC0
* Source: folklore and Linear-time longest-common-prefix
* computation in suffix arrays and its applications (2001).
* Description: Builds suffix array for a string. $a[i]$ is
* the starting index of the suffix which is $i$-th in the
* sorted suffix array. The {\tt lcp} function calculates
* longest common prefixes for indices.
* Can also sort cyclic permutations
* Time: $O(N \log N)$ where $N$ is the length of the string
* for creation of the SA. $O(log)$ for LCP.
* Memory: $O(N) / O(N \log N)$
* Status: NOT TESTED
*/
#pragma once
struct SuffixArray {
int n, csz;
vector<vector<int>> classes;
vector<int> cnt, order, oldc, newc, left;
string str;
SuffixArray(string s, bool cyclic) :
n(s.size() + !cyclic), csz(max(n, 256)), cnt(csz),
order(n), oldc(n), newc(n), left(n), str(s) {
if (!cyclic) str += '\0';
}
vector<int> Build() {
for (int i = 0; i < n; ++i) {
oldc[i] = newc[i] = str[i];
order[i] = left[i] = i;
}
for (int step = 1; step <= 2 * n; step *= 2) {
// Counting sort (can be replaced by sort with left)
// although not trivial
fill(cnt.begin(), cnt.end(), 0);
for (int i = 0; i < n; ++i) ++cnt[oldc[left[i]]];
for (int i = 1; i < csz; ++i) cnt[i] += cnt[i - 1];
for (int i = n - 1; i >= 0; --i)
order[--cnt[oldc[left[i]]]] = left[i];
newc[order[0]] = 0;
for (int i = 1; i < n; ++i) {
int now1 = order[i], last1 = order[i - 1],
now2 = (now1 + step / 2) % n,
last2 = (last1 + step / 2) % n;
newc[now1] = newc[last1] + (oldc[now1] != oldc[last1]
or oldc[now2] != oldc[last2]);
}
classes.push_back(newc);
swap(oldc, newc);
for (int i = 0; i < n; ++i) {
left[i] = (order[i] + n - step) % n;
}
}
return order;
}
int Compare(int i, int j, int len) {
for (int step = 0; len; ++step, len /= 2) {
if (len % 2 == 0) continue;
int ret = classes[step][i] - classes[step][j];
if (ret != 0) return ret < 0 ? -1 : 1;
i = (i + (1 << step)) % n;
j = (j + (1 << step)) % n;
}
return 0;
}
int GetLCP(int i, int j) {
if (i == j) return str.back() == '\0' ? n - i - 1 : n;
int ans = 0;
for (int step = classes.size() - 1; step >= 0; --step) {
if (classes[step][i] == classes[step][j]) {
i = (i + (1 << step)) % n;
j = (j + (1 << step)) % n;
ans += (1 << step);
}
}
return min(ans, n); // if cyclic
}
};
/** FFT
* Author: Lucian Bicsi
* Date: 2015-06-25
* License: GNU Free Documentation License 1.2
* Source: http://rosettacode.org/wiki/Fast_Fourier_transform
Papers about accuracy: http://www.daemonology.net/papers/fft.pdf, http://www.cs.berkeley.edu/~fateman/papers/fftvsothers.pdf
* Description: Fast Fourier transform. Also includes a function for convolution:
\texttt{conv(a, b) = c}, where $c[x] = \sum a[i]b[x-i]$. $a$ and $b$ should be of roughly equal size.
Does about 1.2s for $10^6$ elements.
Rounding the results of conv
works if $(|a| + |b|)\max(a, b) < \mathtt{\sim} 10^9$ (in theory maybe $10^6$);
you may want to use an NTT from the Number Theory chapter instead.
* Time: O(N \log N)
* Status: somewhat tested
*/
#pragma once
#include <bits/stdc++.h>
using namespace std;
struct FFTSolver {
using Complex = complex<double>;
const double kPi = 4.0 * atan(1.0);
vector<int> rev;
int __lg(int n) { return n == 1 ? 0 : 1 + __lg(n / 2); }
void compute_rev(int n, int lg) {
rev.resize(n); rev[0] = 0;
for (int i = 1; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lg - 1));
}
}
vector<Complex> fft(vector<Complex> V, bool invert) {
int n = V.size(), lg = __lg(n);
if ((int)rev.size() != n) compute_rev(n, lg);
for (int i = 0; i < n; ++i) {
if (i < rev[i])
swap(V[i], V[rev[i]]);
}
for (int step = 2; step <= n; step *= 2) {
const double ang = 2 * kPi / step;
Complex eps(cos(ang), sin(ang));
if (invert) eps = conj(eps);
for (int i = 0; i < n; i += step) {
Complex w = 1;
for (int a = i, b = i+step/2; b < i+step; ++a, ++b) {
Complex aux = w * V[b];
V[b] = V[a] - aux;
V[a] = V[a] + aux;
w *= eps;
}
}
}
return V;
}
vector<Complex> transform(vector<Complex> V) {
int n = V.size();
vector<Complex> ret(n);
Complex div_x = Complex(0, 1) * (4.0 * n);
for (int i = 0; i < n; ++i) {
int j = (n - i) % n;
ret[i] = (V[i] + conj(V[j]))
* (V[i] - conj(V[j])) / div_x;
}
return ret;
}
vector<int> Multiply(vector<int> A, vector<int> B) {
int n = A.size() + B.size() - 1;
vector<int> ret(n);
while (n != (n & -n)) ++n;
A.resize(n); B.resize(n);
vector<Complex> V(n);
for (int i = 0; i < n; ++i) {
V[i] = Complex(A[i], B[i]);
}
V = fft(move(V), false);
V = transform(move(V));
V = fft(move(V), true);
for (int i = 0; i < (int)ret.size(); ++i)
ret[i] = round(real(V[i]));
return ret;
}
};
/** Matrix Determinant
* Author: Unknown
* Date: 2014-11-27
* Source: somewhere on github
* Description: Calculates determinant using modular arithmetics.
* Modulos can also be removed to get a pure-integer version.
* Status: bruteforce-tested for N <= 3, mod <= 7
* Time: $O(N^3)$
*/
#pragma once
using int64 = int64_t;
const int64 kMod = 12345;
int64 IntDeterminant(vector<vector<int64>>& M) {
int n = M.size(); int64 ans = 1;
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j) {
while (M[j][i] != 0) { // gcd step
int64 t = M[i][i] / M[j][i];
if (t) for (int k = i; k < n; ++k)
M[i][k] = (M[i][k] - M[j][k] * t) % kMod;
swap(M[i], M[j]);
ans *= -1;
}
}
ans = ans * a[i][i] % mod;
if (!ans) return 0;
}
return (ans + kMod) % kMod;
}
/** Matrix Inverse
* Author: Max Bennedich
* Date: 2004-02-08
* Description: Invert matrix $A$. Returns rank; result is stored in $A$ unless singular (rank < n).
* Can easily be extended to prime moduli; for prime powers, repeatedly
* set $A^{-1} = A^{-1} (2I - AA^{-1})\ (\text{mod }p^k)$ where $A^{-1}$ starts as
* the inverse of A mod p, and k is doubled in each step.
* Time: O(n^3)
* Status: Slightly tested
*/
#pragma once
int matInv(vector<vector<double>>& A) {
int n = sz(A); vi col(n);
vector<vector<double>> tmp(n, vector<double>(n));
rep(i,0,n) tmp[i][i] = 1, col[i] = i;
rep(i,0,n) {
int r = i, c = i;
rep(j,i,n) rep(k,i,n)
if (fabs(A[j][k]) > fabs(A[r][c]))
r = j, c = k;
if (fabs(A[r][c]) < 1e-12) return i;
A[i].swap(A[r]); tmp[i].swap(tmp[r]);
rep(j,0,n)
swap(A[j][i], A[j][c]), swap(tmp[j][i], tmp[j][c]);
swap(col[i], col[c]);
double v = A[i][i];
rep(j,i+1,n) {
double f = A[j][i] / v;
A[j][i] = 0;
rep(k,i+1,n) A[j][k] -= f*A[i][k];
rep(k,0,n) tmp[j][k] -= f*tmp[i][k];
}
rep(j,i+1,n) A[i][j] /= v;
rep(j,0,n) tmp[i][j] /= v;
A[i][i] = 1;
}
/// forget A at this point, just eliminate tmp backward
for (int i = n-1; i > 0; --i) rep(j,0,i) {
double v = A[j][i];
rep(k,0,n) tmp[j][k] -= v*tmp[i][k];
}
rep(i,0,n) rep(j,0,n) A[col[i]][col[j]] = tmp[i][j];
return n;
}
/** Simplex
* Author: Stanford
* Source: Stanford Notebook
* Description: Solves a general linear maximization problem: maximize $c^T x$ subject to $Ax \le b$, $x \ge 0$.
* Returns -inf if there is no solution, inf if there are arbitrarily good solutions, or the maximum value of $c^T x$ otherwise.
* The input vector is set to an optimal $x$ (or in the unbounded case, an arbitrary solution fulfilling the constraints).
* Numerical stability is not guaranteed. For better performance, define variables such that $x = 0$ is viable.
* Usage:
* A = {{1,-1}, {-1,1}, {-1,-2}};
* b = {1,1,-4}, c = {-1,-1}, x;
* T val = LPSolver(A, b, c).solve(x);
* Time: O(NM * \#pivots), where a pivot may be e.g. an edge relaxation. O(2^n) in the general case.
* Status: seems to work?
*/
#pragma once
typedef double T; // long double, Rational, double + mod<P>...
typedef vector<T> vd;
typedef vector<vd> vvd;
const T eps = 1e-8, inf = 1/.0;
#define MP make_pair
#define ltj(X) if(s == -1 || MP(X[j],N[j]) < MP(X[s],N[s])) s=j
#define rep(i, a, b) for(int i = a; i < (b); ++i)
#define sz(x) (int)(x).size()
struct LPSolver {
int m, n; vi N, B; vvd D;
LPSolver(const vvd& A, const vd& b, const vd& c) :
m(sz(b)), n(sz(c)), N(n+1), B(m), D(m+2, vd(n+2)) {
rep(i,0,m) rep(j,0,n) D[i][j] = A[i][j];
rep(i,0,m) { B[i] = n+i; D[i][n] = -1; D[i][n+1] = b[i];}
rep(j,0,n) { N[j] = j; D[m][j] = -c[j]; }
N[n] = -1; D[m+1][n] = 1;
}
void pivot(int r, int s) {
T *a = D[r].data(), inv = 1 / a[s];
rep(i,0,m+2) if (i != r && abs(D[i][s]) > eps) {
T *b = D[i].data(), inv2 = b[s] * inv;
rep(j,0,n+2) b[j] -= a[j] * inv2;
b[s] = a[s] * inv2;
}
rep(j,0,n+2) if (j != s) D[r][j] *= inv;
rep(i,0,m+2) if (i != r) D[i][s] *= -inv;
D[r][s] = inv;
swap(B[r], N[s]);
}
bool simplex(int phase) {
int x = m + phase - 1;
for (;;) {
int s = -1;
rep(j,0,n+1) if (N[j] != -phase) ltj(D[x]);
if (D[x][s] >= -eps) return true;
int r = -1;
rep(i,0,m) {
if (D[i][s] <= eps) continue;
if (r == -1 || MP(D[i][n+1] / D[i][s], B[i])
< MP(D[r][n+1] / D[r][s], B[r])) r = i;
}
if (r == -1) return false;
pivot(r, s);
}
}
T Solve(vd &x) {
int r = 0;
rep(i,1,m) if (D[i][n+1] < D[r][n+1]) r = i;
if (D[r][n+1] < -eps) {
pivot(r, n);
if (!simplex(2) || D[m+1][n+1] < -eps) return -inf;
rep(i,0,m) if (B[i] == -1) {
int s = 0;
rep(j,1,n+1) ltj(D[i]);
pivot(i, s);
}
}
bool ok = simplex(1); x = vd(n);
rep(i,0,m) if (B[i] < n) x[B[i]] = D[i][n+1];
return ok ? D[m][n+1] : inf;
}
};
/** Miller-Rabin prime tester
* Author: Lukas Polacek
* Date: 2010-01-26
* License: CC0
* Source: TopCoder tutorial
* Description: Miller-Rabin primality probabilistic test.
* Probability of failing one iteration is at most 1/4. 15 iterations should be
* enough for 50-bit numbers.
* Time: 15 times the complexity of $a^b \mod c$.
*/
#pragma once
#include "ModMulLL.h"
using ull = unsigned long long;
bool IsPrime(ull p) {
if (p == 2) return true;
if (p == 1 || p % 2 == 0) return false;
ull s = p - 1;
while (s % 2 == 0) s /= 2;
for (int i = 0; i < 15; ++i) {
ull a = rand() % (p - 1) + 1, tmp = s;
ull mod = ModPow(a, tmp, p);
while (tmp != p - 1 && mod != 1 && mod != p - 1) {
mod = ModMul(mod, mod, p);
tmp *= 2;
}
if (mod != p - 1 && tmp % 2 == 0) return false;
}
return true;
}
/** modular sqrt
* Author: Simon Lindholm
* Date: 2016-08-31
* License: CC0
* Source: http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python/
* Description: Tonelli-Shanks algorithm for modular square roots.
* Time: O(\log^2 p) worst case, often O(\log p)
* Status: Tested for all a,p <= 10000
*/
#pragma once
#include "ModPow.h"
ll sqrt(ll a, ll p) {
a %= p; if (a < 0) a += p;
if (a == 0) return 0;
assert(modpow(a, (p-1)/2, p) == 1);
if (p % 4 == 3) return modpow(a, (p+1)/4, p);
// a^(n+3)/8 or 2^(n+3)/8 * 2^(n-1)/4 works if p % 8 == 5
ll s = p - 1;
int r = 0;
while (s % 2 == 0)
++r, s /= 2;
ll n = 2; // find a non-square mod p
while (modpow(n, (p - 1) / 2, p) != p - 1) ++n;
ll x = modpow(a, (s + 1) / 2, p);
ll b = modpow(a, s, p);
ll g = modpow(n, s, p);
for (;;) {
ll t = b;
int m = 0;
for (; m < r; ++m) {
if (t == 1) break;
t = t * t % p;
}
if (m == 0) return x;
ll gs = modpow(g, 1 << (r - m - 1), p);
g = gs * gs % p;
x = x * gs % p;
b = b * g % p;
r = m;
}
}
/** Matching weighted
* Author: Stanford
* Date: Unknown
* Source: Stanford Notebook
* Description: Min cost perfect bipartite matching. Negate costs for max cost.
* Time: O(N^3)
* Status: tested during ICPC 2015
*/
#pragma once
int sgn(double x) {return abs(x) < 1e-9 ? 0 : x > 0 ? 1 : -1;}
double MinCostMatching(vector<vector<double>> cost,
vector<int>& L, vector<int>& R)
{
int n = cost.size(), mated = 0;
vector<double> dist(n), u(n), v(n);
vector<int> par(n), seen(n);
/// construct dual feasible solution
for (int i = 0; i < n; ++i) {
u[i] = cost[i][0];
for (int j = 1; j < n; ++j)
u[i] = min(u[i], cost[i][j]);
}
for (int j = 0; j < n; ++j) {
v[j] = cost[0][j] - u[0];
for (int i = 1; i < n; ++i)
v[j] = min(v[j], cost[i][j] - u[i]);
}
/// find primal solution satisfying complementary slackness
L = R = vector<int>(n, -1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j) {
if (R[j] != -1) continue;
if (sgn(cost[i][j] - u[i] - v[j]) == 0) {
L[i] = j; R[j] = i; mated++; break;
}
}
for (; mated < n; mated++) { // until solution is feasible
int s = 0;
while (L[s] != -1) s++;
fill(par.begin(), par.end(), -1);
fill(seen.begin(), seen.end(), 0);
for (int k = 0; k < n; ++k)
dist[k] = cost[s][k] - u[s] - v[k];
int j = 0;
while (true) { /// find closest
j = -1;
for (int k = 0; k < n; ++k) {
if (seen[k]) continue;
if (j == -1 || dist[k] < dist[j]) j = k;
}
seen[j] = 1;
int i = R[j];
if (i == -1) break;
for (int k = 0; k < n; ++k) { /// relax neighbors
if (seen[k]) continue;
auto new_dist = dist[j] + cost[i][k] - u[i] - v[k];
if (dist[k] > new_dist) {
dist[k] = new_dist;
par[k] = j;
}
}
}
/// update dual variables
for (int k = 0; k < n; ++k) {
if (k == j || !seen[k]) continue;
auto w = dist[k] - dist[j];
v[k] += w, u[R[k]] -= w;
}
u[s] += dist[j];
/// augment along path
while (par[j] != -1) {
int p = par[j];
R[j] = R[p]; L[R[j]] = j; j = p;
}
R[j] = s; L[s] = j;
}
double ret = 0;
for (int i = 0; i < n; ++i)
ret += cost[i][L[i]];
return ret;
}
/** Dinic flow
* Author: Lucian Bicsi
* License: CC0
* Description: Quick flow algorithm.
* Time: $O(V^2 * E)$ or $O(E * sqrt(E))$ on unit graphs
* Status: Tested on kattis and SPOJ
*/
#pragma once
struct Dinic {
struct Edge { int to, cap, flow, nxt; };
vector<Edge> edges;
vector<int> graph, at, dist;
int src = 0, dest = 1;
Dinic(int n) : graph(n, -1), dist(n, -1) {}
void add_edge(int from, int to, int cap) {
edges.push_back(Edge {to, cap, 0, graph[from]});
graph[from] = edges.size() - 1;
}
void AddEdge(int from, int to, int cap) {
add_edge(from, to, cap);
add_edge(to, from, 0);
}
bool bfs() {
queue<int> q;
fill(dist.begin(), dist.end(), -1);
dist[src] = 0; q.push(src);
while (!q.empty()) {
int node = q.front(); q.pop();
for (int i = graph[node]; i >= 0; i = edges[i].nxt) {
const auto &e = edges[i];
if (dist[e.to] == -1 && e.flow < e.cap) {
dist[e.to] = dist[node] + 1;
q.push(e.to);
}
}
}
return dist[dest] != -1;
}
int dfs(int node, int flow) {
if (flow == 0) return 0;
if (node == dest) return flow;
while (at[node] != -1) {
int eid = at[node]; const auto &e = edges[eid];
if (dist[e.to] == dist[node] + 1) {
if (int ret = dfs(e.to, min(flow, e.cap - e.flow))) {
edges[ eid ].flow += ret;
edges[eid^1].flow -= ret;
return ret;
}
}
at[node] = e.nxt;
}
return 0;
}
int Compute(int src, int dest) {
this->src = src; this->dest = dest; int ret = 0;
while (bfs()) {
at = graph;
while (int flow = dfs(src, 2e9))
ret += flow;
}
return ret;
}
};
//// Geometry
/** Point struct
* Author: Lucian Bicsi
* License: CC0
* Description: Point declaration, and basic operations.
* Status: Works fine, used a lot
*/
#pragma once
using Point = complex<double>;
const double kPi = 4.0 * atan(1.0);
const double kEps = 1e-9; // Good eps for long double is ~1e-11
#define X() real()
#define Y() imag()
double dot(Point a, Point b) { return (conj(a) * b).X(); }
double cross(Point a, Point b) { return (conj(a) * b).Y(); }
double dist(Point a, Point b) { return abs(b - a); }
Point perp(Point a) { return Point{-a.Y(), a.X()}; } // +90deg
double rotateCCW(Point a, double theta) {
return a * polar(1.0, theta); }
double det(Point a, Point b, Point c) {
return cross(b - a, c - a); }
// abs() is norm (length) of vector
// norm() is square of abs()
// arg() is angle of vector
// det() is twice the signed area of the triangle abc
// and is > 0 iff c is to the left as viewed from a towards b.
// polar(r, theta) gets a vector from abs() and arg()
void ExampleUsage() {
Point a{1.0, 1.0}, b{2.0, 3.0};
cerr << a << " " << b << endl;
cerr << "Length of ab is: " << dist(a, b) << endl;
cerr << "Angle of a is: " << arg(a) << endl;
cerr << "axb is: " << cross(a, b) << endl;
}
/** circle intersection
* Author: Simon Lindholm
* Date: 2015-09-01
* License: CC0
* Description: Computes the intersection between two circles
* Returns duplicate point when tangent
* Status: somewhat tested
*/
#pragma once
#include "Point.h"
#include <bits/stdc++.h>
struct Circle { Point c; double r; };
vector<Point> CircleIntersection(Circle& A, Circle& B) {
Point a = A.c, b = B.c;
double r1 = A.r, r2 = B.r;
Point delta = b - a;
if (abs(delta) < kEps) return {};
double r = r1 + r2, d2 = norm(delta);
double p = (d2 + r1*r1 - r2*r2) / (2.0 * d2);
double h2 = r1*r1 - p*p*d2;
if (d2 > r*r + kEps || h2 < -kEps) return {};
Point mid = a + delta * p,
per = perp(delta) * sqrt(abs(h2) / d2);
return {mid - per, mid + per};
}
/** distance to segment
* Author: Lucian Bicsi
* License: CC0
* Source:
* Description:\\
\begin{minipage}{75mm}
Returns the shortest distance between point p and the line segment from point s to e.
\end{minipage}
\begin{minipage}{15mm}
\vspace{-10mm}
\includegraphics[width=\textwidth]{../content/geometry/SegmentDistance}
\end{minipage}
* Status: NOT TESTED
* Usage:
* Point a{0, 0}, b{2, 2}, p{1, 1};
* bool onSegment = SegmentDistance(p, a, b) < kEps;
*/
#pragma once
#include "Point.h"
double SegmentDistance(Point p, Point a, Point b) {
if (a == b) return abs(p - a); // Beware of precision!!!
double d = norm(b - a);
double t = min(d, max(.0, dot(p - a, b - a)));
return abs((p - a) * d - (b - a) * t) / d;
}
// Projects point p on segment [a, b]
Point ProjectPointOnSegment(Point p, Point a, Point b) {
double d = norm(b - a);
if (d == 0) return a; // Beware of precision!!!
double r = dot(p - a, b - a) / d;
return (r < 0) ? a : (r > 1) ? b : (a + (b - a) * r);
}
Theodor Moroianu theodor.moroianu