Maximum Triangle Area (Part 3 and final)

After second part, our algorithm became something like this -
for each 0 <= i < n, initialize A = P_i, B = P_{(i + 1) % n}, C = P_{(i + 2) % n}, then move C counter clockwise until area non-decreasing, move B to its adjacent point counter clockwise, if area non-decreasing, move C again until area non decreasing and so on. You will get the maximum area for each A and output the maximum of these areas. For a given point starting point P_i, if the maximum area is ABC. We can say ABC is 'stable' triangle for A, in other word s(p_i, P_{i + 1}, P_{i + 2}) = ABC, that means we cannot move B or C any further by not decreasing the area.

Now, after getting a maximum ABC for a fixed first point P_i, we will start with P_{i + 1}, P_{i + 2} and P_{i + 3}. We get A'B'C' as a stable triangle i.e. s(p_i, p_{i + 1}, p_{i + 2}) = A'B'C'. Hence our algorithm is O(n ^ 2). But it can be proved that, rather than, start from P_{i + 1} as a second point and P_{i + 2} as a third point, we can start with B and C as a second and third point and get the same result, i.e. s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C), where s(P_i, P_{i + 1}, P_{i + 2}) = ABC.

Lets assume that we are wrong, and s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = A'B'C' > s(p_i, P_{i + 1}, P_{i + 2}).

Case 1:
B' < B and C' = C


We know that,
ABC >= AB'C (otherwise largest triangle for A would be AB'C)
=> ABC + BB'A >= AB'C + BB'A
=> AB'BC >= AB'C + BB'A
=> AB'C + BB'C >= AB'C + BB'A
=> BB'C >= BB'A
=> BB'C >= BB'A' (for base BB' area is non decreasing from BB'C to BB'A, so it
will be non decreasing from BB'A to BB'C)
=> BB'C + A'B'C >= BB'A' + A'B'C
=> A'B'BC >= BB'A' + A'B'C
=> A'BC + BB'A' >= BB'A' + A'B'C
=> A'BC >= A'B'C
=> A'BC >= A'B'C' ...... [1] (as C = C')

So, for case 1, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

Case 2:
Let's say, A'B'C' is the largest triangle for first point A'.
Where, B' < B and C' > C. See the following image.



For first point A, ABC is the largest trianlge,
So, ABC >= AB'C
=> ABC + BB'C >= ABC + BB'C
=> ABC + BB'C >= AB'BC
=> ABC + BB'C >= ABC + BB'A
=> BB'C >= BB'A
i.e. BB'A <= BB'C Hence, BB'A <= BB'C' (otherwise, BB'A > BB'C and BB'A <= BB'C, which is not possible for a convex polygon) So, BB'A' <= BB'C' ... [2] (as the triangle area non increasing by moving C' to A, it will be non increasing from A to A') Now as we are considering, for first point A', the largest area triangle is A'B'C'. So, A'B'C' > A'BC'
=> A'B'C' + BB'A' > A'BC' + BB'A
=> A'B'C' + BB'A' > A'B'BC'
=> A'B'C' + BB'A' > A'B'C' + BB'C'
=> BB'A' > BB'C', this contradicts to [2].

So it is not possible to have, A'B'C' > A'BC' i.e. for case 2, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

Case 3:
Let's say, A'B'C' is the largest triangle for fixed first point A', where,
B' < B and C' < C

As ABC is the largest triangle having first point A,
either ABC' >= AB'C' (by moving B' to B, area increases) .. (3A)
or AB'C >= AB'C' (by moving C' to C, area increases) .. (3B)

Case 3A
ABC' >= AB'C'
=> ABC' + BB'A >= AB'C' + BB'A
=> AB'BC' >= AB'C' + BB'A
=> AB'C' + BB'C' >= AB'C' + BB'A
=> BB'C' >= BB'A
=> BB'C' >= BB'A' (As by moving C' to A area non increasing, moving A to A' will area be non increasing)
=> BB'C + A'B'C' >= BB'A' + A'B'C'
=> BB'A'C' >= BB'A' + A'B'C'
=> BB'A' + A'BC' >= BB'A' + A'B'C'
=> A'BC >= A'B'C'

So it is not possible to have, A'B'C' > A'BC' i.e. for case 3A, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

Case 3B
AB'C >= AB'C'
=> AB'C + CC'B' >= AB'C' + CC'B'
=> AB'C'C >= AB'C' + CC'B'
=> AB'C' + CC'A >= AB'C' + CC'B'
=> CC'A >= CC'B'
=> CC'A' >= CC'B' (As moving A to B', area decreases, it's not possible moving A' to B' area increases)
=> CC'A' + A'B'C' >= CC'B' + A'B'C'
=> A'B'C'C >= CC'B' + A'B'C'
=> CC'B' + A'B'C >= CC'B' + A'B'C'
=> A'B'C >= A'B'C'

