## Saturday, August 20, 2011

### Maximum Triangle Area (Part 3 and final)

Part 1
Part 2

After second part, our algorithm became something like this -
for each 0 <= i < n, initialize A = Pi, 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 Pi, if the maximum area is ABC. We can say ABC is 'stable' triangle for A, in other word s(pi, Pi + 1, Pi + 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 Pi, we will start with Pi + 1, Pi + 2 and Pi + 3. We get A'B'C' as a stable triangle i.e. s(pi, pi + 1, pi + 2) = A'B'C'. Hence our algorithm is O(n ^ 2). But it can be proved that, rather than, start from Pi + 1 as a second point and Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C), where s(Pi, Pi + 1, Pi + 2) = ABC.

Lets assume that we are wrong, and s(Pi + 1, Pi + 2, Pi + 3) = A'B'C' > s(pi, Pi + 1, Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

So, from all of the 5 possible cases, we found that, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 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 equalvector<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;}