Problem: Problem D Making Perimeter of the Convex Hull Shortest

The convex hull of a set of three or more planar points is, when not all of them are on one line, the convex polygon with the smallest area that has all the points of the set on its boundary or in its inside. Your task is, given positions of the points of a set, to find how much shorter the perimeter of the convex hull can be made by excluding two points from the set.

The figures below correspond to the three cases given as Sample Input 1 to 3. Encircled points are excluded to make the shortest convex hull depicted as thick dashed lines.

<image> | <image> | <image>
---|---|---
Sample Input 1 | Sample Input 2 | Sample Input 3

Input

The input consists of a single test case in the following format.


$n$
$x_1$ $y_1$
...
$x_n$ $y_n$


Here, $n$ is the number of points in the set satisfying $5 \leq n \leq 10^5$. For each $i$, ($x_i, y_i$) gives the coordinates of the position of the $i$-th point in the set. $x_i$ and $y_i$ are integers between $-10^6$ and $10^6$, inclusive. All the points in the set are distinct, that is, $x_j \ne x_k$ or $y_j \ne y_k$ holds when $j \ne k$. It is guaranteed that no single line goes through $n - 2$ or more points in the set.

Output

Output the difference of the perimeter of the convex hull of the original set and the shortest of the perimeters of the convex hulls of the subsets with two points excluded from the original set. The output should not have an error greater than $10^{-4}$.

Sample Input 1


10
-53 62
-19 58
-11 11
-9 -22
45 -7
37 -39
47 -58
-2 41
-37 10
13 42


Sample Output 1


72.96316928


Sample Input 2


10
-53 62
-19 58
-11 11
-9 -22
45 -7
43 -47
47 -58
-2 41
-37 10
13 42


Sample Output 2


62.62947992


Sample Input 3


10
-53 62
-35 47
-11 11
-9 -22
45 -7
43 -47
47 -58
-2 41
-37 10
13 42


Sample Output 3


61.58166534






Example

Input

10
-53 62
-19 58
-11 11
-9 -22
45 -7
37 -39
47 -58
-2 41
-37 10
13 42


Output

72.96316928

Solution (in CPP):

Solution:
#include <bits/stdc++.h>

using namespace std;
using uint = unsigned int;
using ll = long long;
using ull = unsigned long long;
constexpr ll TEN(int n) { return (n==0) ? 1 : 10*TEN(n-1); }
template<class T> using V = vector<T>;
template<class T> using VV = V<V<T>>;

using D = double;
const D PI = acos(D(-1)), EPS = 1e-10;

struct Pt2 {
    D x, y;
    Pt2() {}
    Pt2(D _x, D _y) : x(_x), y(_y) {}
    Pt2 operator+(const Pt2 &r) const { return Pt2(x+r.x, y+r.y); }
    Pt2 operator-(const Pt2 &r) const { return Pt2(x-r.x, y-r.y); }
    Pt2 operator*(const Pt2 &r) const { return Pt2(x*r.x-y*r.y, x*r.y+y*r.x); }
    Pt2& operator+=(const Pt2 &r) { return *this=*this+r; }
    Pt2& operator-=(const Pt2 &r) { return *this=*this-r; }
    Pt2& operator*=(const Pt2 &r) { return *this=*this*r; }
    
    Pt2 operator-() const { return Pt2(-x, -y); }

    D abs() const { return sqrt(x*x + y*y); }
    D rabs() const { return max(x, y); } // robust abs
    D arg() const { return atan2(y, x); }

    pair<D, D> to_pair() const { return make_pair(x, y); }
    static Pt2 polar(D le, D th) { return Pt2(le*cos(th), le*sin(th)); }
};
template<class T> ostream& operator<<(ostream& os, const Pt2 &p) {
    os << "(" << p.x << ", " << p.y << ")";
    return os;
}
using P = Pt2;


int sgn(D a) {
    if (abs(a) <= EPS) return 0;
    return (a < 0) ? -1 : 1;
}
//relative sign
int rsgn(D a, D f) {
    if (abs(a) <= f*EPS) return 0;
    return (a < 0) ? -1 : 1;
}
bool near(P a, P b) { return !sgn((a-b).abs()); }

bool lessP(P l, P r) {
    if (sgn(r.x-l.x)) return l.x < r.x;
    if (sgn(r.y-l.y)) return l.y < r.y;
    return false;
}

D cross(P a, P b) { return a.x*b.y - a.y*b.x; }
D dot(P a, P b) { return a.x*b.x + a.y*b.y; }

// -2, -1, 0, 1, 2 : front, clock, on, cclock, back
int ccw(P b, P c) {
    int s = rsgn(cross(b, c), b.rabs());
    if (s) return s;
    if (!sgn(c.rabs()) || !sgn((c-b).rabs())) return 0;
    if (dot(b, c) < 0) return 2;
    if (dot(-b, c-b) < 0) return -2;
    return 0;
}
int ccw(P a, P b, P c) { return ccw(b-a, c-a); }


using Pol = V<P>;

Pol convex_up(const Pol &p, int l, int r, int e1 = -1, int e2 = -1) {
    Pol up;
    for (int i = l; i <= r; i++) {
        if (i == e1 || i == e2) continue;
        P d = p[i];
        while (up.size() > 1 && ccw(up[up.size()-2], up[up.size()-1], d) == 1) up.pop_back();
        up.push_back(d);
    }
    return up;
}