So it is not possible to have, A'B'C' > A'B'C i.e. for case 3B, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

Case 4:
C' < C and B' = B

Here,
ABC >= ABC'
=> ABC + CC'A >= ABC' + CC'A
=> ABC + CC'A >= ABC'C
=> ABC + CC'A >= ABC + CC'B
=> CC'A >= CC'B
=> CC'A' >= CC'B
=> CC'A' + A'BC' >= CC'B + A'BC'
=> A'BC'C >= CC'B + A'BC'
=> A'BC + CC'B >= CC'B + A'BC'
=> A'BC >= A'BC'
=> A'BC >= A'B'C' (as B' = B)

So it is not possible to have, A'B'C' > A'BC i.e. for case 3B, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

Case 5:
B' > B and C' < C


Now,
A'BC >= A'BC' (From case 4)
=> A'BC + CC'B >= A'BC' + CC'B
=> A'BC'C >= A'BC' + CC'B
=> A'BC' + CC'A' >= A'BC' + CC'B
=> CC'A' >= CC'B
=> CC'A' >= CC'B'
=> CC'A' + A'B'C' >= CC'B' + A'B'C'
=> CC'B'A' >= CC'B' + A'B'C'
=> CC'B' + A'B'C >= CC'B' + A'B'C'
=> A'B'C >= A'B'C'

So it is not possible to have, A'B'C' > A'B'C i.e. for case 5, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C).

So, from all of the 5 possible cases, we found that, s(P_{i + 1}, P_{i + 2}, P_{i + 3}) = s(P_{i + 1}, B, C). So our algorithm will still work, if we start our iteration for first point i = i + 1, without changing the position 2nd and 3rd point. We will move first point from 0 to n - 1, and 2nd and 3rd point will also not move more than n times. So our algorithm in worst case is O (3 * n) = O(n)

My simple implementation in C++.

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<vector>
#include<string>
#include<algorithm>
#include<set>
#include<map>
#include<cassert>
#include<sstream>
#include<queue>
#include<stack>
#include<iostream>
using namespace std;

struct Point {
    int x;
    int y;
    Point(int _x, int _y) : x(_x), y(_y) {}
    Point() : x(0), y(0) {}
    bool operator <(const Point& p) const {
        return x < p.x || (x == p.x && y < p.y);
    }
};

bool left_turn(const Point& p1, const Point& p2, const Point& p3) {
    return (p2.y - p1.y) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.y - p1.y) > 0;
}

// Returns list of points of convex hull in counter clockwise
// The last point and first point will be equal
vector<Point> convex_hull(vector<Point> p) {
    vector<Point> ret;
    sort(p.begin(), p.end());
    // build lower hull
    for (int i = 0; i < p.size(); ++i) {
        while (ret.size() >= 2 &&
               left_turn(ret[ret.size() - 2], ret[ret.size() - 1], p[i])) {
            ret.pop_back();
        }
        ret.push_back(p[i]);
    }
    int lower_hull_size = ret.size();
    // build upper hull
    for (int i = p.size() - 2; i >= 0; --i) {
        while (ret.size() - lower_hull_size >= 1 &&
               left_turn(ret[ret.size() - 2], ret[ret.size() - 1], p[i])) {
            ret.pop_back();
        }
        ret.push_back(p[i]);
    }
    return ret;
}

double triangle_area (const Point& p1, const Point& p2, const Point& p3) {
    return abs((double) p1.x * p2.y +
               (double) p2.x * p3.y +
               (double) p3.x * p1.y -
               (double) p1.y * p2.x -
               (double) p2.y * p3.x -
               (double) p3.y * p1.x) / 2.0;
}

int main (void)
{
    int n;
    while (scanf("%d", &n) && n != -1) {
        vector<Point> p;
        for (int i = 0; i < n; ++i) {
            int x, y;
            scanf("%d%d", &x, &y);
            p.push_back(Point(x, y));
        }
        p = convex_hull(p); p.pop_back();
        int i = 0;
        int j = i + 1;
        int k = j + 1;
        double res = 0.;
        while (true) {
            double area = triangle_area(p[i], p[j], p[k]);
            while (true) {
                while (true) {
                    int nk = (k + 1) % n;
                    double narea = triangle_area(p[i], p[j], p[nk]);
                    if (narea >= area) {
                        k = nk;
                        area = narea;
                    } else {
                        break;
                    }
                }
                int nj = (j + 1) % n;
                double narea = triangle_area(p[i], p[nj], p[k]);
                if (narea >= area) {
                    j = nj;
                    area = narea;
                } else {
                    break;
                }
            }
            if (area > res) res = area;
            i++;
            if (i == j) j = (j + 1) % n;
            if (j == k) k = (k + 1) % n;
            if (i == n) break;
        }
        printf("%.2lf\n", res);
    }
    return 0;
}