Pol convex_down(const Pol &p, int l, int r, int e1 = -1, int e2 = -1) {
    Pol down;
    for (int i = l; i <= r; i++) {
        if (i == e1 || i == e2) continue;
        P d = p[i];
        while (down.size() > 1 && ccw(down[down.size()-2], down[down.size()-1], d) == -1) down.pop_back();
        down.push_back(d);
    }
    return down;
}

D line_len(const Pol &p) {
    int n = int(p.size());
    D sm = 0;
    for (int i = 0; i < n-1; i++) {
        sm += (p[i+1]-p[i]).abs();
    }
    return sm;
}

bool operator<(const P &x, const P &y) { return lessP(x, y); }


D solve(V<P> v, bool one = false) {
    int n = int(v.size());
    map<P, int> mp;
    for (int i = 0; i < n; i++) mp[v[i]] = i;
    
    auto up = convex_up(v, 0, n-1);
    auto down = convex_down(v, 0, n-1);

    int un = int(up.size()), dn = int(down.size());
    V<D> bufu, bufd;
    for (int i = 1; i < un-1; i++) {
        D base = (up[i-1]-up[i]).abs() + (up[i]-up[i+1]).abs();
        bufu.push_back(base - line_len(convex_up(v, mp[up[i-1]], mp[up[i+1]], mp[up[i]])));
    }

    for (int i = 1; i < dn-1; i++) {
        D base = (down[i-1]-down[i]).abs() + (down[i]-down[i+1]).abs();
        bufd.push_back(base - line_len(convex_down(v, mp[down[i-1]], mp[down[i+1]], mp[down[i]])));
    }

    D ans = -1e100;
    D uma = -1e100, dma = -1e100;
    for (int i = 1; i < un-2; i++) {
        ans = max(ans, uma + bufu[i]);
        uma = max(uma, bufu[i-1]);
    }
    for (int i = 1; i < dn-2; i++) {
        ans = max(ans, dma + bufd[i]);
        dma = max(dma, bufd[i-1]);
    }
    for (auto d: bufu) uma = max(uma, d);
    for (auto d: bufd) dma = max(dma, d);

    if (one) return max(uma, dma);

    ans = max(ans, uma + dma);

    for (int i = 1; i < un-1; i++) {
        D base = (up[i-1]-up[i]).abs() + (up[i]-up[i+1]).abs();
        auto npo = convex_up(v, mp[up[i-1]], mp[up[i+1]], mp[up[i]]);
        base -= line_len(npo);
        int m = int(npo.size());
        for (int j = 1; j < m-1; j++) {
            D base2 = (npo[j-1]-npo[j]).abs() + (npo[j]-npo[j+1]).abs();
            ans = max(ans, base + base2 - line_len(convex_up(v, mp[npo[j-1]], mp[npo[j+1]], mp[npo[j]], mp[up[i]])));
        }
    }

    for (int i = 1; i < dn-1; i++) {
        D base = (down[i-1]-down[i]).abs() + (down[i]-down[i+1]).abs();
        auto npo = convex_down(v, mp[down[i-1]], mp[down[i+1]], mp[down[i]]);
        base -= line_len(npo);
        int m = int(npo.size());
        for (int j = 1; j < m-1; j++) {
            D base2 = (npo[j-1]-npo[j]).abs() + (npo[j]-npo[j+1]).abs();
            ans = max(ans, base + base2 - line_len(convex_down(v, mp[npo[j-1]], mp[npo[j+1]], mp[npo[j]], mp[down[i]])));
        }
    }

    for (int i = 1; i < un-2; i++) {
        D base = (up[i-1]-up[i]).abs() + (up[i]-up[i+1]).abs() + (up[i+1]-up[i+2]).abs();
        ans = max(ans, base - line_len(convex_up(v, mp[up[i-1]], mp[up[i+2]], mp[up[i]], mp[up[i+1]])));
    }
    for (int i = 1; i < dn-2; i++) {
        D base = (down[i-1]-down[i]).abs() + (down[i]-down[i+1]).abs() + (down[i+1]-down[i+2]).abs();
        ans = max(ans, base - line_len(convex_down(v, mp[down[i-1]], mp[down[i+2]], mp[down[i]], mp[down[i+1]])));
    }

    return ans;
}


D naive(V<P> v) {
    int n = int(v.size());
    auto up = convex_up(v, 0, n-1);
    auto down = convex_down(v, 0, n-1);

    return line_len(up) + line_len(down);
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    cout << setprecision(20);

    int n;
    cin >> n;
    V<P> v(n);
    for (int i = 0; i < n; i++) {
        D x, y;
        cin >> x >> y;
        v[i] = P(x, y);
        v[i] *= P::polar(1, 13.4514);
    }
    sort(begin(v), end(v));
    
    D ma = -1;
    D base = naive(v);
    ma = max(ma, base-naive(V<P>(begin(v), end(v)-2)));
    ma = max(ma, base-naive(V<P>(begin(v)+1, end(v)-1)));
    ma = max(ma, base-naive(V<P>(begin(v)+2, end(v))));

    ma = max(ma, base-naive(V<P>(begin(v), end(v)-1))+solve(V<P>(begin(v), end(v)-1), true));
    ma = max(ma, base-naive(V<P>(begin(v)+1, end(v)))+solve(V<P>(begin(v)+1, end(v)), true));

    ma = max(ma, solve(v));
    cout << ma << endl;

    return 0;
}